Coverage for src/gwtransport/fronttracking/math.py: 97%
352 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""
2Mathematical Foundation for Front Tracking with Nonlinear Sorption.
4This module provides exact analytical computations for:
6- Freundlich, Langmuir, and constant retardation models
7- Brooks-Corey and van Genuchten-Mualem unsaturated conductivity models
8 (for Kinematic-Wave percolation, see :mod:`gwtransport.percolation`)
9- Shock velocities via Rankine-Hugoniot condition
10- Characteristic velocities and positions
11- First arrival time calculations
12- Entropy condition verification
14All sorption-class computations are exact analytical formulas; the
15van Genuchten-Mualem class uses ``scipy.optimize.brentq`` for the two
16inversions that have no closed form.
18Available functions:
20- :class:`NonlinearSorption` - Abstract base for concentration-dependent isotherms. Subclasses supply
21 ``retardation(c)``, ``total_concentration(c)`` and ``concentration_from_retardation(r)``; the base derives
22 the Rankine-Hugoniot ``shock_speed(c_left, c_right)``, the Lax ``check_entropy_condition``, the paired
23 ``c_and_total_from_retardation`` and ``fan_converges_at_infinity`` from them.
25- :class:`FreundlichSorption` - Power-law isotherm ``s(C) = k_f · C^(1/n)`` with closed-form
26 ``R(C) = 1 + (ρ_b · k_f)/(n_por · n) · C^(1/n − 1)``. ``n > 1`` makes higher concentrations travel faster,
27 ``n < 1`` mirrors that, and ``n = 1`` is rejected (use :class:`ConstantRetardation`). Because ``R → ∞`` as
28 ``C → 0`` for ``n > 1``, ``retardation`` and ``concentration_from_retardation`` clamp ``C`` at ``c_min``.
30- :class:`ConstantRetardation` - Linear sorption ``s(C) = K_d · C`` reduced to one ``retardation_factor``:
31 every concentration moves at ``dV/dθ = 1/R``, so no rarefaction ever forms and the solution is a pure
32 θ-shift with shocks only at inlet discontinuities.
34- :class:`LangmuirSorption` - Favorable isotherm ``s(C) = s_max · C/(K_L + C)`` with
35 ``R(C) = 1 + (ρ_b · s_max · K_L)/(n_por · (K_L + C)²)``, which is finite at ``C = 0``, so no concentration
36 floor is needed; ``R`` decreases with ``C``, giving shocks on concentration rises and fans on falls.
38- :class:`BrooksCoreyConductivity` - Brooks-Corey unsaturated conductivity ``K(θ) = K_s · Θ^a`` recast into
39 the ``(C, C_T)`` variables of the sorption interface by identifying ``C ≡ K`` (flux) and ``C_T ≡ θ − θ_r``
40 (storage), which lets :mod:`gwtransport.percolation` run Kinematic-Wave percolation on this solver. All
41 three curve methods are closed form; ``C`` is clamped at a dry-soil floor in ``retardation`` and
42 ``concentration_from_retardation`` but not in ``total_concentration``, keeping the ``c_R = 0``
43 wetting-front shock speed exact.
45- :class:`VanGenuchtenMualemConductivity` - Mualem conductivity for the van Genuchten retention curve, recast
46 the same way. ``total_concentration`` and the flux derivative are closed form; the two inversions ``S_e(C)``
47 and ``S_e(R)`` use ``scipy.optimize.brentq``. The retention parameter ``α_vG`` is not needed because the
48 Kinematic-Wave approximation drops capillary suction.
50- :func:`characteristic_speed` - Characteristic speed ``dV/dθ = 1/R(c)`` in (V, θ) coordinates, a property of
51 the isotherm alone. Its ``sorption`` argument is any of the classes above (the ``SorptionModel`` alias is
52 their union).
54- :func:`characteristic_position` - Position of a characteristic at cumulative flow ``theta``, evaluated as
55 ``v_start + (θ − θ_start)/R(c)``, or ``None`` when ``theta`` lies before ``theta_start``.
57- :func:`compute_first_front_arrival_theta` - Cumulative flow at which the first injected concentration level
58 is fully present at the outlet: the θ-edge of the first bin that both carries water (positive θ-width) and
59 has ``cin > 0``, plus ``V · R(c_first)`` for Freundlich ``n < 1`` and ``V · C_T(c_first)/c_first``
60 otherwise. These are *tail*-arrival semantics — a rarefaction head can reach the outlet much earlier — and
61 the result is ``inf`` when no water-carrying bin injects concentration.
63This file is part of gwtransport which is released under AGPL-3.0 license.
64See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
65"""
67import math
68from abc import ABC, abstractmethod
69from dataclasses import dataclass, field
71import numpy as np
72import numpy.typing as npt
73from scipy.optimize import brentq
75# Numerical tolerance constants
76EPSILON_FREUNDLICH_N = 1e-10 # Tolerance for checking if n ≈ 1.0 (Freundlich constructor rejects this)
77EPSILON_DENOMINATOR = 1e-15 # Tolerance for near-zero denominators in shock velocity
78_C_MIN = 1e-12 # Shared dry-soil singularity floor for Freundlich n>1, Brooks-Corey, vG-Mualem.
79BRENTQ_XTOL = 1e-14 # brentq absolute tolerance for vG-Mualem inversions; matches _invert_freundlich_cr_zero.
82class NonlinearSorption(ABC):
83 """Abstract base for concentration-dependent sorption models.
85 Subclasses must implement `retardation`, `total_concentration`, and
86 `concentration_from_retardation`. Shock velocity and entropy checking
87 are provided generically via the Rankine-Hugoniot and Lax conditions.
89 See Also
90 --------
91 FreundlichSorption : Freundlich isotherm implementation.
92 LangmuirSorption : Langmuir isotherm implementation.
93 ConstantRetardation : Linear (constant R) retardation model.
94 """
96 @abstractmethod
97 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
98 """Compute retardation factor R(C).
100 R fixes the characteristic celerity in (V, θ) coordinates: ``dV/dθ = 1/R(C)``.
102 Parameters
103 ----------
104 c : float or numpy.ndarray
105 Dissolved concentration [mass/volume]. Non-negative.
107 Returns
108 -------
109 float or numpy.ndarray
110 Retardation factor [-].
111 """
113 @abstractmethod
114 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
115 """Compute total concentration (dissolved + sorbed per unit pore volume).
117 ``C_T`` is the conserved quantity of ``∂C_T/∂θ + ∂C/∂V = 0``; the flux carries
118 only the dissolved part because sorbed mass is immobile.
120 Parameters
121 ----------
122 c : float or numpy.ndarray
123 Dissolved concentration [mass/volume]. Non-negative.
125 Returns
126 -------
127 float or numpy.ndarray
128 Total concentration [mass/volume].
129 """
131 @abstractmethod
132 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
133 """Invert retardation factor to obtain concentration.
135 Used by rarefaction fans, where the self-similar solution gives R as a
136 function of position and cumulative flow.
138 Parameters
139 ----------
140 r : float or numpy.ndarray
141 Retardation factor [-].
143 Returns
144 -------
145 float or numpy.ndarray
146 Dissolved concentration [mass/volume]. Non-negative.
147 """
149 def shock_speed(self, c_left: float, c_right: float) -> float:
150 """Compute shock speed dV/dθ via Rankine-Hugoniot in (V, θ) coordinates.
152 With cumulative-flow coordinate θ = ∫flow(t') dt', the PDE
153 ``∂C_T/∂t + flow·∂C/∂V = 0`` becomes ``∂C_T/∂θ + ∂C/∂V = 0``, and
154 Rankine-Hugoniot reduces to::
156 dV_s/dθ = (C_R - C_L) / (C_T(C_R) - C_T(C_L))
158 Flow drops out entirely; the result is a property of the sorption
159 isotherm alone.
161 Parameters
162 ----------
163 c_left : float
164 Concentration upstream (behind) shock [mass/volume].
165 c_right : float
166 Concentration downstream (ahead of) shock [mass/volume].
168 Returns
169 -------
170 shock_speed : float
171 Shock speed dV/dθ [m³ / m³ flow = dimensionless].
172 """
173 c_total_left = self.total_concentration(c_left)
174 c_total_right = self.total_concentration(c_right)
175 denom = c_total_right - c_total_left
177 if abs(denom) < EPSILON_DENOMINATOR:
178 avg_retardation = 0.5 * float(self.retardation(c_left) + self.retardation(c_right))
179 # Degenerate (zero-strength) shock: its speed is the characteristic speed 1/R. A pair of
180 # saturated states (R = 0, e.g. Mualem-vG at S_e = 1) gives +∞, matching characteristic_speed.
181 return float("inf") if avg_retardation == 0.0 else 1.0 / avg_retardation
183 return float((c_right - c_left) / denom)
185 def c_and_total_from_retardation(self, r: float) -> tuple[float, float]:
186 """Return ``(c, C_T(c))`` at a given retardation ``r``.
188 Default implementation calls ``concentration_from_retardation(r)`` then
189 ``total_concentration(c)`` — two independent root-finds for sorptions
190 where both routes back-solve the same equation (e.g. vG-Mualem with
191 ``L ≠ 0``). Subclasses for which both can be computed from a single
192 root-find should override this for ~2× speedup of the IBP fan
193 integrators.
194 """
195 c = float(self.concentration_from_retardation(r))
196 ct = float(self.total_concentration(c))
197 return c, ct
199 def fan_converges_at_infinity(self) -> bool: # noqa: PLR6301
200 """Whether a ``c_apex=0`` fan's ``∫ c dθ`` converges as ``θ → +∞``.
202 True when ``c → 0`` as ``R → ∞`` (so ``base·c → 0`` faster than ``base → ∞``):
203 Brooks-Corey, van Genuchten-Mualem, Langmuir, and Freundlich ``n > 1``. The
204 only divergent case is Freundlich ``n < 1`` (``c → ∞`` as ``R → ∞``), which
205 overrides this to ``False``. Used by the universal temporal fan integrator to
206 reject a ``+∞`` upper bound when the integral diverges.
207 """
208 return True
210 def check_entropy_condition(self, c_left: float, c_right: float, shock_speed: float) -> bool:
211 """Verify Lax entropy condition in (V, θ) coordinates.
213 In θ-space, characteristic speeds are ``λ_θ(C) = 1 / R(C)``, and the
214 Lax condition for a physical shock is::
216 λ_θ(C_L) ≥ dV_s/dθ ≥ λ_θ(C_R)
218 Parameters
219 ----------
220 c_left : float
221 Concentration upstream of shock [mass/volume].
222 c_right : float
223 Concentration downstream of shock [mass/volume].
224 shock_speed : float
225 Shock speed dV/dθ.
227 Returns
228 -------
229 satisfies : bool
230 True if shock satisfies entropy condition (is physical).
231 """
232 r_left = float(self.retardation(c_left))
233 r_right = float(self.retardation(c_right))
234 lambda_left = float("inf") if r_left == 0.0 else 1.0 / r_left
235 lambda_right = float("inf") if r_right == 0.0 else 1.0 / r_right
237 # A saturated upstream state (λ_left = +∞, e.g. a Mualem-vG wetting front at S_e = 1) is
238 # physical; reject only a non-finite shock speed or downstream characteristic, where the
239 # Lax test itself is ill-posed.
240 if not np.isfinite(shock_speed) or not np.isfinite(lambda_right):
241 return False
243 finite_left = abs(lambda_left) if np.isfinite(lambda_left) else 0.0
244 tolerance = 1e-14 * max(finite_left, abs(lambda_right), abs(shock_speed))
246 return bool((lambda_left > shock_speed - tolerance) and (shock_speed > lambda_right - tolerance))
249@dataclass
250class FreundlichSorption(NonlinearSorption):
251 """
252 Freundlich sorption isotherm with exact analytical methods.
254 The Freundlich isotherm is: s(C) = k_f * C^(1/n)
256 where:
257 - s is sorbed concentration [mass/mass of solid]
258 - C is dissolved concentration [mass/volume of water]
259 - k_f is Freundlich coefficient [(volume/mass)^(1/n)]
260 - n is Freundlich exponent (dimensionless)
262 For n > 1: Higher C travels faster
263 For n < 1: Higher C travels slower
264 For n = 1: linear (not supported, use ConstantRetardation instead)
266 Notes
267 -----
268 The retardation factor is defined as:
269 R(C) = 1 + (rho_b/n_por) * ds/dC
270 = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)
272 For Freundlich sorption, R depends on C, which creates nonlinear wave behavior.
274 For n>1 (higher C travels faster), R(C)→∞ as C→0, which can cause extremely slow
275 wave propagation. The c_min parameter prevents this by enforcing a minimum
276 concentration, making R(C) finite for all C≥0.
278 Examples
279 --------
280 >>> sorption = FreundlichSorption(
281 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
282 ... )
283 >>> r = sorption.retardation(5.0)
284 >>> c_back = sorption.concentration_from_retardation(r)
285 >>> bool(np.isclose(c_back, 5.0))
286 True
287 """
289 k_f: float
290 """Freundlich coefficient [(m³/kg)^(1/n)]. Positive."""
291 n: float
292 """Freundlich exponent [-]. Positive and != 1."""
293 bulk_density: float
294 """Bulk density of porous medium [kg/m³]. Positive."""
295 porosity: float
296 """Porosity [-]. In (0, 1)."""
297 c_min: float = 1e-12
298 """Dry-soil singularity floor [mass/volume]; keeps ``R(C)`` finite as ``C → 0`` for n>1."""
299 _ret_coefficient: float = field(init=False, repr=False, compare=False)
300 """Cached ``(rho_b*k_f)/(n_por*n)`` — shared by the scalar and array paths."""
301 _ret_exponent: float = field(init=False, repr=False, compare=False)
302 """Cached ``(1/n) - 1`` retardation exponent."""
303 _ct_coefficient: float = field(init=False, repr=False, compare=False)
304 """Cached ``(rho_b/n_por)*k_f`` sorbed-mass coefficient."""
305 _cfr_inv_exponent: float = field(init=False, repr=False, compare=False)
306 """Cached ``1/((1/n) - 1)`` inversion exponent for concentration_from_retardation."""
308 def __post_init__(self):
309 """Validate parameters after initialization.
311 Raises
312 ------
313 ValueError
314 If any parameter is outside its valid range: ``k_f`` <= 0,
315 ``n`` <= 0, ``n`` == 1, ``bulk_density`` <= 0, ``porosity``
316 outside (0, 1), or ``c_min`` < 0.
317 """
318 if self.k_f <= 0:
319 msg = f"k_f must be positive, got {self.k_f}"
320 raise ValueError(msg)
321 if self.n <= 0:
322 msg = f"n must be positive, got {self.n}"
323 raise ValueError(msg)
324 if abs(self.n - 1.0) < EPSILON_FREUNDLICH_N:
325 msg = "n = 1 (linear case) not supported, use ConstantRetardation instead"
326 raise ValueError(msg)
327 if self.bulk_density <= 0:
328 msg = f"bulk_density must be positive, got {self.bulk_density}"
329 raise ValueError(msg)
330 if not 0 < self.porosity < 1:
331 msg = f"porosity must be in (0, 1), got {self.porosity}"
332 raise ValueError(msg)
333 if self.c_min < 0:
334 msg = f"c_min must be non-negative, got {self.c_min}"
335 raise ValueError(msg)
337 self._ret_exponent = (1.0 / self.n) - 1.0
338 self._ret_coefficient = (self.bulk_density * self.k_f) / (self.porosity * self.n)
339 self._ct_coefficient = (self.bulk_density / self.porosity) * self.k_f
340 self._cfr_inv_exponent = 1.0 / self._ret_exponent
342 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
343 """``R(C) = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)``. See :meth:`NonlinearSorption.retardation`.
345 Notes
346 -----
347 R decreases with C for n > 1 (higher C travels faster) and increases for n < 1.
348 With n<1 and ``c_min=0``, ``R(0) = 1`` (no sorption at zero) because clamping to
349 ``c_min=0`` leaves ``C^((1/n)-1) = 0^positive = 0``; otherwise ``c`` is clamped to
350 ``c_min`` before evaluation.
352 Clamping with ``np.maximum`` before the power keeps a single general path
353 for every ``(n, c_min)`` combination and avoids raising the base to a
354 fractional power on negative ``c``.
356 A pure-float fast path handles the (dominant) scalar case, bit-identical to
357 the array path; it falls through to the array expression only when the
358 clamped ``c`` is ``0`` (``c_min == 0`` and ``c <= 0``), where numpy's
359 ``0.0**exponent`` yields the documented ``inf`` (n>1) / ``0`` (n<1) that a
360 pure ``0.0**neg`` would instead raise on.
361 """
362 if not isinstance(c, np.ndarray):
363 cf = float(c)
364 c_eff = max(self.c_min, cf)
365 if c_eff > 0.0:
366 return 1.0 + self._ret_coefficient * (c_eff**self._ret_exponent)
367 c_eff = np.maximum(np.asarray(c), self.c_min)
368 result = 1.0 + self._ret_coefficient * (c_eff**self._ret_exponent)
369 return result if isinstance(c, np.ndarray) else float(result)
371 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
372 """``C_T = C + (rho_b/n_por)·k_f·C^(1/n)``. See :meth:`NonlinearSorption.total_concentration`.
374 Notes
375 -----
376 For ``c = 0``, ``c^(1/n) = 0`` exactly (no singularity for any
377 ``n > 0``), so ``C_T(0) = 0`` is physically correct and no ``c_min``
378 clamp is needed here. ``c_min`` is only required to keep
379 :meth:`retardation` finite as ``c -> 0`` for ``n > 1``; clamping
380 ``total_concentration`` to ``c_min`` would bias Rankine-Hugoniot
381 shock speeds when ``c_R = 0`` (e.g. the canonical 0->c->0 pulse).
382 Negative ``c`` is clamped to ``0`` defensively.
383 """
384 if not isinstance(c, np.ndarray):
385 cf = float(c)
386 c_eff = max(0.0, cf)
387 return c_eff + self._ct_coefficient * (c_eff ** (1.0 / self.n))
388 c_arr = np.maximum(np.asarray(c), 0.0)
389 return c_arr + self._ct_coefficient * (c_arr ** (1.0 / self.n))
391 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
392 """``C = [(R-1)·n_por·n/(rho_b·k_f)]^(n/(1-n))``. See :meth:`NonlinearSorption.concentration_from_retardation`.
394 Notes
395 -----
396 The exponent ``n/(1-n)`` is undefined at n = 1, which is why linear sorption
397 must use :class:`ConstantRetardation` instead. Results below ``c_min`` are
398 clamped to it.
399 """
400 # __post_init__ rejects |n-1| < EPSILON_FREUNDLICH_N, so _cfr_inv_exponent is finite.
401 if not isinstance(r, np.ndarray):
402 base = (float(r) - 1.0) / self._ret_coefficient
403 if base > 0.0:
404 return max(base**self._cfr_inv_exponent, self.c_min)
405 return self.c_min
407 base = (np.asarray(r) - 1.0) / self._ret_coefficient
408 # Mask base to a safe placeholder before exponentiation; NumPy emits
409 # RuntimeWarning otherwise for base <= 0 with a fractional exponent.
410 safe_base = np.where(base > 0, base, 1.0)
411 c = safe_base**self._cfr_inv_exponent
412 return np.where(base > 0, np.maximum(c, self.c_min), self.c_min)
414 def fan_converges_at_infinity(self) -> bool:
415 """Freundlich ``n > 1``: ``c → 0`` as ``R → ∞`` (converges). ``n < 1``: ``c → ∞`` (diverges)."""
416 return self.n > 1.0
419@dataclass
420class ConstantRetardation:
421 """
422 Constant (linear) retardation model.
424 For linear sorption: s(C) = K_d * C
425 This gives constant retardation: R(C) = 1 + (rho_b/n_por) * K_d = constant
427 This is a special case where concentration-dependent behavior disappears.
428 Used for conservative tracers or as approximation for weak sorption.
430 Notes
431 -----
432 With constant retardation:
433 - All concentrations travel at same speed in (V, θ): dV/dθ = 1/R
434 - No rarefaction waves form (all concentrations travel together)
435 - Shocks occur only at concentration discontinuities at inlet
436 - Solution reduces to simple θ-shifting (and then t-shifting via the θ↔t map)
438 This is equivalent to a single-pore-volume advective time-shift (the deterministic limit of
439 :func:`gwtransport.advection.infiltration_to_extraction`) in the gwtransport package.
441 Examples
442 --------
443 >>> sorption = ConstantRetardation(retardation_factor=2.0)
444 >>> sorption.retardation(5.0)
445 2.0
446 >>> sorption.retardation(10.0)
447 2.0
448 """
450 retardation_factor: float
451 """Constant retardation factor [-]. At least 1.0; ``R = 1`` is a conservative tracer."""
453 def __post_init__(self):
454 """Validate parameters after initialization.
456 Raises
457 ------
458 ValueError
459 If ``retardation_factor`` is less than 1.0.
460 """
461 if self.retardation_factor < 1.0:
462 msg = f"retardation_factor must be >= 1.0, got {self.retardation_factor}"
463 raise ValueError(msg)
465 def retardation(self, c: float) -> float: # noqa: ARG002
466 """Constant retardation factor, independent of ``c``."""
467 return self.retardation_factor
469 def total_concentration(self, c: float) -> float:
470 """``C_T = C·R`` for linear sorption."""
471 return c * self.retardation_factor
473 def concentration_from_retardation(self, r: float) -> float:
474 """Not applicable: with constant R the inversion is not meaningful.
476 Raises
477 ------
478 NotImplementedError
479 Always.
480 """
481 msg = "concentration_from_retardation not applicable for ConstantRetardation (R is independent of C)"
482 raise NotImplementedError(msg)
484 def shock_speed(self, c_left: float, c_right: float) -> float: # noqa: ARG002
485 """``dV/dθ = 1/R`` for any concentration pair — identical to every characteristic speed."""
486 return 1.0 / self.retardation_factor
488 def check_entropy_condition(self, c_left: float, c_right: float, shock_speed: float) -> bool: # noqa: PLR6301
489 """Check the Lax entropy condition; always satisfied, since with constant R it holds as an equality.
491 Returns
492 -------
493 satisfies : bool
494 Always True.
495 """
496 del c_left, c_right, shock_speed
497 return True
500@dataclass
501class LangmuirSorption(NonlinearSorption):
502 """
503 Langmuir sorption isotherm with exact analytical methods.
505 The Langmuir isotherm is: s(C) = s_max * C / (K_L + C)
507 where:
508 - s is sorbed concentration [mass/mass of solid]
509 - C is dissolved concentration [mass/volume of water]
510 - s_max is maximum sorption capacity [mass/mass of solid]
511 - K_L is half-saturation constant [mass/volume]
513 Retardation always decreases with C (favorable isotherm), and R(0) is
514 finite — unlike Freundlich with n > 1, no minimum concentration threshold
515 is needed.
517 See Also
518 --------
519 FreundlichSorption : Freundlich isotherm (unbounded sorption).
520 ConstantRetardation : Linear (constant R) retardation model.
521 :ref:`concept-nonlinear-sorption` : Background on nonlinear sorption.
523 Notes
524 -----
525 The retardation factor is defined as:
526 R(C) = 1 + (rho_b * s_max * K_L) / (n_por * (K_L + C)^2)
528 Key properties:
530 - R(0) = 1 + rho_b * s_max / (n_por * K_L) -- finite for all parameters
531 - R -> 1 as C -> infinity (all sorption sites saturated)
532 - R always decreases with increasing C (higher C travels faster)
533 - Shocks form on concentration increases, rarefaction fans on decreases
535 Examples
536 --------
537 >>> sorption = LangmuirSorption(
538 ... s_max=0.1, k_l=5.0, bulk_density=1500.0, porosity=0.3
539 ... )
540 >>> r = sorption.retardation(5.0)
541 >>> c_back = sorption.concentration_from_retardation(r)
542 >>> bool(np.isclose(c_back, 5.0))
543 True
544 """
546 s_max: float
547 """Maximum sorption capacity [mass/mass of solid]. Positive."""
548 k_l: float
549 """Half-saturation constant [mass/volume] — the C at which ``s = s_max/2``. Positive."""
550 bulk_density: float
551 """Bulk density of porous medium [kg/m³]. Positive."""
552 porosity: float
553 """Porosity [-]. In (0, 1)."""
555 def __post_init__(self):
556 """Validate parameters after initialization.
558 Raises
559 ------
560 ValueError
561 If any parameter is outside its valid range: ``s_max`` <= 0,
562 ``k_l`` <= 0, ``bulk_density`` <= 0, or ``porosity``
563 outside (0, 1).
564 """
565 if self.s_max <= 0:
566 msg = f"s_max must be positive, got {self.s_max}"
567 raise ValueError(msg)
568 if self.k_l <= 0:
569 msg = f"k_l must be positive, got {self.k_l}"
570 raise ValueError(msg)
571 if self.bulk_density <= 0:
572 msg = f"bulk_density must be positive, got {self.bulk_density}"
573 raise ValueError(msg)
574 if not 0 < self.porosity < 1:
575 msg = f"porosity must be in (0, 1), got {self.porosity}"
576 raise ValueError(msg)
578 self.a_coeff: float = self.bulk_density * self.s_max * self.k_l / self.porosity
579 """Lumped retardation constant rho_b * s_max * K_L / n_por."""
580 self._ct_coefficient: float = (self.bulk_density / self.porosity) * self.s_max
581 """Cached ``(rho_b/n_por)*s_max`` sorbed-mass coefficient (scalar + array paths)."""
583 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
584 """``R(C) = 1 + A/(K_L + C)²`` with ``A = rho_b·s_max·K_L/n_por``.
586 See :meth:`NonlinearSorption.retardation`. ``R(0)`` is finite, R decreases with
587 C, and ``R → 1`` as ``C → ∞`` (all sorption sites saturated).
588 """
589 if not isinstance(c, np.ndarray):
590 cf = float(c)
591 c_eff = max(0.0, cf)
592 return 1.0 + self.a_coeff / (self.k_l + c_eff) ** 2
593 c_eff = np.maximum(np.asarray(c), 0.0)
594 return 1.0 + self.a_coeff / (self.k_l + c_eff) ** 2
596 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
597 """``C_T = C + (rho_b/n_por)·s_max·C/(K_L + C)``. See :meth:`NonlinearSorption.total_concentration`."""
598 if not isinstance(c, np.ndarray):
599 cf = float(c)
600 c_eff = max(0.0, cf)
601 return cf + self._ct_coefficient * c_eff / (self.k_l + c_eff)
602 c_arr = np.asarray(c)
603 c_eff = np.maximum(c_arr, 0.0)
604 return c_arr + self._ct_coefficient * c_eff / (self.k_l + c_eff)
606 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
607 """``C = sqrt(A/(R - 1)) - K_L``. See :meth:`NonlinearSorption.concentration_from_retardation`.
609 Notes
610 -----
611 Returns ``0.0`` for ``R <= 1`` (unphysical) and for ``R >= R(0) = 1 + A/K_L²``
612 (at or below zero concentration).
613 """
614 if not isinstance(r, np.ndarray):
615 r_minus_1 = float(r) - 1.0
616 if r_minus_1 > 0.0:
617 return max(math.sqrt(self.a_coeff / r_minus_1) - self.k_l, 0.0)
618 return 0.0
620 r_minus_1 = np.asarray(r) - 1.0
621 # Mask r_minus_1 to a safe placeholder before division to avoid the
622 # RuntimeWarning emitted by np.where's eager evaluation when r == 1.
623 safe_r_minus_1 = np.where(r_minus_1 > 0, r_minus_1, 1.0)
624 c = np.where(r_minus_1 > 0, np.sqrt(self.a_coeff / safe_r_minus_1) - self.k_l, 0.0)
625 return np.maximum(c, 0.0)
628@dataclass
629class BrooksCoreyConductivity(NonlinearSorption):
630 r"""Brooks-Corey unsaturated conductivity recast as a NonlinearSorption.
632 Used by :mod:`gwtransport.percolation` to model gravity-driven percolation
633 through a thick unsaturated zone via the Kinematic-Wave method. The
634 closed-form conductivity curve
636 .. math::
637 K(\\theta) = K_s \\cdot \\Theta^a, \\qquad
638 \\Theta = (\\theta - \\theta_r)/(\\theta_s - \\theta_r), \\qquad
639 a = 3 + 2/\\lambda \\;(\\text{Burdine})
641 is recast in the framework's ``(C, C_T)`` variables by identifying
642 ``C ≡ K`` (the flux variable) and ``C_T ≡ θ - θ_r`` (the conserved
643 storage). All three abstract methods have closed forms; ``shock_speed``
644 and ``check_entropy_condition`` are inherited unchanged from
645 :class:`NonlinearSorption`.
647 See Also
648 --------
649 VanGenuchtenMualemConductivity : Van Genuchten variant with brentq inversions.
650 FreundlichSorption : Power-law sorption isotherm (closed form, analogous shape).
651 gwtransport.percolation.root_zone_to_water_table_kinematic_wave : The public wrapper.
653 Notes
654 -----
655 The retardation factor and total-concentration relation are:
657 .. math::
658 C_T(C) = \\Delta\\theta \\cdot (C/K_s)^{1/a}, \\qquad
659 R(C) = (\\Delta\\theta / (a K_s)) \\cdot (C/K_s)^{1/a - 1},
661 with ``Δθ = θ_s − θ_r``. Since ``1/a − 1 < 0`` always (``a > 3``),
662 ``R(C) → ∞`` as ``C → 0`` (dry-soil singularity). The class clamps ``C``
663 to a small floor in ``retardation`` and ``concentration_from_retardation``
664 (the same pattern as :class:`FreundlichSorption` with ``n > 1``);
665 ``total_concentration`` and the inherited ``shock_speed`` do **not**
666 clamp, so the canonical wetting-front shock ``c_R = 0`` produces the
667 correct Rankine-Hugoniot velocity.
669 Examples
670 --------
671 >>> sorption = BrooksCoreyConductivity(
672 ... theta_r=0.01, theta_s=0.337, k_s=0.174, brooks_corey_lambda=0.25
673 ... )
674 >>> r = sorption.retardation(0.05)
675 >>> c = sorption.concentration_from_retardation(r)
676 >>> bool(np.isclose(c, 0.05, rtol=1e-13))
677 True
678 """
680 theta_r: float
681 """Residual volumetric moisture content [-]; ``0 <= theta_r < theta_s``."""
682 theta_s: float
683 """Saturated volumetric moisture content [-]; ``theta_r < theta_s < 1``. The porosity for typical soils."""
684 k_s: float
685 """Saturated hydraulic conductivity [length/time]. Positive."""
686 brooks_corey_lambda: float
687 """Pore-size distribution index ``λ`` [-]. Positive.
689 The exponent ``a = 3 + 2/λ`` is the Burdine pore-connectivity result. The Mualem
690 variant (``L = 0.5``) gives ``a = 2.5 + 2/λ`` and is not implemented; a user wanting
691 it can re-derive ``λ`` so the Burdine ``a`` matches the desired Mualem exponent."""
692 a: float = field(init=False)
693 """Exponent ``a = 3 + 2/λ`` (Burdine); set in ``__post_init__``."""
694 delta_theta: float = field(init=False)
695 """``θ_s − θ_r``; set in ``__post_init__``."""
696 _inv_a: float = field(init=False, repr=False, compare=False)
697 """Cached ``1/a`` total-concentration exponent (scalar + array paths)."""
698 _ret_coefficient: float = field(init=False, repr=False, compare=False)
699 """Cached ``Δθ/(a·K_s)`` retardation coefficient."""
700 _ret_exponent: float = field(init=False, repr=False, compare=False)
701 """Cached ``1/a − 1`` retardation exponent."""
702 _cfr_exponent: float = field(init=False, repr=False, compare=False)
703 """Cached ``−a/(a−1)`` inversion exponent for concentration_from_retardation."""
705 def __post_init__(self) -> None:
706 """Validate parameters and derive ``a``, ``delta_theta``.
708 Raises
709 ------
710 ValueError
711 If any parameter is outside its valid range.
712 """
713 if not 0.0 <= self.theta_r < self.theta_s:
714 msg = f"theta_r must satisfy 0 <= theta_r < theta_s, got theta_r={self.theta_r}, theta_s={self.theta_s}"
715 raise ValueError(msg)
716 if not self.theta_s < 1.0:
717 msg = f"theta_s must be < 1, got {self.theta_s}"
718 raise ValueError(msg)
719 if self.k_s <= 0.0:
720 msg = f"k_s must be positive, got {self.k_s}"
721 raise ValueError(msg)
722 if self.brooks_corey_lambda <= 0.0:
723 msg = f"brooks_corey_lambda must be positive, got {self.brooks_corey_lambda}"
724 raise ValueError(msg)
725 self.a = 3.0 + 2.0 / self.brooks_corey_lambda
726 self.delta_theta = self.theta_s - self.theta_r
727 self._inv_a = 1.0 / self.a
728 self._ret_coefficient = self.delta_theta / (self.a * self.k_s)
729 self._ret_exponent = 1.0 / self.a - 1.0
730 self._cfr_exponent = -self.a / (self.a - 1.0)
732 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
733 """``C_T(C) = Δθ · (C/K_s)^(1/a)``. Returns 0 at C=0 (no clamp)."""
734 if not isinstance(c, np.ndarray):
735 cf = float(c)
736 c_eff = max(0.0, cf)
737 return self.delta_theta * (c_eff / self.k_s) ** self._inv_a
738 c_arr = np.maximum(np.asarray(c, dtype=float), 0.0)
739 return self.delta_theta * (c_arr / self.k_s) ** self._inv_a
741 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
742 """``R(C) = (Δθ / (a·K_s)) · (C/K_s)^(1/a − 1)``. Clamped at ``_C_MIN``."""
743 if not isinstance(c, np.ndarray):
744 cf = float(c)
745 c_eff = max(_C_MIN, cf)
746 return self._ret_coefficient * (c_eff / self.k_s) ** self._ret_exponent
747 c_eff = np.maximum(np.asarray(c, dtype=float), _C_MIN)
748 return self._ret_coefficient * (c_eff / self.k_s) ** self._ret_exponent
750 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
751 """``C = K_s · (R · a · K_s / Δθ)^{−a/(a−1)}``. Result clamped at ``_C_MIN``."""
752 if not isinstance(r, np.ndarray):
753 base = float(r) * self.a * self.k_s / self.delta_theta
754 if base > 0.0:
755 return max(self.k_s * base**self._cfr_exponent, _C_MIN)
756 return _C_MIN
757 base = np.asarray(r, dtype=float) * self.a * self.k_s / self.delta_theta
758 safe_base = np.where(base > 0, base, 1.0)
759 ratio = safe_base**self._cfr_exponent
760 return np.where(base > 0, np.maximum(self.k_s * ratio, _C_MIN), _C_MIN)
763@dataclass
764class VanGenuchtenMualemConductivity(NonlinearSorption):
765 r"""Mualem prediction for the van Genuchten retention curve, recast as NonlinearSorption.
767 Used by :mod:`gwtransport.percolation` for Kinematic-Wave percolation
768 with the standard Mualem-van Genuchten conductivity curve
770 .. math::
771 K(\\theta) = K_s \\cdot S_e^L \\cdot
772 \\left[1 - \\left(1 - S_e^{1/m}\\right)^m\\right]^2, \\qquad
773 S_e = (\\theta - \\theta_r)/(\\theta_s - \\theta_r), \\qquad
774 m = 1 - 1/n_\\text{vG}.
776 The retention parameter ``α_vG`` is *not* needed for ``K(θ)`` — the
777 Kinematic-Wave approximation drops capillary suction, so only the
778 ``K(S_e)`` curve matters. The two inversions ``S_e(C)`` and
779 ``S_e(R)`` have no closed form; both use ``scipy.optimize.brentq``
780 with ``xtol = BRENTQ_XTOL = 1e-14``.
782 See Also
783 --------
784 BrooksCoreyConductivity : Brooks-Corey closed-form variant.
785 gwtransport.percolation.root_zone_to_water_table_kinematic_wave : The public wrapper.
787 Notes
788 -----
789 The closed-form derivative is
791 .. math::
792 \\frac{dK_M}{dS_e} = K_s \\cdot S_e^{L-1} \\cdot U \\cdot
793 \\left[L \\cdot U + 2 \\cdot S_e^{1/m} \\cdot T^{m-1}\\right],
795 with ``T = 1 - S_e^{1/m}`` and ``U = 1 - T^m``. It is evaluated by
796 ``_dk_dse`` and used for ``retardation(C)`` (after solving ``S_e(C)``)
797 and as the brentq objective in ``_se_from_retardation``.
799 ``dK_M/dS_e`` is strictly increasing for every ``n_vG > 1`` and
800 ``L ≥ 0`` (the conductivity flux is convex, ``d²K/dS_e² > 0``; proved at
801 200-digit precision in ``docs/theory/front_tracking_interactions.md`` §2),
802 so the brentq inversions are always well-posed — no monotonicity guard is
803 needed. A convex flux also means the transport admits only shocks and
804 rarefactions, never compound waves.
806 Examples
807 --------
808 >>> sorption = VanGenuchtenMualemConductivity(
809 ... theta_r=0.01, theta_s=0.337, k_s=0.174, van_genuchten_n=2.28
810 ... )
811 >>> r = sorption.retardation(0.05)
812 >>> c = sorption.concentration_from_retardation(r)
813 >>> bool(np.isclose(c, 0.05, rtol=1e-12))
814 True
815 """
817 theta_r: float
818 """Residual volumetric moisture content [-]; ``0 <= theta_r < theta_s``."""
819 theta_s: float
820 """Saturated volumetric moisture content [-]; ``theta_r < theta_s < 1``."""
821 k_s: float
822 """Saturated hydraulic conductivity [length/time]. Positive."""
823 van_genuchten_n: float
824 """vG shape parameter ``n_vG > 1``; ``m = 1 − 1/n_vG`` is derived from it."""
825 mualem_l: float = 0.5
826 """Mualem pore-connectivity ``L >= 0``. Default 0.5 (standard Mualem).
828 ``L = 0`` (Burdine variant) gives a closed-form ``S_e(C)`` inverse; ``L != 0``
829 requires ``brentq``."""
830 m: float = field(init=False)
831 """Derived ``m = 1 − 1/n_vG``; set in ``__post_init__``."""
832 delta_theta: float = field(init=False)
833 """``θ_s − θ_r``; set in ``__post_init__``."""
835 def __post_init__(self) -> None:
836 """Validate parameters and derive ``m``, ``delta_theta``.
838 No convexity/monotonicity guard is needed: ``dK_M/dS_e`` is strictly
839 increasing (the flux ``K(S_e)`` is convex, ``d²K/dS_e² > 0``) for every
840 ``n_vG > 1`` and ``L ≥ 0`` — proved at 200-digit precision in
841 ``docs/theory/front_tracking_interactions.md`` §2 — so the brentq
842 inversions are always well-posed.
844 Raises
845 ------
846 ValueError
847 If any parameter is outside its valid range.
848 """
849 if not 0.0 <= self.theta_r < self.theta_s:
850 msg = f"theta_r must satisfy 0 <= theta_r < theta_s, got theta_r={self.theta_r}, theta_s={self.theta_s}"
851 raise ValueError(msg)
852 if not self.theta_s < 1.0:
853 msg = f"theta_s must be < 1, got {self.theta_s}"
854 raise ValueError(msg)
855 if self.k_s <= 0.0:
856 msg = f"k_s must be positive, got {self.k_s}"
857 raise ValueError(msg)
858 if self.van_genuchten_n <= 1.0:
859 msg = f"van_genuchten_n must be > 1, got {self.van_genuchten_n}"
860 raise ValueError(msg)
861 if self.mualem_l < 0.0:
862 msg = f"mualem_l must be >= 0, got {self.mualem_l}"
863 raise ValueError(msg)
864 self.m = 1.0 - 1.0 / self.van_genuchten_n
865 self.delta_theta = self.theta_s - self.theta_r
867 def _k_se(self, s: float) -> float:
868 """``K_M(S_e)`` evaluated at a scalar ``S_e``. Returns 0 at ``S_e = 0``."""
869 if s <= 0.0:
870 return 0.0
871 if s >= 1.0:
872 return self.k_s
873 t = 1.0 - s ** (1.0 / self.m)
874 u = 1.0 - t**self.m
875 return self.k_s * s**self.mualem_l * u * u
877 def _dk_dse(self, s: float) -> float:
878 """Closed-form ``dK_M/dS_e`` at scalar ``S_e``.
880 At ``s → 1`` (saturation), ``dK/dS_e`` diverges because ``t^(m-1) → ∞``
881 for ``m < 1``. The function returns ``+∞`` at and above ``s = 1`` so that
882 ``brentq`` can use ``s = 1`` as a closed upper bracket endpoint.
883 """
884 if s <= 0.0:
885 # Limit form: K vanishes as S^(L + 2/m), so derivative is 0 at S=0.
886 return 0.0
887 s_pow_inv_m = s ** (1.0 / self.m)
888 t = 1.0 - s_pow_inv_m
889 if t <= 0.0:
890 # Numerical underflow or s ≥ 1 — the dK/dS_e singularity at saturation.
891 return float("inf")
892 u = 1.0 - t**self.m
893 return self.k_s * s ** (self.mualem_l - 1.0) * u * (self.mualem_l * u + 2.0 * s_pow_inv_m * t ** (self.m - 1.0))
895 def _se_from_c(self, c: float) -> float:
896 """Invert ``K_M(S_e) = c`` for ``S_e``. Closed form for ``mualem_l = 0``; brentq otherwise.
898 For the Burdine variant (``L = 0``), ``K_M(S_e) = K_s · [1 − (1 − S_e^{1/m})^m]^2``
899 is invertible as ``S_e = (1 − (1 − √(K/K_s))^{1/m})^m`` — completely closed
900 form. For ``L ≠ 0`` (default Mualem ``L = 0.5``), no closed-form inverse
901 exists; ``scipy.optimize.brentq`` with ``xtol = BRENTQ_XTOL = 1e-14`` is
902 used. The brentq call is unavoidable in the Mualem case because the
903 ``K_M(S_e)`` function is transcendental.
904 """
905 c_eff = max(float(c), _C_MIN)
906 if c_eff >= self.k_s:
907 return 1.0
908 if self.mualem_l == 0.0:
909 u = (c_eff / self.k_s) ** 0.5 # U = 1 − (1−S_e^{1/m})^m
910 one_minus_u = 1.0 - u
911 one_minus_q_to_inv_m = one_minus_u ** (1.0 / self.m)
912 q = 1.0 - one_minus_q_to_inv_m
913 return float(q**self.m)
914 return float(brentq(lambda s: self._k_se(s) - c_eff, _C_MIN, 1.0, xtol=BRENTQ_XTOL)) # type: ignore[arg-type]
916 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
917 """``C_T = Δθ · S_e(C)``. Returns 0 at C=0 (no clamp)."""
918 is_array = isinstance(c, np.ndarray)
919 c_arr = np.maximum(np.asarray(c, dtype=float), 0.0)
920 flat = c_arr.ravel()
921 se = np.fromiter(
922 (self._se_from_c(ci) if ci > 0.0 else 0.0 for ci in flat), dtype=float, count=flat.size
923 ).reshape(c_arr.shape)
924 result = self.delta_theta * se
925 return result if is_array else float(result)
927 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
928 """``R = Δθ / (dK_M/dS_e)|_{S_e(C)}`` via ``_dk_dse``; clamps C at ``_C_MIN``."""
929 is_array = isinstance(c, np.ndarray)
930 c_arr = np.maximum(np.asarray(c, dtype=float), _C_MIN)
931 flat = c_arr.ravel()
932 out = np.empty(flat.size, dtype=float)
933 for i, ci in enumerate(flat):
934 s = self._se_from_c(ci)
935 out[i] = self.delta_theta / self._dk_dse(s)
936 result = out.reshape(c_arr.shape)
937 return result if is_array else float(result)
939 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
940 """Invert ``R(C) = r``. Solve ``dK_M/dS_e(S_e) = Δθ/r`` via brentq, then ``C = K_M(S_e)``."""
941 is_array = isinstance(r, np.ndarray)
942 r_arr = np.asarray(r, dtype=float)
943 flat = r_arr.ravel()
944 out = np.empty(flat.size, dtype=float)
945 for i, ri in enumerate(flat):
946 s = self._se_from_retardation(float(ri))
947 out[i] = max(self._k_se(s), _C_MIN)
948 result = out.reshape(r_arr.shape)
949 return result if is_array else float(result)
951 def _se_from_retardation(self, r: float) -> float:
952 """Invert ``dK_M/dS_e(S_e) = Δθ/r`` for ``S_e`` via brentq.
954 Single root-find for vG-Mualem; shared by ``concentration_from_retardation``
955 and ``c_and_total_from_retardation`` to avoid duplicate brentq calls.
956 """
957 if r <= 0.0:
958 return _C_MIN
959 target = self.delta_theta / r
960 try:
961 return float(brentq(lambda s, tgt=target: self._dk_dse(s) - tgt, _C_MIN, 1.0, xtol=BRENTQ_XTOL)) # type: ignore[arg-type]
962 except ValueError:
963 return _C_MIN
965 def c_and_total_from_retardation(self, r: float) -> tuple[float, float]:
966 """Return ``(c, C_T)`` at retardation ``r`` from a SINGLE brentq call.
968 Overrides the default base-class implementation (which calls
969 ``concentration_from_retardation`` and ``total_concentration``
970 separately and ends up doing two independent brentq solves on the same
971 underlying equation). Halves the iterative-solver cost in the IBP fan
972 integrators.
973 """
974 s = self._se_from_retardation(r)
975 c = max(self._k_se(s), _C_MIN)
976 ct = self.delta_theta * s
977 return c, ct
980SorptionModel = NonlinearSorption | ConstantRetardation
981"""Type alias for all sorption models accepted by the front-tracking solver."""
984def characteristic_speed(c: float, sorption: SorptionModel) -> float:
985 """Compute characteristic speed dV/dθ = 1/R(C).
987 In (V, θ) coordinates, every characteristic propagates at a flow-free
988 speed determined solely by the local concentration and the sorption
989 isotherm.
991 Parameters
992 ----------
993 c : float
994 Dissolved concentration [mass/volume].
995 sorption : SorptionModel
996 Sorption model.
998 Returns
999 -------
1000 speed : float
1001 Characteristic speed dV/dθ.
1003 Examples
1004 --------
1005 >>> sorption = FreundlichSorption(
1006 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
1007 ... )
1008 >>> s = characteristic_speed(c=5.0, sorption=sorption)
1009 >>> s > 0
1010 True
1011 """
1012 r = float(sorption.retardation(c))
1013 return float("inf") if r == 0.0 else 1.0 / r
1016def characteristic_position(
1017 c: float,
1018 sorption: SorptionModel,
1019 theta_start: float,
1020 v_start: float,
1021 theta: float,
1022) -> float | None:
1023 """Compute position of a characteristic at cumulative flow θ.
1025 Characteristics propagate linearly in θ::
1027 V(θ) = v_start + characteristic_speed(C) * (θ - θ_start)
1029 Parameters
1030 ----------
1031 c : float
1032 Concentration carried by characteristic [mass/volume].
1033 sorption : SorptionModel
1034 Sorption model.
1035 theta_start : float
1036 Cumulative flow at which the characteristic starts [m³].
1037 v_start : float
1038 Starting position [m³].
1039 theta : float
1040 Cumulative flow at which to evaluate position [m³].
1042 Returns
1043 -------
1044 position : float or None
1045 Position at θ [m³], or None if θ < θ_start.
1047 Examples
1048 --------
1049 >>> sorption = ConstantRetardation(retardation_factor=2.0)
1050 >>> v = characteristic_position(
1051 ... c=5.0, sorption=sorption, theta_start=0.0, v_start=0.0, theta=1000.0
1052 ... )
1053 >>> bool(np.isclose(v, 500.0)) # v = (1/2) * 1000 = 500
1054 True
1055 """
1056 if theta < theta_start:
1057 return None
1059 return v_start + characteristic_speed(c, sorption) * (theta - theta_start)
1062def compute_first_front_arrival_theta(
1063 cin: npt.NDArray[np.floating],
1064 theta_edges: npt.NDArray[np.floating],
1065 aquifer_pore_volume: float,
1066 sorption: SorptionModel,
1067) -> float:
1068 """Cumulative-flow θ at which ``c_first`` arrives at the outlet (end of spin-up).
1070 "Arrival" means the θ at which the ``c_first`` *level* is fully present at
1071 the outlet, ``θ_emit + V·R(c_first)`` for ``n<1`` and
1072 ``θ_emit + V·C_T(c_first)/c_first`` for ``n>1``/constant retardation.
1074 .. warning::
1076 For ``n<1`` with ``c_min > 0`` (default ``c_min = 1e-12`` in
1077 :class:`FreundlichSorption`), the actual wave emitted is a
1078 :class:`~gwtransport.fronttracking.waves.RarefactionWave` whose head (``c = c_min ≈ 0``) reaches the
1079 outlet at θ ≈ ``V·R(c_min) ≈ V`` — *much* earlier than the value this
1080 function returns (which is the *tail* arrival ``V·R(c_first)``).
1081 The function returns "tail arrival" semantics: the returned θ is a
1082 conservative end-of-spin-up where c ≤ c_first everywhere before it.
1083 Consult the solver event log for the true rarefaction head crossing.
1085 Parameters
1086 ----------
1087 cin : numpy.ndarray
1088 Inlet concentration [mass/volume].
1089 theta_edges : numpy.ndarray
1090 Cumulative-flow edges; length ``len(cin) + 1``.
1091 aquifer_pore_volume : float
1092 Total pore volume [m³]. Must be positive.
1093 sorption : SorptionModel
1094 Sorption model.
1096 Returns
1097 -------
1098 theta_first_arrival : float
1099 Cumulative-flow θ at which ``c_first`` is fully present at the outlet
1100 [m³]. Returns ``np.inf`` only if no water-carrying (positive θ-width)
1101 bin has nonzero ``cin``.
1103 Examples
1104 --------
1105 >>> cin = np.array([0.0, 10.0] + [10.0] * 10)
1106 >>> theta_edges = np.arange(0.0, 1300.0, 100.0) # constant flow=100, dt=1
1107 >>> sorption = ConstantRetardation(retardation_factor=2.0)
1108 >>> theta_first = compute_first_front_arrival_theta(
1109 ... cin, theta_edges, 500.0, sorption
1110 ... )
1111 >>> bool(np.isclose(theta_first, 100.0 + 500.0 * 2.0)) # θ_emit + V·R
1112 True
1113 """
1114 # Zero-flow (pump-off) bins have zero θ-width and emit no wave (the solver
1115 # skips them), so they cannot be the first arrival either.
1116 nonzero_indices = np.where((cin > 0) & (np.diff(theta_edges) > 0))[0]
1117 if len(nonzero_indices) == 0:
1118 return float(np.inf)
1120 idx_first = int(nonzero_indices[0])
1121 c_first = float(cin[idx_first])
1123 if isinstance(sorption, FreundlichSorption) and sorption.n < 1.0:
1124 # n<1: the 0→c_first step emits a rarefaction; its tail (c=c_first)
1125 # reaches the outlet after V·R(c_first) units of cumulative flow.
1126 target_volume = aquifer_pore_volume * float(sorption.retardation(c_first))
1127 else:
1128 # n>1 or constant: R-H shock with speed = c / (C_T(c) - C_T(0));
1129 # target volume = V · C_T(c_first) / c_first.
1130 target_volume = aquifer_pore_volume * float(sorption.total_concentration(c_first)) / c_first
1132 return float(theta_edges[idx_first]) + target_volume