Coverage for src/gwtransport/residence_time.py: 95%
328 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1"""
2Residence Time Calculations for Retarded Compound Transport.
4This module provides functions to compute residence times for compounds traveling through
5aquifer systems, accounting for flow variability, pore volume, and retardation due to
6physical or chemical interactions with the aquifer matrix. Residence time represents the
7duration a compound spends traveling from infiltration to extraction points, depending on
8flow rate (higher flow yields shorter residence time), pore volume (larger volume yields
9longer residence time), and retardation factor (interaction with matrix yields longer
10residence time).
12Available functions:
14- :func:`full` - Compute the flow-weighted mean residence time over
15 output bins, per pore volume (full ``(n_pore_volumes, n_bins)`` array). Follows the package's
16 bin-edge convention and is the form consumed elsewhere in the package. Supports both forward
17 (infiltration to extraction) and reverse (extraction to infiltration) directions.
19- :func:`mean` - Compute the mean residence time over output bins for a discrete
20 aquifer pore-volume distribution (an array of equally-weighted pore volumes). Collapses the
21 pore-volume axis to a single per-bin series. The ``spinup`` policy (default ``"constant"``)
22 warm-starts the spin-up by extrapolating the boundary flow.
24- :func:`gamma` - Compute the closed-form mean residence time over output bins for a
25 (shifted) gamma aquifer pore-volume distribution, with no pore-volume discretization. The
26 ``spinup`` policy (default ``"constant"``) warm-starts the spin-up; ``spinup=0.0`` instead
27 renormalizes over the covered sub-mass exactly.
29- :func:`fraction_explained_full`, :func:`fraction_explained_mean`,
30 :func:`fraction_explained_gamma` - Compute the **advective** fraction of each output bin that is
31 explained by the flow record: the flow-weighted share of the bin whose retarded advective parcel
32 was infiltrated/extracted inside the record. ``full`` returns one row per pore volume, ``mean``
33 the equal-weight discrete-APVD mean, and ``gamma`` the closed-form (shifted) gamma-APVD value,
34 mirroring :func:`full` / :func:`mean` / :func:`gamma`. These are **purely advective** -- molecular
35 diffusion and microdispersion spread each bin over a range of infiltration times that is
36 not captured here, so no bin is fully informed once dispersion is present (for that dispersive
37 informed fraction use the captured kernel mass of the diffusion coefficient matrix).
39- :func:`freundlich_retardation` - Compute concentration-dependent retardation factors from a
40 Freundlich isotherm, for use as the ``retardation_factor`` input to the transport functions.
42Spin-up period
43--------------
44The spin-up **region** is determined entirely by the supplied flow record (``tedges``, which
45fixes the cumulative throughflow volume ``V`` from ``0`` at the record start to ``V_end`` at the
46record end) together with the retarded pore volume ``retardation_factor * V_p`` -- it is not a
47length you set. A residence time for an output time needs the corresponding parcel to stay inside
48the flow record:
50* ``direction='extraction_to_infiltration'`` looks **back** to the infiltration event, so the
51 spin-up sits at the **start** of the output record: the residence time of a pore volume ``V_p``
52 needs ``V(t) >= retardation_factor * V_p`` (the extracted water was infiltrated before the record
53 began otherwise).
54* ``direction='infiltration_to_extraction'`` looks **forward** to the extraction event, so the
55 spin-up sits at the **end** of the output record: it needs ``V_end - V(t) >= retardation_factor *
56 V_p`` (the infiltrated water is extracted after the record ends otherwise).
58The spin-up therefore lengthens with both the pore volume and the retardation factor, and is
59longest for the largest pore volumes of a distribution.
61What happens in that region is governed by a ``spinup`` policy, following the package convention
62(see :mod:`gwtransport.advection`); :func:`full`, :func:`mean` and
63:func:`gamma` all share the contract ``spinup={'constant'} | None | float in
64[0, 1]`` and the default is ``"constant"`` everywhere:
66* ``"constant"`` (default) **warm-starts** by extrapolating the boundary flow (flow held constant
67 at its first/last value), so no in-record output is ``NaN``.
68* ``None`` is strict (no extrapolation), marking a pore volume ``NaN`` for any output bin its parcel
69 leaves the record within. Where the pore-volume axis is collapsed -- :func:`mean` over a
70 discrete set, :func:`gamma` over the continuum -- the bin mean then **renormalizes**
71 over the covered streamtubes / sub-mass, emitted wherever any coverage remains.
72* a ``float`` covered-fraction threshold is the strict mode with a minimum coverage gate: the
73 renormalized mean is emitted only where the covered streamtube fraction / sub-mass fraction is at
74 least ``spinup`` (``0.0`` matches ``None``; larger values demand more coverage). For the
75 per-pore-volume :func:`full` there is no axis to collapse, so the ``float`` behaves
76 exactly like ``None``.
78Output bins lying wholly outside ``tedges`` are ``NaN`` under every policy.
80The :func:`fraction_explained_full` / :func:`fraction_explained_mean` /
81:func:`fraction_explained_gamma` diagnostics report, per output bin, the advective fraction of the
82pore-volume distribution that is out of spin-up (``1.0`` = advectively fully informed, ``0.0`` =
83entirely in spin-up) and are the way to locate the spin-up region when the means warm-start over it.
85This file is part of gwtransport which is released under AGPL-3.0 license.
86See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
87"""
89import numpy as np
90import numpy.typing as npt
91import pandas as pd
92from scipy.stats import gamma as gamma_dist
94from gwtransport._time import tedges_to_days
95from gwtransport.gamma import parse_parameters
96from gwtransport.utils import cumulative_flow_volume, linear_interpolate
98# Relative slack on the covered-fraction spin-up gate. The covered sub-mass ``den`` equals its
99# fully-covered reference only up to summation-reassociation ulps, so an exact ``den >= threshold *
100# reference`` comparison spuriously rejects fully-covered bins at the strictest threshold (1.0). The
101# slack is far above that float noise (band widths reach ~1e4 pieces, ~1e-12 relative) yet negligible
102# as a physical coverage tolerance.
103_SPINUP_GATE_RELTOL = 1e-9
106def _boundary_extrapolated_map(
107 flow: npt.NDArray[np.floating],
108 flow_cum: npt.NDArray[np.floating],
109 tedges_days: npt.NDArray[np.floating],
110 pad: float,
111) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
112 """Extend the cumulative-volume -> time map past the record for the ``"constant"`` warm-start.
114 The boundary extrapolation slope is ``1/Q``, anchored on the nearest strictly-positive flow so a
115 zero-flow boundary bin (whose ``flow_cum`` step is only the strictly-monotone ulp bump) does not
116 give a ``1/0`` extrapolation. Bit-identical to ``1/Q_boundary`` when the boundary already carries
117 flow.
119 Returns
120 -------
121 volume_map : ndarray
122 ``flow_cum`` padded by ``pad`` on each end.
123 time_map : ndarray
124 ``tedges_days`` padded by ``pad * (1/Q)`` on each end, aligned with ``volume_map``.
125 Shared by :func:`full` and :func:`gamma`.
126 """
127 positive_flow = flow[flow > 0.0]
128 inv_q_first = 1.0 / positive_flow[0]
129 inv_q_last = 1.0 / positive_flow[-1]
130 volume_map = np.concatenate([[flow_cum[0] - pad], flow_cum, [flow_cum[-1] + pad]])
131 time_map = np.concatenate([
132 [tedges_days[0] - pad * inv_q_first],
133 tedges_days,
134 [tedges_days[-1] + pad * inv_q_last],
135 ])
136 return volume_map, time_map
139def _resolve_spinup(spinup: str | float | None) -> tuple[bool, float]:
140 """Normalize the residence-time ``spinup`` policy to ``(extrapolate, threshold)``.
142 The three sibling residence-time functions share one spin-up contract,
143 ``{'constant'} | None | float in [0, 1]``, matching the package convention (see
144 :mod:`gwtransport.advection`):
146 * ``'constant'`` -> ``(True, 0.0)`` -- warm-start by extrapolating the boundary flow.
147 * ``None`` -> ``(False, 0.0)`` -- strict map (no extrapolation); a collapsed mean is emitted
148 wherever any covered sub-mass / valid streamtube remains.
149 * ``float`` in ``[0, 1]`` -> ``(False, float)`` -- strict map, with the covered-fraction
150 threshold gating a collapsed mean. ``0.0`` matches ``None``; larger values require a larger
151 covered fraction before a bin is emitted, and ``1.0`` is strictest (full coverage required).
153 Parameters
154 ----------
155 spinup : {'constant'}, None, or float in [0, 1]
156 Public spin-up policy.
158 Returns
159 -------
160 extrapolate : bool
161 Whether to warm-start by extrapolating the boundary flow.
162 threshold : float
163 Covered-fraction gate applied where the pore-volume axis is collapsed (ignored by the
164 per-pore-volume :func:`full`, which is strict per bin).
166 Raises
167 ------
168 ValueError
169 If ``spinup`` is not ``'constant'``, ``None``, or a float in ``[0, 1]``.
170 """
171 if spinup == "constant":
172 return True, 0.0
173 if spinup is None:
174 return False, 0.0
175 if isinstance(spinup, int | float) and not isinstance(spinup, bool) and 0.0 <= spinup <= 1.0:
176 return False, float(spinup)
177 msg = "spinup should be 'constant', None, or a float in [0, 1]"
178 raise ValueError(msg)
181def _phi_setup(
182 flow: npt.NDArray[np.floating],
183 flow_cum: npt.NDArray[np.floating],
184 tedges_days: npt.NDArray[np.floating],
185 *,
186 extrapolate: bool,
187 pad: float,
188) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating], npt.NDArray[np.floating], npt.NDArray[np.floating]]:
189 """Build the antiderivative ``phi(x) = int_0^x T(w) dw`` of the cumulative-volume -> time map ``T``.
191 ``T`` is piecewise-linear (knots at the cumulative-volume edges), so ``phi`` is piecewise-quadratic
192 with the same knots. With ``extrapolate`` the map is extended one anchor past each end of the record
193 at the boundary flow rate (padded by ``pad``) for the ``"constant"`` warm-start; otherwise it is the
194 raw record. Shared by :func:`full` and :func:`gamma`.
196 Returns
197 -------
198 phi_v : ndarray
199 Cumulative-volume knots of the map (extended by ``pad`` at each end when ``extrapolate``).
200 phi_t : ndarray
201 Time knots aligned with ``phi_v``.
202 phi_knot : ndarray
203 ``phi`` evaluated at each volume knot.
204 phi_rate : ndarray
205 Per-segment ``dV/dt`` (flow) rate between consecutive knots.
206 """
207 # pad == 0 (no spin-up reach, e.g. retardation_factor or all pore volumes 0) carries no
208 # extrapolation and would only add zero-width boundary segments (0/0 -> NaN phi_rate), so skip it.
209 if extrapolate and pad > 0.0 and np.any(flow > 0.0):
210 phi_v, phi_t = _boundary_extrapolated_map(flow, flow_cum, tedges_days, pad)
211 else:
212 phi_v, phi_t = flow_cum, tedges_days
213 phi_dv = phi_v[1:] - phi_v[:-1]
214 phi_rate = phi_dv / (phi_t[1:] - phi_t[:-1])
215 phi_knot = np.concatenate([[0.0], np.cumsum(phi_t[:-1] * phi_dv + phi_dv**2 / (2 * phi_rate))])
216 return phi_v, phi_t, phi_knot, phi_rate
219def _eval_phi(
220 x: npt.NDArray[np.floating],
221 phi_v: npt.NDArray[np.floating],
222 phi_t: npt.NDArray[np.floating],
223 phi_knot: npt.NDArray[np.floating],
224 phi_rate: npt.NDArray[np.floating],
225 *,
226 strict_nan: bool = False,
227) -> npt.NDArray[np.floating]:
228 """Evaluate the piecewise-quadratic antiderivative ``phi`` from :func:`_phi_setup` at ``x``.
230 ``x`` is clipped to the map range ``[phi_v[0], phi_v[-1]]`` (the warm-start extrapolation lives in
231 that range when padded). With ``strict_nan`` any ``x`` outside the range returns ``NaN`` instead --
232 used by :func:`full` so an output bin whose parcel leaves the record is ``NaN``.
234 Returns
235 -------
236 ndarray
237 ``phi`` evaluated at ``x``, with the same shape as ``x``.
238 """
239 x = np.asarray(x, dtype=float)
240 xc = np.clip(x, phi_v[0], phi_v[-1])
241 j = np.clip(np.searchsorted(phi_v, xc, side="right") - 1, 0, len(phi_rate) - 1)
242 dv = xc - phi_v[j]
243 out = phi_knot[j] + phi_t[j] * dv + dv * dv / (2 * phi_rate[j])
244 if strict_nan:
245 out = np.where((x < phi_v[0]) | (x > phi_v[-1]), np.nan, out)
246 return out
249def full(
250 *,
251 flow: npt.ArrayLike,
252 tedges: pd.DatetimeIndex | np.ndarray,
253 cout_tedges: pd.DatetimeIndex | np.ndarray,
254 aquifer_pore_volumes: npt.ArrayLike,
255 direction: str = "extraction_to_infiltration",
256 retardation_factor: float = 1.0,
257 spinup: str | float | None = "constant",
258) -> npt.NDArray[np.floating]:
259 r"""
260 Compute the mean residence time over output bins, per pore volume.
262 The flow-weighted mean residence time is computed over each output interval
263 ``[cout_tedges[i], cout_tedges[i + 1])`` and returned as the full
264 ``(n_pore_volumes, n_output_bins)`` array -- one row per entry in
265 ``aquifer_pore_volumes``, without collapsing the pore-volume axis. The average is uniform in
266 cumulative throughflow volume, matching the package's bin-edge convention (and what the diffusion
267 modules consume to compute a per-bin retarded velocity).
269 Parameters
270 ----------
271 flow : array-like
272 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one.
273 tedges : pandas.DatetimeIndex
274 Time edges for the flow data, as datetime64 objects, defining the flow intervals.
275 cout_tedges : pandas.DatetimeIndex
276 Output time edges as datetime64 objects; ``n + 1`` edges define ``n`` output bins.
277 aquifer_pore_volumes : float or array-like
278 Pore volume(s) of the aquifer [m³]. A single value or an array of pore volumes
279 representing different flow paths.
280 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
281 Direction of the flow calculation:
283 * 'extraction_to_infiltration':
284 Extraction to infiltration modeling - how many days ago was the extracted water infiltrated.
285 * 'infiltration_to_extraction':
286 Infiltration to extraction modeling - how many days until the infiltrated water is extracted.
288 Default is 'extraction_to_infiltration'.
289 retardation_factor : float, optional
290 Retardation factor of the compound in the aquifer [dimensionless]. A value greater
291 than 1.0 indicates the compound moves slower than water. Default is 1.0.
292 spinup : {'constant'}, None, or float in [0, 1], optional
293 How to treat the spin-up zone, where a pore volume's retarded look-back/forward parcel
294 leaves the flow record. Matches the package convention (see :mod:`gwtransport.advection`).
296 * ``'constant'`` (default): warm-start -- extrapolate the cumulative-volume-to-time map past
297 the record at the boundary flow rates (flow held constant at its first/last value), so
298 the residence time stays finite. No left-edge (extraction) or right-edge (infiltration)
299 spin-up ``NaN``.
300 * ``None`` or a ``float`` in ``[0, 1]``: strict -- a pore volume whose parcel leaves the
301 record at any point within an output bin is ``NaN`` for that bin (all-or-nothing per bin),
302 with no extrapolation. This function returns the full per-pore-volume array, so there is no
303 pore-volume axis to collapse; the ``float`` covered-fraction threshold therefore behaves
304 identically to ``None`` here and only takes effect once the axis is collapsed in
305 :func:`mean` / :func:`gamma`.
307 Output bins lying wholly outside ``tedges`` are ``NaN`` under either policy.
309 Returns
310 -------
311 numpy.ndarray
312 Mean residence time [days], shape ``(n_pore_volumes, n_output_bins)``. The first
313 dimension corresponds to the pore volumes and the second to the ``cout_tedges`` bins.
314 Negative or ``NaN`` ``flow`` makes the cumulative-volume map non-monotone or undefined; the
315 whole array is returned as ``NaN`` (the function refuses rather than raising).
317 Raises
318 ------
319 ValueError
320 If ``tedges`` does not have exactly one more element than ``flow``. If
321 ``direction`` is not ``'extraction_to_infiltration'`` or
322 ``'infiltration_to_extraction'``. If ``spinup`` is not ``'constant'``, ``None``, or a float
323 in ``[0, 1]``.
325 See Also
326 --------
327 fraction_explained_full : Advective fraction of each output bin explained, per pore volume
328 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
329 :ref:`concept-transport-equation` : Flow-weighted averaging convention
331 Notes
332 -----
333 With the default ``spinup='constant'`` the spin-up zone is warm-started by extrapolating the
334 boundary flow, so no in-record bin is ``NaN``; use :func:`fraction_explained_mean` (or
335 ``spinup=None``) to locate the spin-up region. See the module docstring (``Spin-up period``)
336 for the full rule.
338 The single-streamtube residence time :math:`\tau(V) = \mathrm{sign}\,[T(V + \mathrm{sign}\,R V_p)
339 - T(V)]` is piecewise-linear in cumulative throughflow volume :math:`V` (:math:`T` is the
340 volume :math:`\to` time map, :math:`\mathrm{sign} = -1` for ``extraction_to_infiltration`` and
341 :math:`+1` for ``infiltration_to_extraction``). Its flow-weighted bin average is therefore a
342 closed-form difference of the antiderivative :math:`\Phi(x) = \int_0^x T(w)\,dw` (piecewise-
343 quadratic), evaluated at four points per pore volume and output bin:
345 .. math::
347 \bar\tau
348 = \frac{1}{\Delta V}\int_{V_\mathrm{lo}}^{V_\mathrm{hi}} \tau(V)\,dV
349 = \frac{\mathrm{sign}}{\Delta V}\bigl[
350 \Phi(V_\mathrm{hi} + \mathrm{sign}\,R V_p) - \Phi(V_\mathrm{lo} + \mathrm{sign}\,R V_p)
351 - \Phi(V_\mathrm{hi}) + \Phi(V_\mathrm{lo})\bigr],
353 where :math:`V` is cumulative throughflow volume (:math:`dV = Q\,dt`). This avoids materialising a
354 per-streamtube integration grid, so memory and time scale as the output size
355 :math:`O(n_\mathrm{pore\ volumes}\cdot n_\mathrm{bins})`. A zero-throughflow output bin
356 (:math:`\Delta V \to 0`) has a fixed volume while output time advances, so it degenerates to the
357 pointwise residence time at the bin's time midpoint.
359 Examples
360 --------
361 >>> import pandas as pd
362 >>> import numpy as np
363 >>> from gwtransport.residence_time import full
364 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
365 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day
366 >>> mean_times = full(
367 ... flow=flow_values,
368 ... tedges=flow_dates,
369 ... cout_tedges=flow_dates,
370 ... aquifer_pore_volumes=200.0,
371 ... direction="extraction_to_infiltration",
372 ... )
373 >>> # 200 m³ / 100 m³/day = 2 days residence time; the default constant warm-start
374 >>> # extrapolates the boundary flow, so the left-edge spin-up bins are also 2 days
375 >>> print(mean_times) # doctest: +NORMALIZE_WHITESPACE
376 [[2. 2. 2. 2. 2. 2. 2. 2. 2.]]
377 """
378 if direction not in {"extraction_to_infiltration", "infiltration_to_extraction"}:
379 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
380 raise ValueError(msg)
381 extrapolate, _ = _resolve_spinup(spinup)
383 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes)
384 tedges = pd.DatetimeIndex(tedges)
385 cout_tedges = pd.DatetimeIndex(cout_tedges)
386 flow = np.asarray(flow, dtype=float)
387 n_pv = len(aquifer_pore_volumes)
388 n_out = len(cout_tedges) - 1
390 if len(tedges) != len(flow) + 1:
391 msg = "tedges must have one more element than flow"
392 raise ValueError(msg)
393 if np.any(flow < 0) or np.any(np.isnan(flow)):
394 return np.full((n_pv, n_out), np.nan)
396 tedges_days = tedges_to_days(tedges)
397 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0])
398 # Plateaus in flow_cum from Q = 0 bins make the V -> t inversion multi-valued; bump duplicates by
399 # the smallest representable amount so the inverse map is single-valued.
400 flow_cum = cumulative_flow_volume(flow, np.diff(tedges_days), strictly_monotone=True)
402 # Sign convention: sign = -1 for extraction_to_infiltration, +1 for infiltration_to_extraction;
403 # the look-back/forward parcel sits at volume V + shift and tau(V) = sign * (T(V + shift) - T(V)) is
404 # piecewise-linear in V (T is the volume -> time map). Its flow-weighted bin average is a closed-
405 # form difference of the antiderivative phi(x) = int_0^x T(w) dw, so no per-streamtube integration
406 # grid is built (memory/time O(n_pore_volumes * n_bins), not O(n_pore_volumes^2 * n_flow)).
407 sign = -1.0 if direction == "extraction_to_infiltration" else 1.0
408 shift = sign * retardation_factor * aquifer_pore_volumes # (n_pv,)
410 # phi over the cumulative-volume -> time map. With spinup="constant" the map is extrapolated past
411 # the record at the boundary flow (padded by the largest reach R * max(V_p)) so phi warm-starts the
412 # spin-up; otherwise phi is NaN outside the record so a bin whose parcel leaves it becomes NaN.
413 pad = retardation_factor * float(aquifer_pore_volumes.max()) if aquifer_pore_volumes.size else 0.0
414 phi_v, phi_t, phi_knot, phi_rate = _phi_setup(flow, flow_cum, tedges_days, extrapolate=extrapolate, pad=pad)
415 # The map is only actually extended when there is a positive boundary flow to extrapolate; with
416 # all-zero flow (or spinup=None) it stays the raw record, so out-of-record look-backs are NaN.
417 strict = not (extrapolate and bool(np.any(flow > 0.0)))
419 vol_out = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days, left=np.nan, right=np.nan)
420 v_lo = vol_out[:-1]
421 v_hi = vol_out[1:]
422 dvol = v_hi - v_lo
423 bins_within = np.isfinite(v_lo) & np.isfinite(v_hi)
425 phi_base = _eval_phi(vol_out, phi_v, phi_t, phi_knot, phi_rate, strict_nan=strict) # (n_out + 1,)
426 phi_shift = _eval_phi(vol_out[None, :] + shift[:, None], phi_v, phi_t, phi_knot, phi_rate, strict_nan=strict)
428 tol = 1e6 * np.spacing(float(np.max(np.abs(flow_cum)))) if flow_cum.size else 0.0
429 with np.errstate(divide="ignore", invalid="ignore"):
430 result = (
431 sign * ((phi_shift[:, 1:] - phi_shift[:, :-1]) - (phi_base[1:] - phi_base[:-1])[None, :]) / dvol[None, :]
432 )
433 # Zero-throughflow output bin (dvol -> 0): the volume is fixed while output time advances, so the
434 # flow-weighted average degenerates to the pointwise tau at the bin time midpoint, sign*(T(V_lo +
435 # shift) - t_mid). A direct ratio there would catastrophically cancel. Only in-record zero-flow
436 # bins reach it, so skip the interp entirely when none are present.
437 degenerate = bins_within & (dvol <= tol)
438 if np.any(degenerate):
439 t_mid = 0.5 * (cout_tedges_days[:-1] + cout_tedges_days[1:])
440 nan_outside = np.nan if strict else None
441 with np.errstate(divide="ignore", invalid="ignore"):
442 t_lookback = linear_interpolate(
443 x_ref=phi_v, y_ref=phi_t, x_query=v_lo[None, :] + shift[:, None], left=nan_outside, right=nan_outside
444 )
445 point = sign * (t_lookback - t_mid[None, :])
446 result = np.where(dvol[None, :] > tol, result, point)
447 return np.where(bins_within[None, :], result, np.nan)
450def mean(
451 *,
452 flow: npt.ArrayLike,
453 tedges: pd.DatetimeIndex | np.ndarray,
454 cout_tedges: pd.DatetimeIndex | np.ndarray,
455 aquifer_pore_volumes: npt.ArrayLike,
456 direction: str = "extraction_to_infiltration",
457 retardation_factor: float = 1.0,
458 spinup: str | float | None = "constant",
459) -> npt.NDArray[np.floating]:
460 r"""
461 Compute the mean residence time over output bins for a discrete APVD.
463 The mean is taken over a **discrete** set of equally-weighted aquifer pore volumes -- one
464 streamtube per entry in ``aquifer_pore_volumes``. Each streamtube's flow-weighted bin average
465 is computed with :func:`full` and the pore-volume axis is then collapsed to a
466 single per-output-bin series by averaging over the streamtubes that are valid in each bin. For
467 a continuous (shifted) gamma pore-volume distribution evaluated in closed form, use
468 :func:`gamma`.
470 The mean is over the valid streamtubes,
472 .. math::
474 \bar\tau_b = \frac{1}{|V_b|}\sum_{i \in V_b} \tau_{i,b},
475 \qquad V_b = \{\, i : \tau_{i,b}\ \mathrm{finite} \,\}.
477 With the default ``spinup='constant'`` every streamtube is finite within the flow record
478 (the boundary flow is extrapolated), so this is simply the mean over all pore volumes; with
479 ``spinup=None`` it renormalizes over the streamtubes that have broken through.
481 Parameters
482 ----------
483 flow : array-like
484 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one.
485 tedges : pandas.DatetimeIndex
486 Time edges for the flow data, as datetime64 objects, defining the flow intervals.
487 cout_tedges : pandas.DatetimeIndex
488 Output time edges as datetime64 objects; ``n + 1`` edges define ``n`` output bins.
489 aquifer_pore_volumes : array-like
490 Discrete pore volumes [m³], one per (equally-weighted) streamtube. A single value
491 collapses to the per-streamtube mean of :func:`full`.
492 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
493 Direction of the flow calculation:
494 * 'extraction_to_infiltration': how many days ago was the extracted water infiltrated
495 * 'infiltration_to_extraction': how many days until the infiltrated water is extracted
496 Default is 'extraction_to_infiltration'.
497 retardation_factor : float, optional
498 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
499 spinup : {'constant'}, None, or float in [0, 1], optional
500 Spin-up policy, sharing the contract of :func:`gamma`. ``'constant'``
501 (default) warm-starts by extrapolating the boundary flow so no in-record bin is ``NaN``;
502 ``None`` leaves spin-up streamtubes ``NaN`` and the mean renormalizes over those that have
503 broken through (emitted wherever at least one streamtube is valid). A ``float`` in
504 ``[0, 1]`` is the covered-fraction threshold: the renormalized mean is emitted only where
505 the fraction of valid streamtubes is at least ``spinup`` (``0.0`` matches ``None``; ``1.0``
506 demands every streamtube; larger values demand more streamtubes to have broken through). Use
507 :func:`fraction_explained_mean` to
508 locate the spin-up region.
510 Returns
511 -------
512 numpy.ndarray
513 Mean residence time [days], shape ``(n_output_bins,)``. Output bins with no valid
514 streamtube (outside the flow record, or -- with ``spinup=None`` -- fully in the spin-up
515 zone) are NaN; with a ``float`` ``spinup`` so are bins whose valid-streamtube fraction is
516 below the threshold. Negative or ``NaN`` ``flow`` makes the cumulative-volume map non-monotone
517 or undefined; the whole series is returned as ``NaN`` (the function refuses rather than raising).
519 See Also
520 --------
521 gamma : Exact closed-form mean for a continuous (shifted) gamma APVD
522 full : Per-pore-volume mean residence time over output bins
523 fraction_explained_mean : Advective fraction of each output bin explained by the record
524 gwtransport.gamma.bins : Discretize a gamma APVD into pore-volume bins
525 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
527 Notes
528 -----
529 With ``spinup=None`` the spin-up is **all-or-nothing per streamtube**: a streamtube whose
530 look-back/forward parcel leaves the flow record part-way through an output bin has a ``NaN`` bin
531 average (inherited from :func:`full`) and is dropped from that bin's mean
532 entirely, rather than contributing its partially-covered share; the bin is ``NaN`` only once
533 every streamtube is in spin-up. In that mode the discrete mean differs from
534 :func:`gamma`, which renormalizes over the covered sub-mass exactly. See the
535 module docstring (``Spin-up period``) for the full rule.
537 Examples
538 --------
539 >>> import pandas as pd
540 >>> import numpy as np
541 >>> from gwtransport.residence_time import mean
542 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-02-10", freq="D")
543 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # 100 m³/day
544 >>> tau_bar = mean(
545 ... flow=flow_values,
546 ... tedges=flow_dates,
547 ... cout_tedges=flow_dates,
548 ... aquifer_pore_volumes=[400.0, 600.0], # two equally-weighted streamtubes
549 ... )
550 >>> # Deep in the record: mean pore volume 500 / 100 m³/day = 5 days
551 >>> float(np.round(tau_bar[-1], 6))
552 5.0
553 """
554 _, threshold = _resolve_spinup(spinup)
555 rt = full(
556 flow=flow,
557 tedges=tedges,
558 cout_tedges=cout_tedges,
559 aquifer_pore_volumes=aquifer_pore_volumes,
560 direction=direction,
561 retardation_factor=retardation_factor,
562 spinup=spinup,
563 )
565 # Mean over the streamtubes that are valid (non-NaN) in each output bin; bins with no valid
566 # streamtube reduce to 0/0 and are NaN. With spinup='constant' every in-record streamtube is
567 # finite, so this is the plain mean; otherwise it renormalizes over the broken-through set and
568 # a float covered-fraction threshold further NaNs bins where too few streamtubes have arrived.
569 n_streamtubes = rt.shape[0]
570 valid_count = np.isfinite(rt).sum(axis=0)
571 with np.errstate(invalid="ignore"):
572 mean = np.nansum(rt, axis=0) / valid_count
573 if threshold > 0.0:
574 mean = np.where(valid_count >= threshold * n_streamtubes, mean, np.nan)
575 return mean
578def gamma(
579 *,
580 flow: npt.ArrayLike,
581 tedges: pd.DatetimeIndex | np.ndarray,
582 cout_tedges: pd.DatetimeIndex | np.ndarray,
583 mean: float | None = None,
584 std: float | None = None,
585 loc: float = 0.0,
586 alpha: float | None = None,
587 beta: float | None = None,
588 direction: str = "extraction_to_infiltration",
589 retardation_factor: float = 1.0,
590 spinup: str | float | None = "constant",
591 _max_tile_elements: int = 1_000_000,
592) -> npt.NDArray[np.floating]:
593 r"""
594 Compute the mean residence time over output bins for a (shifted) gamma APVD.
596 The expectation over a (shifted) gamma aquifer pore-volume distribution (APVD),
597 parameterized by either ``(mean, std, loc)`` or ``(alpha, beta, loc)``, is taken in closed
598 form -- no pore-volume binning and no ``n_bins`` accuracy/cost knob. The bin mean is
599 flow-weighted (uniform in cumulative volume), matching the bin-edge convention of the
600 package, and a single per-output-bin series is returned.
602 The single-streamtube residence time is piecewise-linear in the pore volume :math:`V_p`, so
603 its per-bin time integral :math:`G_b(V_p) = \int_{\mathrm{bin}} \tau\,dV` is piecewise-
604 quadratic in :math:`V_p` and the covered length :math:`L_b(V_p)` piecewise-linear. The bin
605 mean is the ratio of two closed-form integrals against the gamma density -- its zeroth,
606 first and second partial moments (regularized incomplete gamma) -- formed once after
607 integrating. The ``spinup`` policy sets what happens where part of the APVD lacks flow
608 history: ``'constant'`` (default) extrapolates the boundary flow over the full distribution
609 (the package default warm-start), while a ``float`` threshold renormalizes the mean over the
610 covered sub-mass (``0.0`` reproduces the exact covered-sub-mass conditional mean).
612 Parameters
613 ----------
614 flow : array-like
615 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one.
616 tedges : pandas.DatetimeIndex
617 Time edges for the flow data, as datetime64 objects, defining the flow intervals.
618 cout_tedges : pandas.DatetimeIndex
619 Output time edges as datetime64 objects; ``n + 1`` edges define ``n`` output bins.
620 mean : float, optional
621 Mean of the gamma APVD [m³]. Must be strictly greater than ``loc``. Provide either
622 ``(mean, std)`` or ``(alpha, beta)``.
623 std : float, optional
624 Standard deviation of the gamma APVD [m³]. Must be positive.
625 loc : float, optional
626 Location (lower bound of support) of the gamma APVD [m³]; a guaranteed minimum pore
627 volume. Must satisfy ``0 <= loc < mean``. Default is 0.0.
628 alpha : float, optional
629 Shape parameter of the gamma APVD (must be > 0).
630 beta : float, optional
631 Scale parameter of the gamma APVD (must be > 0).
632 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
633 Direction of the flow calculation:
634 * 'extraction_to_infiltration': how many days ago was the extracted water infiltrated
635 * 'infiltration_to_extraction': how many days until the infiltrated water is extracted
636 Default is 'extraction_to_infiltration'.
637 retardation_factor : float, optional
638 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
639 spinup : {'constant'}, None, or float in [0, 1], optional
640 How to treat the spin-up zone, where part of the gamma APVD lacks flow history. Matches
641 the package convention (see :mod:`gwtransport.advection`).
643 * ``'constant'`` (default): warm-start -- extrapolate the cumulative-volume-to-time map past
644 the record at the boundary flow rates (flow held constant at its first/last value) and
645 integrate the full distribution, so no in-record bin is ``NaN``.
646 * ``None`` or a ``float`` in ``[0, 1]``: renormalize the mean over the covered sub-mass,
647 emitting a bin only where the covered fraction of the distribution is at least the
648 threshold. ``None`` and ``0.0`` both give the exact covered-sub-mass conditional mean
649 (emit whenever any sub-mass is covered); larger values demand a larger covered fraction,
650 and ``1.0`` requires the full distribution to be covered.
652 Output bins lying wholly outside ``tedges`` are ``NaN`` under either policy.
654 Returns
655 -------
656 numpy.ndarray
657 APVD-mean residence time [days], shape ``(n_output_bins,)``. Output bins outside the flow
658 record are NaN; with a ``float`` ``spinup`` so are bins whose covered fraction is below the
659 threshold. Negative or ``NaN`` ``flow`` makes the cumulative-volume map non-monotone or
660 undefined; the whole series is returned as ``NaN`` (the function refuses rather than raising).
662 Raises
663 ------
664 ValueError
665 If ``tedges`` does not have exactly one more element than ``flow``. If ``direction``
666 is not ``'extraction_to_infiltration'`` or ``'infiltration_to_extraction'``. If ``spinup``
667 is not ``'constant'``, ``None``, or a float in ``[0, 1]``. Gamma parameter validation is
668 delegated to :func:`gwtransport.gamma.parse_parameters`.
670 See Also
671 --------
672 mean : Equally-weighted mean for a discrete set of pore volumes
673 full : Per-pore-volume mean residence time over output bins
674 fraction_explained_mean : Advective fraction of each output bin explained by the record
675 gwtransport.gamma.bins : Discretize a gamma APVD into pore-volume bins
676 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
677 :ref:`concept-gamma-distribution` : Two-parameter pore volume model
679 Notes
680 -----
681 With the default ``spinup='constant'`` the spin-up is warm-started exactly as in
682 :func:`mean` (constant-boundary-flow extrapolation), so the two agree everywhere. With
683 ``spinup=0.0`` the spin-up is instead handled by exact covered-sub-mass renormalization: each
684 output bin integrates over only the pore-volume sub-range with sufficient flow history. See the
685 module docstring (``Spin-up period``) for the full rule.
687 Examples
688 --------
689 >>> import pandas as pd
690 >>> import numpy as np
691 >>> from gwtransport.residence_time import gamma
692 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-02-10", freq="D")
693 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # 100 m³/day
694 >>> tau_bar = gamma(
695 ... flow=flow_values,
696 ... tedges=flow_dates,
697 ... cout_tedges=flow_dates,
698 ... mean=500.0,
699 ... std=100.0,
700 ... direction="extraction_to_infiltration",
701 ... )
702 >>> # Deep in the record the mean residence time approaches mean / flow = 5 days
703 >>> float(np.round(tau_bar[-1], 6))
704 5.0
705 """
706 if direction not in {"extraction_to_infiltration", "infiltration_to_extraction"}:
707 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
708 raise ValueError(msg)
709 extrapolate, spinup_threshold = _resolve_spinup(spinup)
711 alpha, beta, loc = parse_parameters(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta)
713 flow = np.asarray(flow, dtype=float)
714 tedges = pd.DatetimeIndex(tedges)
715 cout_tedges = pd.DatetimeIndex(cout_tedges)
716 n_out = len(cout_tedges) - 1
718 if len(tedges) != len(flow) + 1:
719 msg = "tedges must have one more element than flow"
720 raise ValueError(msg)
721 if np.any(flow < 0) or np.any(np.isnan(flow)):
722 return np.full(n_out, np.nan)
724 sign = -1.0 if direction == "extraction_to_infiltration" else 1.0
725 r = retardation_factor
727 tedges_days = tedges_to_days(tedges)
728 flow_cum = cumulative_flow_volume(flow, np.diff(tedges_days), strictly_monotone=True)
729 v_end = flow_cum[-1]
730 n_edges = len(flow_cum)
732 # Finite support: drop the gamma tails (mass ~1e-13, far below the discretization error this
733 # closed form replaces). Restricting the integral to [support_lo, support_hi] is what keeps
734 # the per-bin flow-edge band -- and the gamma CDF evaluations -- bounded.
735 tail = 1e-13
736 support_lo = float(loc + gamma_dist.ppf(tail, alpha, scale=beta))
737 support_hi = float(loc + gamma_dist.ppf(1.0 - tail, alpha, scale=beta))
739 # phi(v) = int_0^v T(w) dw with T the inverse cumulative-volume map (piecewise-linear); phi is
740 # piecewise-quadratic with knots at the cumulative-volume edges, so the per-bin time integral of
741 # tau is a difference of phi at the look-back/forward limits. With spinup="constant" the map is
742 # extended past the record at the boundary flow rates (one anchor each end, padded by the largest
743 # reach r*support_hi) so phi extrapolates the spin-up; with a float spinup it stays clipped to
744 # [0, v_end] (the covered sub-mass only).
745 phi_v, phi_t, phi_knot, phi_rate = _phi_setup(
746 flow, flow_cum, tedges_days, extrapolate=extrapolate, pad=r * support_hi
747 )
749 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0])
750 vol_out = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days, left=np.nan, right=np.nan)
751 good_all = np.isfinite(vol_out[:-1]) & np.isfinite(vol_out[1:]) & (vol_out[1:] > vol_out[:-1])
752 if not np.any(good_all):
753 return np.full(n_out, np.nan)
754 v_lo_all = np.where(good_all, vol_out[:-1], 0.0)
755 v_hi_all = np.where(good_all, vol_out[1:], 1.0)
757 # Fixed band of flow edges each bin's look-back/forward sweep can cross over the supported
758 # pore volumes. band_width is the global maximum so every tile shares one column layout; a
759 # spurious column clips to the support and merely splits a quadratic piece without changing
760 # its integral, so the band is an exact superset.
761 a_min = v_lo_all - r * support_hi if sign < 0 else v_lo_all + r * support_lo
762 a_max = v_hi_all - r * support_lo if sign < 0 else v_hi_all + r * support_hi
763 jlo_all = np.clip(np.searchsorted(flow_cum, a_min, "left") - 1, 0, n_edges - 1)
764 jhi_all = np.clip(np.searchsorted(flow_cum, a_max, "right"), 0, n_edges - 1)
765 band_width = int((jhi_all - jlo_all).max()) + 1
766 band = np.arange(band_width)
768 # Tile over output bins to bound peak memory. Each bin carries ~3*band_width pieces through a
769 # (3 nodes, 3 phi-args) stack of 9 * n_pieces elements; size the tile to that element budget.
770 n_pieces = 3 * band_width + 3
771 tile = max(1, _max_tile_elements // (9 * n_pieces))
773 result = np.full(n_out, np.nan)
774 for t0 in range(0, n_out, tile):
775 t1 = min(t0 + tile, n_out)
776 nt = t1 - t0
777 good = good_all[t0:t1]
778 v_lo = v_lo_all[t0:t1]
779 v_hi = v_hi_all[t0:t1]
780 cols = np.clip(jlo_all[t0:t1, None] + band[None, :], 0, n_edges - 1)
781 fcb = flow_cum[cols]
782 vlo = v_lo[:, None]
783 vhi = v_hi[:, None]
785 # Direction-specific breakpoint families: the V_p values where a phi argument (a clipped
786 # look-back/forward limit) crosses a flow edge or a validity bound. Clip to the support
787 # and sort to form the integration pieces (zero-width pieces contribute nothing).
788 if sign < 0:
789 cand = np.concatenate([(vhi - fcb) / r, (vlo - fcb) / r, fcb / r, vlo / r, vhi / r], axis=1)
790 else:
791 cand = np.concatenate(
792 [(fcb - vlo) / r, (fcb - vhi) / r, (v_end - fcb) / r, (v_end - vlo) / r, (v_end - vhi) / r], axis=1
793 )
794 np.clip(cand, support_lo, support_hi, out=cand)
795 cand.sort(axis=1)
796 edges = np.concatenate([np.full((nt, 1), support_lo), cand, np.full((nt, 1), support_hi)], axis=1)
797 lo = edges[:, :-1]
798 hi = edges[:, 1:]
799 mid = 0.5 * (lo + hi)
800 half = 0.5 * (hi - lo)
801 nodes = np.stack([lo, mid, hi], axis=0) # (3, nt, n_pieces)
803 # G(V_p) at the three quadrature nodes (piece lo/mid/hi) as a difference of phi. The constant
804 # term phi(v_hi)/phi(v_lo) does not depend on V_p, so evaluate it once per bin and batch only
805 # the three V_p-dependent phi arguments. With spinup="constant" the full bin is integrated
806 # against the extrapolated phi (length is the full bin width); with a float spinup the limits
807 # clamp to the streamtube's covered sub-interval and the covered length renormalizes.
808 if sign < 0:
809 a_hi = vhi[None] - r * nodes
810 if extrapolate:
811 v_start = np.broadcast_to(vlo[None], a_hi.shape)
812 a_lo = vlo[None] - r * nodes
813 length = np.broadcast_to(vhi[None] - vlo[None], a_hi.shape)
814 else:
815 v_start = np.maximum(vlo[None], r * nodes)
816 a_lo = np.maximum(vlo[None] - r * nodes, 0.0)
817 length = np.maximum(vhi[None] - v_start, 0.0)
818 phi_const = _eval_phi(v_hi, phi_v, phi_t, phi_knot, phi_rate)
819 phi_stack = _eval_phi(np.stack([v_start, a_hi, a_lo]), phi_v, phi_t, phi_knot, phi_rate)
820 g = phi_const[None, :, None] - phi_stack[0] - phi_stack[1] + phi_stack[2]
821 else:
822 a_lo = vlo[None] + r * nodes
823 if extrapolate:
824 v_stop = np.broadcast_to(vhi[None], a_lo.shape)
825 a_hi = vhi[None] + r * nodes
826 length = np.broadcast_to(vhi[None] - vlo[None], a_lo.shape)
827 else:
828 v_stop = np.minimum(vhi[None], v_end - r * nodes)
829 a_hi = np.minimum(vhi[None] + r * nodes, v_end)
830 length = np.maximum(v_stop - vlo[None], 0.0)
831 phi_const = _eval_phi(v_lo, phi_v, phi_t, phi_knot, phi_rate)
832 phi_stack = _eval_phi(np.stack([a_hi, a_lo, v_stop]), phi_v, phi_t, phi_knot, phi_rate)
833 g = phi_stack[0] - phi_stack[1] - phi_stack[2] + phi_const[None, :, None]
834 g = np.where(length > 0, g, 0.0)
835 length = np.where(length > 0, length, 0.0)
836 g_lo, g_mid, g_hi = g
837 l_lo, l_mid, l_hi = length
839 # Gamma partial moments over each piece: one CDF per shape on the piece edges, then diff.
840 # M1/M2 follow from the shifted-gamma partial-moment identities.
841 cdf_edges = edges - loc
842 f0 = gamma_dist.cdf(cdf_edges, alpha, scale=beta)
843 f1 = gamma_dist.cdf(cdf_edges, alpha + 1, scale=beta)
844 f2 = gamma_dist.cdf(cdf_edges, alpha + 2, scale=beta)
845 m0 = np.diff(f0, axis=1)
846 d1 = np.diff(f1, axis=1)
847 m1 = alpha * beta * d1 + loc * m0
848 d2 = np.diff(f2, axis=1)
849 m2 = alpha * (alpha + 1) * beta**2 * d2 + 2 * loc * alpha * beta * d1 + loc * loc * m0
851 # Piece-centred gamma moments mu_k = int (x - mid)^k f over each piece. They are formed by
852 # shifting the raw partial moments (m1 - mid m0, m2 - 2 mid m1 + mid^2 m0), which for R > 1
853 # subtracts large near-equal terms (|mid| up to support_hi, with m1/m2 carrying the matching
854 # mid powers) and catastrophically cancels on near-degenerate pieces. They are, however,
855 # bounded purely by the piece geometry -- |x - mid| <= half over the piece -- so clip each to
856 # its exact range. This is not a strict no-op on well-conditioned pieces (a raw value may sit
857 # a rounding step outside the tight bound; clamping it back is a negligible correction). It
858 # earns its place on the near-degenerate pieces, where the cancellation noise -- which
859 # otherwise pairs with the 1/half^2 second-difference below to manufacture a spurious
860 # contribution -- is forced back into the physically valid band.
861 b0 = np.maximum(m0, 0.0)
862 mu1 = np.clip(m1 - mid * m0, -half * b0, half * b0)
863 mu2 = np.clip(m2 - 2 * mid * m1 + mid**2 * m0, 0.0, half**2 * b0)
865 # G is quadratic per piece, L linear; reconstruct each from its three nodes (Lagrange-
866 # exact) centred at mid and contract with the moments: int (a(x-mid)^2 + b(x-mid) + c) f
867 # = a mu2 + b mu1 + c mu0. Zero-width pieces (duplicate breakpoints) divide by half == 0;
868 # their masked result is discarded by np.where.
869 safe = half > 0
870 with np.errstate(divide="ignore", invalid="ignore"):
871 g_a = np.where(safe, (g_lo - 2 * g_mid + g_hi) / (2 * half**2), 0.0)
872 g_b = np.where(safe, (g_hi - g_lo) / (2 * half), 0.0)
873 l_b = np.where(safe, (l_hi - l_lo) / (2 * half), 0.0)
874 int_g = g_a * mu2 + g_b * mu1 + g_mid * m0
875 int_l = l_b * mu1 + l_mid * m0
877 num = int_g.sum(axis=1)
878 den = int_l.sum(axis=1)
879 covered = good & (den > 0)
880 if spinup_threshold > 0.0:
881 # float spinup: emit only where the covered fraction of the APVD reaches the threshold.
882 # den is the covered sub-mass; (v_hi - v_lo) * m0.sum is the fully-covered reference (equal
883 # up to reassociation ulps when fully covered). The relative slack keeps the strictest
884 # threshold (1.0) from rejecting fully-covered bins on that float noise.
885 covered &= den >= (spinup_threshold - _SPINUP_GATE_RELTOL) * (v_hi - v_lo) * m0.sum(axis=1)
886 tile_result = np.full(nt, np.nan)
887 tile_result[covered] = num[covered] / den[covered]
888 result[t0:t1] = tile_result
890 # Zero-throughflow output bins (Q = 0 over the bin) have a cumulative-volume window only as wide
891 # as the strictly-monotone ulp bump, so the bin-average num/den ratio above catastrophically
892 # cancels. There the bin-average degenerates to its well-defined zero-width-bin limit: the
893 # pointwise gamma-mean residence time at the bin's cumulative volume (matching full).
894 dvol = vol_out[1:] - vol_out[:-1]
895 tol = 1e6 * np.spacing(float(np.max(np.abs(flow_cum)))) if flow_cum.size else 0.0
896 degenerate = good_all & (dvol <= tol)
897 if np.any(degenerate):
898 v = vol_out[:-1][degenerate] # (k,) bin cumulative volume (constant over a zero-flow bin)
899 # Over a zero-flow bin the volume is fixed but output time advances, so tau ramps linearly with
900 # output time; its bin-average is the value at the bin's time midpoint, not at T(v).
901 t_v = (0.5 * (cout_tedges_days[:-1] + cout_tedges_days[1:]))[degenerate]
902 # Upper V_p integration bound, mirroring the main loop's spin-up handling. With
903 # spinup="constant" the extrapolated phi map warm-starts the full support; otherwise a V_p
904 # whose look-back/forward parcel leaves the record (e2i: v - r*V_p < 0, i2e: v + r*V_p >
905 # v_end) is in spin-up, so integrate only the covered sub-range [support_lo, vp_hi] and let the
906 # covered sub-mass renormalize the mean (den_pt below), keeping the parcel inside the raw phi map.
907 if extrapolate:
908 vp_hi = np.full(v.shape, support_hi)
909 else:
910 vp_hi = np.clip((v if sign < 0 else v_end - v) / r, support_lo, support_hi)
911 # tau(v, V_p) = sign * (T(v + sign*r*V_p) - t_mid) is piecewise-linear in V_p with knots where
912 # v + sign*r*V_p crosses a phi_v edge; integrate it against the gamma density over the sub-range.
913 bp = np.clip(sign * (phi_v[None, :] - v[:, None]) / r, support_lo, vp_hi[:, None])
914 bp.sort(axis=1)
915 edges_pt = np.concatenate([np.full((v.size, 1), support_lo), bp, vp_hi[:, None]], axis=1)
916 lo_pt, hi_pt = edges_pt[:, :-1], edges_pt[:, 1:]
917 nodes_pt = np.stack([lo_pt, 0.5 * (lo_pt + hi_pt), hi_pt]) # (3, k, n_pieces)
918 a_pt = v[None, :, None] + sign * r * nodes_pt
919 tau_lo, tau_mid, tau_hi = sign * (
920 np.interp(a_pt.ravel(), phi_v, phi_t).reshape(a_pt.shape) - t_v[None, :, None]
921 )
922 width = np.where(hi_pt > lo_pt, hi_pt - lo_pt, 1.0)
923 slope = np.where(hi_pt > lo_pt, (tau_hi - tau_lo) / width, 0.0)
924 mid_pt = 0.5 * (lo_pt + hi_pt)
925 cdf_pt = edges_pt - loc
926 m0_pt = np.diff(gamma_dist.cdf(cdf_pt, alpha, scale=beta), axis=1)
927 m1_pt = alpha * beta * np.diff(gamma_dist.cdf(cdf_pt, alpha + 1, scale=beta), axis=1) + loc * m0_pt
928 num_pt = (tau_mid * m0_pt + slope * (m1_pt - mid_pt * m0_pt)).sum(axis=1)
929 den_pt = m0_pt.sum(axis=1) # covered sub-mass over [support_lo, vp_hi]
930 covered_pt = den_pt > 0
931 if spinup_threshold > 0.0:
932 # float spinup: emit only where the covered fraction reaches the threshold. den_pt is the
933 # covered sub-mass; the full-support mass is the fully-covered reference. The same relative
934 # slack as the main loop keeps threshold 1.0 robust to reassociation ulps in den_pt.
935 full_mass = gamma_dist.cdf(support_hi - loc, alpha, scale=beta) - gamma_dist.cdf(
936 support_lo - loc, alpha, scale=beta
937 )
938 covered_pt &= den_pt >= (spinup_threshold - _SPINUP_GATE_RELTOL) * full_mass
939 with np.errstate(divide="ignore", invalid="ignore"):
940 result[degenerate] = np.where(covered_pt, num_pt / den_pt, np.nan)
941 return result
944def fraction_explained_full(
945 *,
946 flow: npt.ArrayLike,
947 tedges: pd.DatetimeIndex | np.ndarray,
948 cout_tedges: pd.DatetimeIndex | np.ndarray,
949 aquifer_pore_volumes: npt.ArrayLike,
950 direction: str = "extraction_to_infiltration",
951 retardation_factor: float = 1.0,
952) -> npt.NDArray[np.floating]:
953 r"""
954 Advective coverage per pore volume: the fraction of each output bin explained by the record.
956 For each streamtube (entry in ``aquifer_pore_volumes``) and each output bin
957 ``[cout_tedges[i], cout_tedges[i + 1])`` this returns the flow-weighted fraction of the bin whose
958 retarded **advective** parcel lies inside the supplied flow record -- the share of the bin's
959 throughflow volume for which the look-back infiltration (``extraction_to_infiltration``) or
960 look-forward extraction (``infiltration_to_extraction``) event is covered by ``cin``. ``1.0``
961 means the whole bin is explained for that pore volume, ``0.0`` that none of it is. The full
962 ``(n_pore_volumes, n_output_bins)`` array is returned -- one row per pore volume, mirroring
963 :func:`full`.
965 .. warning::
967 This is a **purely advective** diagnostic: it uses only the cumulative-volume look-back
968 ``V(t) - retardation_factor * V_p`` and ignores molecular diffusion and longitudinal
969 dispersion. Those spread each output bin over a *range* of infiltration times whose kernel
970 tails extend outside any finite record, so a bin that is advectively "fully explained"
971 (``1.0``) is not fully informed once dispersion is present. For the dispersive informed
972 fraction of an advection-dispersion model use the captured kernel mass (the column sum of the
973 diffusion coefficient matrix), not this function.
975 Parameters
976 ----------
977 flow : array-like
978 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one.
979 tedges : pandas.DatetimeIndex
980 Time edges for the flow data; ``n + 1`` edges for ``n`` flow values.
981 cout_tedges : pandas.DatetimeIndex
982 Output time edges; ``n + 1`` edges define ``n`` output bins.
983 aquifer_pore_volumes : float or array-like
984 Pore volume(s) of the aquifer [m³], one per streamtube.
985 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
986 Direction of the flow calculation. Default is 'extraction_to_infiltration'.
987 retardation_factor : float, optional
988 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
990 Returns
991 -------
992 numpy.ndarray
993 Advective coverage [dimensionless], shape ``(n_pore_volumes, n_output_bins)``, values in
994 ``[0, 1]``. Output bins lying wholly outside ``tedges`` are ``NaN``. Negative or ``NaN``
995 ``flow`` makes the cumulative-volume map non-monotone or undefined; the whole array is
996 returned as ``NaN`` (the function refuses rather than raising).
998 Raises
999 ------
1000 ValueError
1001 If ``tedges`` does not have exactly one more element than ``flow``, or if ``direction`` is not
1002 ``'extraction_to_infiltration'`` or ``'infiltration_to_extraction'``.
1004 See Also
1005 --------
1006 fraction_explained_mean : Equal-weight mean of this over a discrete APVD
1007 fraction_explained_gamma : Closed-form coverage for a (shifted) gamma APVD
1008 full : Per-pore-volume mean residence time over output bins
1009 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
1010 """
1011 if direction not in {"extraction_to_infiltration", "infiltration_to_extraction"}:
1012 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
1013 raise ValueError(msg)
1015 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes)
1016 tedges = pd.DatetimeIndex(tedges)
1017 cout_tedges = pd.DatetimeIndex(cout_tedges)
1018 flow = np.asarray(flow, dtype=float)
1020 if len(tedges) != len(flow) + 1:
1021 msg = "tedges must have one more element than flow"
1022 raise ValueError(msg)
1024 n_out = len(cout_tedges) - 1
1025 # Negative or non-finite flow makes V(t) non-monotone or undefined; refuse to answer (match siblings).
1026 if np.any(flow < 0) or np.any(np.isnan(flow)):
1027 return np.full((len(aquifer_pore_volumes), n_out), np.nan)
1029 tedges_days = tedges_to_days(tedges)
1030 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0])
1031 flow_cum = cumulative_flow_volume(flow, np.diff(tedges_days))
1032 v_total = flow_cum[-1]
1034 vol_out = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days, left=np.nan, right=np.nan)
1035 v_lo = vol_out[:-1]
1036 v_hi = vol_out[1:]
1037 dvol = v_hi - v_lo
1038 r_vp = retardation_factor * aquifer_pore_volumes
1040 # The flow-weighted (uniform-in-volume) coverage of [v_lo, v_hi] is a clipped ramp of the
1041 # in-record indicator: e2i needs V >= R*V_p, i2e needs V <= v_total - R*V_p. A zero-throughflow
1042 # bin (dvol <= tol) has no volume to average over, so use the pointwise indicator at its volume.
1043 tol = 1e6 * np.spacing(float(np.max(np.abs(flow_cum)))) if flow_cum.size else 0.0
1044 with np.errstate(divide="ignore", invalid="ignore"):
1045 if direction == "extraction_to_infiltration":
1046 frac = np.clip((v_hi[None, :] - r_vp[:, None]) / dvol[None, :], 0.0, 1.0)
1047 point = (r_vp[:, None] <= v_lo[None, :]).astype(float)
1048 else:
1049 frac = np.clip(((v_total - r_vp[:, None]) - v_lo[None, :]) / dvol[None, :], 0.0, 1.0)
1050 point = (r_vp[:, None] <= v_total - v_lo[None, :]).astype(float)
1051 out = np.where(dvol[None, :] > tol, frac, point)
1052 out[:, ~(np.isfinite(v_lo) & np.isfinite(v_hi))] = np.nan
1053 return out
1056def fraction_explained_mean(
1057 *,
1058 flow: npt.ArrayLike,
1059 tedges: pd.DatetimeIndex | np.ndarray,
1060 cout_tedges: pd.DatetimeIndex | np.ndarray,
1061 aquifer_pore_volumes: npt.ArrayLike,
1062 direction: str = "extraction_to_infiltration",
1063 retardation_factor: float = 1.0,
1064) -> npt.NDArray[np.floating]:
1065 """
1066 Advective coverage for a discrete APVD: equal-weight mean of :func:`fraction_explained_full`.
1068 Collapses the pore-volume axis of :func:`fraction_explained_full` to a single per-output-bin
1069 series by averaging over the equally-weighted streamtubes in ``aquifer_pore_volumes`` -- the
1070 coverage analogue of :func:`mean`. ``1.0`` means every streamtube fully explains the
1071 bin, ``0.0`` that none do.
1073 .. warning::
1075 Purely advective -- see :func:`fraction_explained_full`. Molecular diffusion and longitudinal
1076 dispersion spreading are not captured, so a value of ``1.0`` is advective coverage, not full
1077 dispersive information.
1079 Parameters
1080 ----------
1081 flow : array-like
1082 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one.
1083 tedges : pandas.DatetimeIndex
1084 Time edges for the flow data; ``n + 1`` edges for ``n`` flow values.
1085 cout_tedges : pandas.DatetimeIndex
1086 Output time edges; ``n + 1`` edges define ``n`` output bins.
1087 aquifer_pore_volumes : array-like
1088 Discrete pore volumes [m³], one per (equally-weighted) streamtube.
1089 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
1090 Direction of the flow calculation. Default is 'extraction_to_infiltration'.
1091 retardation_factor : float, optional
1092 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
1094 Returns
1095 -------
1096 numpy.ndarray
1097 Advective coverage [dimensionless], shape ``(n_output_bins,)``, values in ``[0, 1]``.
1098 Output bins lying wholly outside ``tedges`` are ``NaN``. Negative or ``NaN`` ``flow`` makes
1099 the cumulative-volume map non-monotone or undefined; the whole series is returned as ``NaN``
1100 (the function refuses rather than raising).
1102 See Also
1103 --------
1104 fraction_explained_full : Per-pore-volume coverage (the array this averages)
1105 fraction_explained_gamma : Closed-form coverage for a (shifted) gamma APVD
1106 mean : Equally-weighted mean residence time for a discrete APVD
1108 Examples
1109 --------
1110 >>> import numpy as np
1111 >>> import pandas as pd
1112 >>> from gwtransport.residence_time import fraction_explained_mean
1113 >>> tedges = pd.date_range("2020-01-01", periods=11, freq="D")
1114 >>> flow = np.full(10, 100.0)
1115 >>> fraction_explained_mean(
1116 ... flow=flow,
1117 ... tedges=tedges,
1118 ... cout_tedges=tedges,
1119 ... aquifer_pore_volumes=[200.0, 1500.0],
1120 ... ).tolist()
1121 [0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
1122 """
1123 return fraction_explained_full(
1124 flow=flow,
1125 tedges=tedges,
1126 cout_tedges=cout_tedges,
1127 aquifer_pore_volumes=aquifer_pore_volumes,
1128 direction=direction,
1129 retardation_factor=retardation_factor,
1130 ).mean(axis=0)
1133def fraction_explained_gamma(
1134 *,
1135 flow: npt.ArrayLike,
1136 tedges: pd.DatetimeIndex | np.ndarray,
1137 cout_tedges: pd.DatetimeIndex | np.ndarray,
1138 mean: float | None = None,
1139 std: float | None = None,
1140 loc: float = 0.0,
1141 alpha: float | None = None,
1142 beta: float | None = None,
1143 direction: str = "extraction_to_infiltration",
1144 retardation_factor: float = 1.0,
1145) -> npt.NDArray[np.floating]:
1146 r"""
1147 Closed-form advective coverage for a (shifted) gamma APVD.
1149 The expectation of the advective in-record indicator over a (shifted) gamma aquifer pore-volume
1150 distribution (APVD), parameterized by either ``(mean, std, loc)`` or ``(alpha, beta, loc)``, is
1151 taken in closed form -- the continuum analogue of :func:`fraction_explained_mean`, with no
1152 pore-volume binning. For each output bin it returns the flow-weighted fraction of the bin whose
1153 advective parcel lies inside the flow record.
1155 The flow-weighted bin average :math:`\frac{1}{\Delta V}\int_{V_\mathrm{lo}}^{V_\mathrm{hi}}
1156 F_{V_p}(\mathrm{threshold}(V))\,dV` (with :math:`\mathrm{threshold}(V) = V / R` for
1157 ``extraction_to_infiltration`` and :math:`(V_\mathrm{end} - V) / R` for
1158 ``infiltration_to_extraction``) is evaluated from the antiderivative of the shifted-gamma CDF,
1160 .. math::
1162 \Phi(x) = \int_\mathrm{loc}^{x} F_{V_p}(s)\,ds
1163 = y\,P(\alpha, y/\beta) - \alpha\beta\,P(\alpha + 1, y/\beta),
1164 \qquad y = \max(x - \mathrm{loc},\, 0),
1166 with :math:`P` the regularized lower incomplete gamma function -- two CDF evaluations per output
1167 edge, no quadrature and no pore-volume binning.
1169 .. warning::
1171 Purely advective -- see :func:`fraction_explained_full`. Molecular diffusion and longitudinal
1172 dispersion are not captured; a value of ``1.0`` is advective coverage, not full dispersive
1173 information.
1175 Parameters
1176 ----------
1177 flow : array-like
1178 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one.
1179 tedges : pandas.DatetimeIndex
1180 Time edges for the flow data; ``n + 1`` edges for ``n`` flow values.
1181 cout_tedges : pandas.DatetimeIndex
1182 Output time edges; ``n + 1`` edges define ``n`` output bins.
1183 mean : float, optional
1184 Mean of the gamma APVD [m³]. Must be strictly greater than ``loc``. Provide either
1185 ``(mean, std)`` or ``(alpha, beta)``.
1186 std : float, optional
1187 Standard deviation of the gamma APVD [m³]. Must be positive.
1188 loc : float, optional
1189 Location (lower bound of support) of the gamma APVD [m³]. Must satisfy ``0 <= loc < mean``.
1190 Default is 0.0.
1191 alpha : float, optional
1192 Shape parameter of the gamma APVD (must be > 0).
1193 beta : float, optional
1194 Scale parameter of the gamma APVD (must be > 0).
1195 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
1196 Direction of the flow calculation. Default is 'extraction_to_infiltration'.
1197 retardation_factor : float, optional
1198 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
1200 Returns
1201 -------
1202 numpy.ndarray
1203 Advective coverage [dimensionless], shape ``(n_output_bins,)``, values in ``[0, 1]``.
1204 Output bins lying wholly outside ``tedges`` are ``NaN``. Negative or ``NaN`` ``flow`` makes
1205 the cumulative-volume map non-monotone or undefined; the whole series is returned as ``NaN``
1206 (the function refuses rather than raising).
1208 Raises
1209 ------
1210 ValueError
1211 If ``tedges`` does not have exactly one more element than ``flow``, or if ``direction`` is
1212 not ``'extraction_to_infiltration'`` or ``'infiltration_to_extraction'``. Gamma parameter
1213 validation is delegated to :func:`gwtransport.gamma.parse_parameters`.
1215 See Also
1216 --------
1217 fraction_explained_mean : Discrete equal-weight APVD coverage
1218 fraction_explained_full : Per-pore-volume coverage
1219 gamma : Closed-form mean residence time for a (shifted) gamma APVD
1220 :ref:`concept-gamma-distribution` : Two-parameter pore volume model
1221 """
1222 if direction not in {"extraction_to_infiltration", "infiltration_to_extraction"}:
1223 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
1224 raise ValueError(msg)
1225 alpha, beta, loc = parse_parameters(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta)
1227 tedges = pd.DatetimeIndex(tedges)
1228 cout_tedges = pd.DatetimeIndex(cout_tedges)
1229 flow = np.asarray(flow, dtype=float)
1231 if len(tedges) != len(flow) + 1:
1232 msg = "tedges must have one more element than flow"
1233 raise ValueError(msg)
1235 n_out = len(cout_tedges) - 1
1236 if np.any(flow < 0) or np.any(np.isnan(flow)):
1237 return np.full(n_out, np.nan)
1239 r = retardation_factor
1240 tedges_days = tedges_to_days(tedges)
1241 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0])
1242 flow_cum = cumulative_flow_volume(flow, np.diff(tedges_days))
1243 v_total = flow_cum[-1]
1245 vol_out = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days, left=np.nan, right=np.nan)
1246 v_lo = vol_out[:-1]
1247 v_hi = vol_out[1:]
1248 dvol = v_hi - v_lo
1250 # threshold(V) per output edge, evaluated once on the full edge array so each shared interior
1251 # edge's CDF is computed a single time (the per-bin lo/hi are then slices of these):
1252 # Phi(x) = int_loc^x F_Vp(s) ds = y P(alpha, y/beta) - alpha beta P(alpha+1, y/beta), y = max(x-loc, 0)
1253 # cdf(x) = F_Vp(x) = P(alpha, y/beta)
1254 threshold = (vol_out if direction == "extraction_to_infiltration" else v_total - vol_out) / r
1255 y = np.maximum(threshold - loc, 0.0)
1256 cdf_edge = gamma_dist.cdf(y, alpha, scale=beta)
1257 phi_edge = y * cdf_edge - alpha * beta * gamma_dist.cdf(y, alpha + 1, scale=beta)
1259 # Flow-weighted bin average (1/dvol) int F_Vp(threshold(V)) dV = (R/dvol) [Phi(hi) - Phi(lo)]; the
1260 # i2e threshold decreases in V so its edge order flips. A zero-throughflow bin (dvol <= tol)
1261 # degenerates to the pointwise CDF at the bin's lower-edge volume.
1262 tol = 1e6 * np.spacing(float(np.max(np.abs(flow_cum)))) if flow_cum.size else 0.0
1263 with np.errstate(divide="ignore", invalid="ignore"):
1264 if direction == "extraction_to_infiltration":
1265 ratio = (r / dvol) * (phi_edge[1:] - phi_edge[:-1])
1266 else:
1267 ratio = (r / dvol) * (phi_edge[:-1] - phi_edge[1:])
1268 point = cdf_edge[:-1]
1269 out = np.where(dvol > tol, ratio, point)
1270 out[~(np.isfinite(v_lo) & np.isfinite(v_hi))] = np.nan
1271 return out
1274def freundlich_retardation(
1275 *,
1276 concentration: npt.ArrayLike,
1277 freundlich_k: float,
1278 freundlich_n: float,
1279 bulk_density: float,
1280 porosity: float,
1281) -> npt.NDArray[np.floating]:
1282 """
1283 Compute concentration-dependent retardation factors using Freundlich isotherm.
1285 The Freundlich isotherm relates sorbed concentration s to aqueous concentration C using the
1286 heterogeneity-index convention (matching :class:`gwtransport.fronttracking.math.FreundlichSorption`
1287 and :func:`gwtransport.advection.infiltration_to_extraction_nonlinear_sorption`, so a fitted
1288 ``freundlich_n`` is portable across the package)::
1290 s = k_f * C ^ (1 / n)
1292 The retardation factor is computed as::
1294 R = 1 + (rho_b/θ) * ds/dC = 1 + (rho_b/θ) * k_f * (1/n) * C^(1/n - 1)
1296 Parameters
1297 ----------
1298 concentration : array-like
1299 Concentration of compound in water [mass/volume]. One value per time bin, consistent
1300 with the ``flow`` array passed to the transport function.
1301 freundlich_k : float
1302 Freundlich coefficient [(m³/kg)^(1/n)] (under s = k_f * C^(1/n) with s dimensionless
1303 and C in [kg/m³]).
1304 freundlich_n : float
1305 Freundlich sorption exponent [dimensionless] (heterogeneity index; ``n = 1`` recovers a
1306 linear isotherm).
1307 bulk_density : float
1308 Bulk density of aquifer material [mass/volume].
1309 porosity : float
1310 Porosity of aquifer [dimensionless, 0-1].
1312 Returns
1313 -------
1314 numpy.ndarray
1315 Retardation factors for each flow interval.
1316 Length equals len(concentration) for use as retardation_factor in the transport functions.
1318 Raises
1319 ------
1320 ValueError
1321 If ``porosity`` is not in ``(0, 1)``, if ``bulk_density`` is not positive, if
1322 ``freundlich_k`` is negative, or if any ``concentration`` is non-positive while
1323 ``freundlich_n > 1`` (the retardation factor diverges as ``C -> 0``).
1325 See Also
1326 --------
1327 full : Compute residence times from flow and pore volume
1328 gwtransport.advection.infiltration_to_extraction_nonlinear_sorption : Transport with nonlinear sorption
1329 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and concentration-dependent retardation
1331 Examples
1332 --------
1333 >>> concentration = np.array([0.1, 0.2, 0.3]) # same length as flow
1334 >>> R = freundlich_retardation(
1335 ... concentration=concentration,
1336 ... freundlich_k=0.5,
1337 ... freundlich_n=2.0,
1338 ... bulk_density=1600, # kg/m³
1339 ... porosity=0.35,
1340 ... )
1341 >>> # Use R as retardation_factor in the transport functions
1342 """
1343 concentration = np.asarray(concentration)
1345 if not 0 < porosity < 1:
1346 msg = f"Porosity must be in (0, 1), got {porosity}"
1347 raise ValueError(msg)
1348 if bulk_density <= 0:
1349 msg = f"Bulk density must be positive, got {bulk_density}"
1350 raise ValueError(msg)
1351 if freundlich_k < 0:
1352 msg = f"Freundlich K must be non-negative, got {freundlich_k}"
1353 raise ValueError(msg)
1355 # For n > 1 the Freundlich retardation factor 1 + (rho_b/theta) * k_f * (1/n) * C^(1/n-1)
1356 # diverges as C -> 0 (the exponent 1/n - 1 < 0). Silently clamping concentration would produce
1357 # a very large but finite value that depends on an arbitrary regularization constant; instead,
1358 # refuse the call so the user can decide how to handle non-positive concentrations.
1359 if freundlich_n > 1.0 and np.any(concentration <= 0):
1360 msg = "concentration must be strictly positive when freundlich_n > 1 (retardation diverges as C -> 0)"
1361 raise ValueError(msg)
1363 inv_n = 1.0 / freundlich_n
1364 return 1.0 + (bulk_density / porosity) * freundlich_k * inv_n * concentration ** (inv_n - 1.0)