Coverage for src/gwtransport/fronttracking/math.py: 96%
318 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +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.
18This file is part of gwtransport which is released under AGPL-3.0 license.
19See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
20"""
22from abc import ABC, abstractmethod
23from dataclasses import dataclass, field
25import numpy as np
26import numpy.typing as npt
27from scipy.optimize import brentq
29# Numerical tolerance constants
30EPSILON_FREUNDLICH_N = 1e-10 # Tolerance for checking if n ≈ 1.0 (Freundlich constructor rejects this)
31EPSILON_DENOMINATOR = 1e-15 # Tolerance for near-zero denominators in shock velocity
32_C_MIN = 1e-12 # Shared dry-soil singularity floor for Freundlich n>1, Brooks-Corey, vG-Mualem.
33BRENTQ_XTOL = 1e-14 # brentq absolute tolerance for vG-Mualem inversions; matches _invert_freundlich_cr_zero.
36class NonlinearSorption(ABC):
37 """Abstract base for concentration-dependent sorption models.
39 Subclasses must implement `retardation`, `total_concentration`, and
40 `concentration_from_retardation`. Shock velocity and entropy checking
41 are provided generically via the Rankine-Hugoniot and Lax conditions.
43 See Also
44 --------
45 FreundlichSorption : Freundlich isotherm implementation.
46 LangmuirSorption : Langmuir isotherm implementation.
47 ConstantRetardation : Linear (constant R) retardation model.
48 """
50 @abstractmethod
51 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
52 """Compute retardation factor R(C)."""
54 @abstractmethod
55 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
56 """Compute total concentration (dissolved + sorbed per unit pore volume)."""
58 @abstractmethod
59 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
60 """Invert retardation factor to obtain concentration."""
62 def shock_speed(self, c_left: float, c_right: float) -> float:
63 """Compute shock speed dV/dθ via Rankine-Hugoniot in (V, θ) coordinates.
65 With cumulative-flow coordinate θ = ∫flow(t') dt', the PDE
66 ``∂C_T/∂t + flow·∂C/∂V = 0`` becomes ``∂C_T/∂θ + ∂C/∂V = 0``, and
67 Rankine-Hugoniot reduces to::
69 dV_s/dθ = (C_R - C_L) / (C_T(C_R) - C_T(C_L))
71 Flow drops out entirely; the result is a property of the sorption
72 isotherm alone.
74 Parameters
75 ----------
76 c_left : float
77 Concentration upstream (behind) shock [mass/volume].
78 c_right : float
79 Concentration downstream (ahead of) shock [mass/volume].
81 Returns
82 -------
83 shock_speed : float
84 Shock speed dV/dθ [m³ / m³ flow = dimensionless].
85 """
86 c_total_left = self.total_concentration(c_left)
87 c_total_right = self.total_concentration(c_right)
88 denom = c_total_right - c_total_left
90 if abs(denom) < EPSILON_DENOMINATOR:
91 avg_retardation = 0.5 * float(self.retardation(c_left) + self.retardation(c_right))
92 # Degenerate (zero-strength) shock: its speed is the characteristic speed 1/R. A pair of
93 # saturated states (R = 0, e.g. Mualem-vG at S_e = 1) gives +∞, matching characteristic_speed.
94 return float("inf") if avg_retardation == 0.0 else 1.0 / avg_retardation
96 return float((c_right - c_left) / denom)
98 def c_and_total_from_retardation(self, r: float) -> tuple[float, float]:
99 """Return ``(c, C_T(c))`` at a given retardation ``r``.
101 Default implementation calls ``concentration_from_retardation(r)`` then
102 ``total_concentration(c)`` — two independent root-finds for sorptions
103 where both routes back-solve the same equation (e.g. vG-Mualem with
104 ``L ≠ 0``). Subclasses for which both can be computed from a single
105 root-find should override this for ~2× speedup of the IBP fan
106 integrators.
107 """
108 c = float(self.concentration_from_retardation(r))
109 ct = float(self.total_concentration(c))
110 return c, ct
112 def fan_converges_at_infinity(self) -> bool: # noqa: PLR6301
113 """Whether a ``c_apex=0`` fan's ``∫ c dθ`` converges as ``θ → +∞``.
115 True when ``c → 0`` as ``R → ∞`` (so ``base·c → 0`` faster than ``base → ∞``):
116 Brooks-Corey, van Genuchten-Mualem, Langmuir, and Freundlich ``n > 1``. The
117 only divergent case is Freundlich ``n < 1`` (``c → ∞`` as ``R → ∞``), which
118 overrides this to ``False``. Used by the universal temporal fan integrator to
119 reject a ``+∞`` upper bound when the integral diverges.
120 """
121 return True
123 def check_entropy_condition(self, c_left: float, c_right: float, shock_speed: float) -> bool:
124 """Verify Lax entropy condition in (V, θ) coordinates.
126 In θ-space, characteristic speeds are ``λ_θ(C) = 1 / R(C)``, and the
127 Lax condition for a physical shock is::
129 λ_θ(C_L) ≥ dV_s/dθ ≥ λ_θ(C_R)
131 Parameters
132 ----------
133 c_left : float
134 Concentration upstream of shock [mass/volume].
135 c_right : float
136 Concentration downstream of shock [mass/volume].
137 shock_speed : float
138 Shock speed dV/dθ.
140 Returns
141 -------
142 satisfies : bool
143 True if shock satisfies entropy condition (is physical).
144 """
145 r_left = float(self.retardation(c_left))
146 r_right = float(self.retardation(c_right))
147 lambda_left = float("inf") if r_left == 0.0 else 1.0 / r_left
148 lambda_right = float("inf") if r_right == 0.0 else 1.0 / r_right
150 # A saturated upstream state (λ_left = +∞, e.g. a Mualem-vG wetting front at S_e = 1) is
151 # physical; reject only a non-finite shock speed or downstream characteristic, where the
152 # Lax test itself is ill-posed.
153 if not np.isfinite(shock_speed) or not np.isfinite(lambda_right):
154 return False
156 finite_left = abs(lambda_left) if np.isfinite(lambda_left) else 0.0
157 tolerance = 1e-14 * max(finite_left, abs(lambda_right), abs(shock_speed))
159 return bool((lambda_left > shock_speed - tolerance) and (shock_speed > lambda_right - tolerance))
162@dataclass
163class FreundlichSorption(NonlinearSorption):
164 """
165 Freundlich sorption isotherm with exact analytical methods.
167 The Freundlich isotherm is: s(C) = k_f * C^(1/n)
169 where:
170 - s is sorbed concentration [mass/mass of solid]
171 - C is dissolved concentration [mass/volume of water]
172 - k_f is Freundlich coefficient [(volume/mass)^(1/n)]
173 - n is Freundlich exponent (dimensionless)
175 For n > 1: Higher C travels faster
176 For n < 1: Higher C travels slower
177 For n = 1: linear (not supported, use ConstantRetardation instead)
179 Parameters
180 ----------
181 k_f : float
182 Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive.
183 n : float
184 Freundlich exponent [-]. Must be positive and != 1.
185 bulk_density : float
186 Bulk density of porous medium [kg/m³]. Must be positive.
187 porosity : float
188 Porosity [-]. Must be in (0, 1).
189 c_min : float, optional
190 Minimum concentration threshold (the dry-soil singularity floor). For
191 n>1, prevents infinite retardation as C→0. Default ``1e-12`` for all n.
193 Notes
194 -----
195 The retardation factor is defined as:
196 R(C) = 1 + (rho_b/n_por) * ds/dC
197 = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)
199 For Freundlich sorption, R depends on C, which creates nonlinear wave behavior.
201 For n>1 (higher C travels faster), R(C)→∞ as C→0, which can cause extremely slow
202 wave propagation. The c_min parameter prevents this by enforcing a minimum
203 concentration, making R(C) finite for all C≥0.
205 Examples
206 --------
207 >>> sorption = FreundlichSorption(
208 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
209 ... )
210 >>> r = sorption.retardation(5.0)
211 >>> c_back = sorption.concentration_from_retardation(r)
212 >>> bool(np.isclose(c_back, 5.0))
213 True
214 """
216 k_f: float
217 """Freundlich coefficient [(m³/kg)^(1/n)]."""
218 n: float
219 """Freundlich exponent [-]."""
220 bulk_density: float
221 """Bulk density of porous medium [kg/m³]."""
222 porosity: float
223 """Porosity [-]."""
224 c_min: float = 1e-12
225 """Minimum concentration threshold to prevent infinite retardation."""
227 def __post_init__(self):
228 """Validate parameters after initialization.
230 Raises
231 ------
232 ValueError
233 If any parameter is outside its valid range: ``k_f`` <= 0,
234 ``n`` <= 0, ``n`` == 1, ``bulk_density`` <= 0, ``porosity``
235 outside (0, 1), or ``c_min`` < 0.
236 """
237 if self.k_f <= 0:
238 msg = f"k_f must be positive, got {self.k_f}"
239 raise ValueError(msg)
240 if self.n <= 0:
241 msg = f"n must be positive, got {self.n}"
242 raise ValueError(msg)
243 if abs(self.n - 1.0) < EPSILON_FREUNDLICH_N:
244 msg = "n = 1 (linear case) not supported, use ConstantRetardation instead"
245 raise ValueError(msg)
246 if self.bulk_density <= 0:
247 msg = f"bulk_density must be positive, got {self.bulk_density}"
248 raise ValueError(msg)
249 if not 0 < self.porosity < 1:
250 msg = f"porosity must be in (0, 1), got {self.porosity}"
251 raise ValueError(msg)
252 if self.c_min < 0:
253 msg = f"c_min must be non-negative, got {self.c_min}"
254 raise ValueError(msg)
256 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
257 """
258 Compute retardation factor R(C).
260 The retardation factor relates concentration speed to pore water speed in
261 (V, θ) coordinates::
263 dV/dθ = 1 / R(C)
265 For Freundlich sorption::
267 R(C) = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)
269 Parameters
270 ----------
271 c : float or array-like
272 Dissolved concentration [mass/volume]. Non-negative.
274 Returns
275 -------
276 r : float or numpy.ndarray
277 Retardation factor [-]. Always >= 1.0.
279 Notes
280 -----
281 - For n > 1: R decreases with increasing C (higher C travels faster)
282 - For n < 1: R increases with increasing C (higher C travels slower)
283 - n<1 with c_min=0: R(0)=1 (no sorption at zero, physically correct)
284 because clamping to ``c_min=0`` leaves ``C^((1/n)-1) = 0^positive = 0``.
285 - Otherwise: ``c`` is clamped to ``c_min`` before evaluation. This pairs with
286 :meth:`total_concentration`, which also clamps to ``c_min``.
288 Clamping with ``np.maximum`` before the power keeps a single general path
289 for every ``(n, c_min)`` combination and avoids raising the base to a
290 fractional power on negative ``c``.
291 """
292 is_array = isinstance(c, np.ndarray)
293 c_eff = np.maximum(np.asarray(c), self.c_min)
294 exponent = (1.0 / self.n) - 1.0
295 coefficient = (self.bulk_density * self.k_f) / (self.porosity * self.n)
296 result = 1.0 + coefficient * (c_eff**exponent)
297 return result if is_array else float(result)
299 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
300 """
301 Compute total concentration (dissolved + sorbed per unit pore volume).
303 Total concentration includes both dissolved and sorbed mass:
304 C_total = C + (rho_b/n_por) * s(C)
305 = C + (rho_b/n_por) * k_f * C^(1/n)
307 Parameters
308 ----------
309 c : float or array-like
310 Dissolved concentration [mass/volume]. Non-negative.
312 Returns
313 -------
314 c_total : float or numpy.ndarray
315 Total concentration [mass/volume]. Always >= c.
317 Notes
318 -----
319 This is the conserved quantity in the transport equation:
320 ∂C_total/∂t + ∂(flow*C)/∂v = 0
322 The flux term only includes dissolved concentration because sorbed mass
323 is immobile.
325 For ``c = 0``, ``c^(1/n) = 0`` exactly (no singularity for any
326 ``n > 0``), so ``C_T(0) = 0`` is physically correct and no ``c_min``
327 clamp is needed here. ``c_min`` is only required to keep
328 :meth:`retardation` finite as ``c -> 0`` for ``n > 1``; clamping
329 ``total_concentration`` to ``c_min`` would bias Rankine-Hugoniot
330 shock speeds when ``c_R = 0`` (e.g. the canonical 0->c->0 pulse).
331 Negative ``c`` is clamped to ``0`` defensively.
332 """
333 is_array = isinstance(c, np.ndarray)
334 c_arr = np.maximum(np.asarray(c), 0.0)
335 sorbed = (self.bulk_density / self.porosity) * self.k_f * (c_arr ** (1.0 / self.n))
336 result = c_arr + sorbed
337 return result if is_array else float(result)
339 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
340 """
341 Invert retardation factor to obtain concentration analytically.
343 Given R, solves R = retardation(C) for C. This is used in rarefaction waves
344 where the self-similar solution gives R as a function of position and time.
346 Parameters
347 ----------
348 r : float or array-like
349 Retardation factor [-]. Must be >= 1.0.
351 Returns
352 -------
353 c : float or numpy.ndarray
354 Dissolved concentration [mass/volume]. Non-negative.
356 Notes
357 -----
358 This inverts the relation:
359 R = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)
361 The analytical solution is:
362 C = [(R-1) * n_por*n / (rho_b*k_f)]^(n/(1-n))
364 For n = 1 (linear sorption), the exponent n/(1-n) is undefined, which is
365 why linear sorption must use ConstantRetardation class instead.
367 Examples
368 --------
369 >>> sorption = FreundlichSorption(
370 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
371 ... )
372 >>> r = sorption.retardation(5.0)
373 >>> c = sorption.concentration_from_retardation(r)
374 >>> bool(np.isclose(c, 5.0, rtol=1e-14))
375 True
376 """
377 is_array = isinstance(r, np.ndarray)
378 r_arr = np.asarray(r)
380 # FreundlichSorption.__post_init__ rejects |n-1| < EPSILON_FREUNDLICH_N,
381 # so the previous n≈1 guard here was unreachable.
382 exponent = (1.0 / self.n) - 1.0
383 coefficient = (self.bulk_density * self.k_f) / (self.porosity * self.n)
384 base = (r_arr - 1.0) / coefficient
385 inversion_exponent = 1.0 / exponent
387 # Mask base to a safe placeholder before exponentiation; NumPy emits
388 # RuntimeWarning otherwise for base <= 0 with a fractional exponent.
389 safe_base = np.where(base > 0, base, 1.0)
390 c = safe_base**inversion_exponent
391 result = np.where(base > 0, np.maximum(c, self.c_min), self.c_min)
393 return result if is_array else float(result)
395 def fan_converges_at_infinity(self) -> bool:
396 """Freundlich ``n > 1``: ``c → 0`` as ``R → ∞`` (converges). ``n < 1``: ``c → ∞`` (diverges)."""
397 return self.n > 1.0
400@dataclass
401class ConstantRetardation:
402 """
403 Constant (linear) retardation model.
405 For linear sorption: s(C) = K_d * C
406 This gives constant retardation: R(C) = 1 + (rho_b/n_por) * K_d = constant
408 This is a special case where concentration-dependent behavior disappears.
409 Used for conservative tracers or as approximation for weak sorption.
411 Parameters
412 ----------
413 retardation_factor : float
414 Constant retardation factor [-]. Must be >= 1.0.
415 R = 1.0 means no retardation (conservative tracer).
417 Notes
418 -----
419 With constant retardation:
420 - All concentrations travel at same speed in (V, θ): dV/dθ = 1/R
421 - No rarefaction waves form (all concentrations travel together)
422 - Shocks occur only at concentration discontinuities at inlet
423 - Solution reduces to simple θ-shifting (and then t-shifting via the θ↔t map)
425 This is equivalent to a single-pore-volume advective time-shift (the deterministic limit of
426 :func:`gwtransport.advection.infiltration_to_extraction`) in the gwtransport package.
428 Examples
429 --------
430 >>> sorption = ConstantRetardation(retardation_factor=2.0)
431 >>> sorption.retardation(5.0)
432 2.0
433 >>> sorption.retardation(10.0)
434 2.0
435 """
437 retardation_factor: float
438 """Constant retardation factor [-]. Must be >= 1.0."""
440 def __post_init__(self):
441 """Validate parameters after initialization.
443 Raises
444 ------
445 ValueError
446 If ``retardation_factor`` is less than 1.0.
447 """
448 if self.retardation_factor < 1.0:
449 msg = f"retardation_factor must be >= 1.0, got {self.retardation_factor}"
450 raise ValueError(msg)
452 def retardation(self, c: float) -> float: # noqa: ARG002
453 """
454 Return constant retardation factor (independent of concentration).
456 Parameters
457 ----------
458 c : float
459 Dissolved concentration (not used for constant retardation).
461 Returns
462 -------
463 r : float
464 Constant retardation factor.
465 """
466 return self.retardation_factor
468 def total_concentration(self, c: float) -> float:
469 """
470 Compute total concentration for linear sorption.
472 For constant retardation:
473 C_total = C * R
475 Parameters
476 ----------
477 c : float
478 Dissolved concentration [mass/volume].
480 Returns
481 -------
482 c_total : float
483 Total concentration [mass/volume].
484 """
485 return c * self.retardation_factor
487 def concentration_from_retardation(self, r: float) -> float:
488 """
489 Not applicable for constant retardation.
491 With constant R, all concentrations have the same retardation, so
492 inversion is not meaningful. This method raises an error.
494 Raises
495 ------
496 NotImplementedError
497 Always raised for constant retardation.
498 """
499 msg = "concentration_from_retardation not applicable for ConstantRetardation (R is independent of C)"
500 raise NotImplementedError(msg)
502 def shock_speed(self, c_left: float, c_right: float) -> float: # noqa: ARG002
503 """Compute shock speed dV/dθ for constant retardation.
505 With constant R, ``dV/dθ = 1 / R`` for any concentration pair —
506 identical to every characteristic speed.
508 Parameters
509 ----------
510 c_left, c_right : float
511 Concentrations (unused — kept for ABC compatibility).
513 Returns
514 -------
515 shock_speed : float
516 Shock speed dV/dθ = 1/R.
517 """
518 return 1.0 / self.retardation_factor
520 def check_entropy_condition(self, c_left: float, c_right: float, shock_speed: float) -> bool: # noqa: PLR6301
521 """Entropy condition for constant retardation: trivially satisfied.
523 With constant R every characteristic speed equals the shock speed in
524 θ-space, so the Lax condition holds as an equality regardless of
525 ``c_left``/``c_right``.
527 Returns
528 -------
529 satisfies : bool
530 Always True.
531 """
532 del c_left, c_right, shock_speed
533 return True
536@dataclass
537class LangmuirSorption(NonlinearSorption):
538 """
539 Langmuir sorption isotherm with exact analytical methods.
541 The Langmuir isotherm is: s(C) = s_max * C / (K_L + C)
543 where:
544 - s is sorbed concentration [mass/mass of solid]
545 - C is dissolved concentration [mass/volume of water]
546 - s_max is maximum sorption capacity [mass/mass of solid]
547 - K_L is half-saturation constant [mass/volume]
549 Retardation always decreases with C (favorable isotherm), and R(0) is
550 finite — unlike Freundlich with n > 1, no minimum concentration threshold
551 is needed.
553 Parameters
554 ----------
555 s_max : float
556 Maximum sorption capacity [mass/mass of solid]. Must be positive.
557 k_l : float
558 Half-saturation constant [mass/volume]. Concentration at which
559 s = s_max / 2. Must be positive.
560 bulk_density : float
561 Bulk density of porous medium [kg/m³]. Must be positive.
562 porosity : float
563 Porosity [-]. Must be in (0, 1).
565 See Also
566 --------
567 FreundlichSorption : Freundlich isotherm (unbounded sorption).
568 ConstantRetardation : Linear (constant R) retardation model.
569 :ref:`concept-nonlinear-sorption` : Background on nonlinear sorption.
571 Notes
572 -----
573 The retardation factor is defined as:
574 R(C) = 1 + (rho_b * s_max * K_L) / (n_por * (K_L + C)^2)
576 Key properties:
578 - R(0) = 1 + rho_b * s_max / (n_por * K_L) -- finite for all parameters
579 - R -> 1 as C -> infinity (all sorption sites saturated)
580 - R always decreases with increasing C (higher C travels faster)
581 - Shocks form on concentration increases, rarefaction fans on decreases
583 Examples
584 --------
585 >>> sorption = LangmuirSorption(
586 ... s_max=0.1, k_l=5.0, bulk_density=1500.0, porosity=0.3
587 ... )
588 >>> r = sorption.retardation(5.0)
589 >>> c_back = sorption.concentration_from_retardation(r)
590 >>> bool(np.isclose(c_back, 5.0))
591 True
592 """
594 s_max: float
595 """Maximum sorption capacity [mass/mass of solid]."""
596 k_l: float
597 """Half-saturation constant [mass/volume]."""
598 bulk_density: float
599 """Bulk density of porous medium [kg/m³]."""
600 porosity: float
601 """Porosity [-]."""
603 def __post_init__(self):
604 """Validate parameters after initialization.
606 Raises
607 ------
608 ValueError
609 If any parameter is outside its valid range: ``s_max`` <= 0,
610 ``k_l`` <= 0, ``bulk_density`` <= 0, or ``porosity``
611 outside (0, 1).
612 """
613 if self.s_max <= 0:
614 msg = f"s_max must be positive, got {self.s_max}"
615 raise ValueError(msg)
616 if self.k_l <= 0:
617 msg = f"k_l must be positive, got {self.k_l}"
618 raise ValueError(msg)
619 if self.bulk_density <= 0:
620 msg = f"bulk_density must be positive, got {self.bulk_density}"
621 raise ValueError(msg)
622 if not 0 < self.porosity < 1:
623 msg = f"porosity must be in (0, 1), got {self.porosity}"
624 raise ValueError(msg)
626 self.a_coeff: float = self.bulk_density * self.s_max * self.k_l / self.porosity
627 """Lumped retardation constant rho_b * s_max * K_L / n_por."""
629 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
630 """
631 Compute retardation factor R(C).
633 For Langmuir sorption:
634 R(C) = 1 + A / (K_L + C)²
636 where A = rho_b * s_max * K_L / n_por.
638 Parameters
639 ----------
640 c : float or array-like
641 Dissolved concentration [mass/volume]. Non-negative.
643 Returns
644 -------
645 r : float or numpy.ndarray
646 Retardation factor [-]. Always >= 1.0.
648 Notes
649 -----
650 - R(0) = 1 + rho_b * s_max / (n_por * K_L) — always finite
651 - R decreases with increasing C (higher C travels faster)
652 - R → 1 as C → ∞ (all sorption sites saturated)
653 """
654 is_array = isinstance(c, np.ndarray)
655 c_arr = np.asarray(c)
656 c_eff = np.maximum(c_arr, 0.0)
657 result = 1.0 + self.a_coeff / (self.k_l + c_eff) ** 2
658 return result if is_array else float(result)
660 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
661 """
662 Compute total concentration (dissolved + sorbed per unit pore volume).
664 For Langmuir sorption:
665 C_total = C + (rho_b / n_por) * s_max * C / (K_L + C)
667 Parameters
668 ----------
669 c : float or array-like
670 Dissolved concentration [mass/volume]. Non-negative.
672 Returns
673 -------
674 c_total : float or numpy.ndarray
675 Total concentration [mass/volume]. Always >= c.
676 """
677 is_array = isinstance(c, np.ndarray)
678 c_arr = np.asarray(c)
679 c_eff = np.maximum(c_arr, 0.0)
680 sorbed = (self.bulk_density / self.porosity) * self.s_max * c_eff / (self.k_l + c_eff)
681 result = c_arr + sorbed
682 return result if is_array else float(result)
684 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
685 """
686 Invert retardation factor to obtain concentration analytically.
688 Given R, solves R = 1 + A / (K_L + C)² for C:
689 C = sqrt(A / (R - 1)) - K_L
691 Parameters
692 ----------
693 r : float or array-like
694 Retardation factor [-]. Must be >= 1.0.
696 Returns
697 -------
698 c : float or numpy.ndarray
699 Dissolved concentration [mass/volume]. Non-negative.
701 Notes
702 -----
703 For R <= 1, returns 0.0 (unphysical region).
704 For R >= R(0) = 1 + A/K_L², returns 0.0 (at or below zero concentration).
706 Examples
707 --------
708 >>> sorption = LangmuirSorption(
709 ... s_max=0.1, k_l=5.0, bulk_density=1500.0, porosity=0.3
710 ... )
711 >>> r = sorption.retardation(5.0)
712 >>> c = sorption.concentration_from_retardation(r)
713 >>> bool(np.isclose(c, 5.0, rtol=1e-14))
714 True
715 """
716 is_array = isinstance(r, np.ndarray)
717 r_arr = np.asarray(r)
719 r_minus_1 = r_arr - 1.0
720 # Mask r_minus_1 to a safe placeholder before division to avoid the
721 # RuntimeWarning emitted by np.where's eager evaluation when r == 1.
722 safe_r_minus_1 = np.where(r_minus_1 > 0, r_minus_1, 1.0)
723 c = np.where(r_minus_1 > 0, np.sqrt(self.a_coeff / safe_r_minus_1) - self.k_l, 0.0)
724 result = np.maximum(c, 0.0)
726 return result if is_array else float(result)
729@dataclass
730class BrooksCoreyConductivity(NonlinearSorption):
731 r"""Brooks-Corey unsaturated conductivity recast as a NonlinearSorption.
733 Used by :mod:`gwtransport.percolation` to model gravity-driven percolation
734 through a thick unsaturated zone via the Kinematic-Wave method. The
735 closed-form conductivity curve
737 .. math::
738 K(\\theta) = K_s \\cdot \\Theta^a, \\qquad
739 \\Theta = (\\theta - \\theta_r)/(\\theta_s - \\theta_r), \\qquad
740 a = 3 + 2/\\lambda \\;(\\text{Burdine})
742 is recast in the framework's ``(C, C_T)`` variables by identifying
743 ``C ≡ K`` (the flux variable) and ``C_T ≡ θ - θ_r`` (the conserved
744 storage). All three abstract methods have closed forms; ``shock_speed``
745 and ``check_entropy_condition`` are inherited unchanged from
746 :class:`NonlinearSorption`.
748 Parameters
749 ----------
750 theta_r : float
751 Residual volumetric moisture content [-]. Must satisfy
752 ``0 <= theta_r < theta_s``.
753 theta_s : float
754 Saturated volumetric moisture content [-]. Equal to the porosity
755 for typical soils. Must satisfy ``theta_r < theta_s < 1``.
756 k_s : float
757 Saturated hydraulic conductivity [length/time]. Positive.
758 brooks_corey_lambda : float
759 Pore-size distribution index ``λ`` [-]. Positive. The exponent
760 ``a = 3 + 2/λ`` is the Burdine pore-connectivity result. The Mualem
761 variant (``L = 0.5``) gives ``a = 2.5 + 2/λ`` and is not implemented;
762 a user wanting it can re-derive ``λ`` so the Burdine ``a`` matches the
763 desired Mualem exponent.
765 See Also
766 --------
767 VanGenuchtenMualemConductivity : Van Genuchten variant with brentq inversions.
768 FreundlichSorption : Power-law sorption isotherm (closed form, analogous shape).
769 gwtransport.percolation.root_zone_to_water_table_kinematic_wave : The public wrapper.
771 Notes
772 -----
773 The retardation factor and total-concentration relation are:
775 .. math::
776 C_T(C) = \\Delta\\theta \\cdot (C/K_s)^{1/a}, \\qquad
777 R(C) = (\\Delta\\theta / (a K_s)) \\cdot (C/K_s)^{1/a - 1},
779 with ``Δθ = θ_s − θ_r``. Since ``1/a − 1 < 0`` always (``a > 3``),
780 ``R(C) → ∞`` as ``C → 0`` (dry-soil singularity). The class clamps ``C``
781 to a small floor in ``retardation`` and ``concentration_from_retardation``
782 (the same pattern as :class:`FreundlichSorption` with ``n > 1``);
783 ``total_concentration`` and the inherited ``shock_speed`` do **not**
784 clamp, so the canonical wetting-front shock ``c_R = 0`` produces the
785 correct Rankine-Hugoniot velocity.
787 Examples
788 --------
789 >>> sorption = BrooksCoreyConductivity(
790 ... theta_r=0.01, theta_s=0.337, k_s=0.174, brooks_corey_lambda=0.25
791 ... )
792 >>> r = sorption.retardation(0.05)
793 >>> c = sorption.concentration_from_retardation(r)
794 >>> bool(np.isclose(c, 0.05, rtol=1e-13))
795 True
796 """
798 theta_r: float
799 """Residual volumetric moisture content [-]."""
800 theta_s: float
801 """Saturated volumetric moisture content [-]."""
802 k_s: float
803 """Saturated hydraulic conductivity [length/time]."""
804 brooks_corey_lambda: float
805 """Pore-size distribution index λ [-]."""
806 a: float = field(init=False)
807 """Exponent ``a = 3 + 2/λ`` (Burdine); set in ``__post_init__``."""
808 delta_theta: float = field(init=False)
809 """``θ_s − θ_r``; set in ``__post_init__``."""
811 def __post_init__(self) -> None:
812 """Validate parameters and derive ``a``, ``delta_theta``.
814 Raises
815 ------
816 ValueError
817 If any parameter is outside its valid range.
818 """
819 if not 0.0 <= self.theta_r < self.theta_s:
820 msg = f"theta_r must satisfy 0 <= theta_r < theta_s, got theta_r={self.theta_r}, theta_s={self.theta_s}"
821 raise ValueError(msg)
822 if not self.theta_s < 1.0:
823 msg = f"theta_s must be < 1, got {self.theta_s}"
824 raise ValueError(msg)
825 if self.k_s <= 0.0:
826 msg = f"k_s must be positive, got {self.k_s}"
827 raise ValueError(msg)
828 if self.brooks_corey_lambda <= 0.0:
829 msg = f"brooks_corey_lambda must be positive, got {self.brooks_corey_lambda}"
830 raise ValueError(msg)
831 self.a = 3.0 + 2.0 / self.brooks_corey_lambda
832 self.delta_theta = self.theta_s - self.theta_r
834 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
835 """``C_T(C) = Δθ · (C/K_s)^(1/a)``. Returns 0 at C=0 (no clamp)."""
836 is_array = isinstance(c, np.ndarray)
837 c_arr = np.maximum(np.asarray(c, dtype=float), 0.0)
838 result = self.delta_theta * (c_arr / self.k_s) ** (1.0 / self.a)
839 return result if is_array else float(result)
841 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
842 """``R(C) = (Δθ / (a·K_s)) · (C/K_s)^(1/a − 1)``. Clamped at ``_C_MIN``."""
843 is_array = isinstance(c, np.ndarray)
844 c_eff = np.maximum(np.asarray(c, dtype=float), _C_MIN)
845 result = (self.delta_theta / (self.a * self.k_s)) * (c_eff / self.k_s) ** (1.0 / self.a - 1.0)
846 return result if is_array else float(result)
848 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
849 """``C = K_s · (R · a · K_s / Δθ)^{−a/(a−1)}``. Result clamped at ``_C_MIN``."""
850 is_array = isinstance(r, np.ndarray)
851 r_arr = np.asarray(r, dtype=float)
852 base = r_arr * self.a * self.k_s / self.delta_theta
853 safe_base = np.where(base > 0, base, 1.0)
854 ratio = safe_base ** (-self.a / (self.a - 1.0))
855 c = self.k_s * ratio
856 result = np.where(base > 0, np.maximum(c, _C_MIN), _C_MIN)
857 return result if is_array else float(result)
860@dataclass
861class VanGenuchtenMualemConductivity(NonlinearSorption):
862 r"""Mualem prediction for the van Genuchten retention curve, recast as NonlinearSorption.
864 Used by :mod:`gwtransport.percolation` for Kinematic-Wave percolation
865 with the standard Mualem-van Genuchten conductivity curve
867 .. math::
868 K(\\theta) = K_s \\cdot S_e^L \\cdot
869 \\left[1 - \\left(1 - S_e^{1/m}\\right)^m\\right]^2, \\qquad
870 S_e = (\\theta - \\theta_r)/(\\theta_s - \\theta_r), \\qquad
871 m = 1 - 1/n_\\text{vG}.
873 The retention parameter ``α_vG`` is *not* needed for ``K(θ)`` — the
874 Kinematic-Wave approximation drops capillary suction, so only the
875 ``K(S_e)`` curve matters. The two inversions ``S_e(C)`` and
876 ``S_e(R)`` have no closed form; both use ``scipy.optimize.brentq``
877 with ``xtol = BRENTQ_XTOL = 1e-14``.
879 Parameters
880 ----------
881 theta_r : float
882 Residual volumetric moisture content [-].
883 theta_s : float
884 Saturated volumetric moisture content [-].
885 k_s : float
886 Saturated hydraulic conductivity [length/time].
887 van_genuchten_n : float
888 Shape parameter ``n_vG > 1``. ``m = 1 − 1/n_vG`` is derived.
889 mualem_l : float, optional
890 Pore-connectivity parameter ``L``. Default 0.5 (standard Mualem).
891 Must satisfy ``L >= 0``. Setting ``L = 0`` (Burdine variant) gives
892 a closed-form ``S_e(C)`` inverse; ``L != 0`` requires ``brentq``.
894 See Also
895 --------
896 BrooksCoreyConductivity : Brooks-Corey closed-form variant.
897 gwtransport.percolation.root_zone_to_water_table_kinematic_wave : The public wrapper.
899 Notes
900 -----
901 The closed-form derivative is
903 .. math::
904 \\frac{dK_M}{dS_e} = K_s \\cdot S_e^{L-1} \\cdot U \\cdot
905 \\left[L \\cdot U + 2 \\cdot S_e^{1/m} \\cdot T^{m-1}\\right],
907 with ``T = 1 - S_e^{1/m}`` and ``U = 1 - T^m``. Used for
908 ``retardation(C)`` (after solving ``S_e(C)``) and for the brentq
909 objective in ``concentration_from_retardation(R)``. The formula is
910 inlined at both call sites, not exposed as a separate method.
912 The class checks monotonicity of ``dK_M/dS_e`` at a single pair of
913 sample points in ``__post_init__`` (cheap directional check). Truly
914 pathological parameter combinations that yield a non-monotone curve
915 surface as a ``brentq`` ValueError at the first inversion call.
917 Examples
918 --------
919 >>> sorption = VanGenuchtenMualemConductivity(
920 ... theta_r=0.01, theta_s=0.337, k_s=0.174, van_genuchten_n=2.28
921 ... )
922 >>> r = sorption.retardation(0.05)
923 >>> c = sorption.concentration_from_retardation(r)
924 >>> bool(np.isclose(c, 0.05, rtol=1e-12))
925 True
926 """
928 theta_r: float
929 """Residual volumetric moisture content [-]."""
930 theta_s: float
931 """Saturated volumetric moisture content [-]."""
932 k_s: float
933 """Saturated hydraulic conductivity [length/time]."""
934 van_genuchten_n: float
935 """vG shape parameter ``n_vG > 1``."""
936 mualem_l: float = 0.5
937 """Mualem pore-connectivity ``L``. Default 0.5."""
938 m: float = field(init=False)
939 """Derived ``m = 1 − 1/n_vG``; set in ``__post_init__``."""
940 delta_theta: float = field(init=False)
941 """``θ_s − θ_r``; set in ``__post_init__``."""
943 def __post_init__(self) -> None:
944 """Validate parameters and run a single-sample monotonicity check.
946 Raises
947 ------
948 ValueError
949 If parameters are outside their valid range, or if the cheap
950 monotonicity sample at ``S_e = 0.5`` vs ``0.99`` indicates
951 ``dK_M/dS_e`` is non-monotone (pathological).
952 """
953 if not 0.0 <= self.theta_r < self.theta_s:
954 msg = f"theta_r must satisfy 0 <= theta_r < theta_s, got theta_r={self.theta_r}, theta_s={self.theta_s}"
955 raise ValueError(msg)
956 if not self.theta_s < 1.0:
957 msg = f"theta_s must be < 1, got {self.theta_s}"
958 raise ValueError(msg)
959 if self.k_s <= 0.0:
960 msg = f"k_s must be positive, got {self.k_s}"
961 raise ValueError(msg)
962 if self.van_genuchten_n <= 1.0:
963 msg = f"van_genuchten_n must be > 1, got {self.van_genuchten_n}"
964 raise ValueError(msg)
965 if self.mualem_l < 0.0:
966 msg = f"mualem_l must be >= 0, got {self.mualem_l}"
967 raise ValueError(msg)
968 self.m = 1.0 - 1.0 / self.van_genuchten_n
969 self.delta_theta = self.theta_s - self.theta_r
970 # Cheap monotonicity sanity check on dK_M/dS_e.
971 s_low, s_high = 0.5, 0.99
972 if self._dk_dse(s_low) >= self._dk_dse(s_high):
973 msg = (
974 f"Non-monotone dK_M/dS_e detected at the sanity-check samples for "
975 f"van_genuchten_n={self.van_genuchten_n}, mualem_l={self.mualem_l}: "
976 f"dK_M/dS_e({s_low})={self._dk_dse(s_low):.6g} should be < "
977 f"dK_M/dS_e({s_high})={self._dk_dse(s_high):.6g}. "
978 f"Brentq inversions in this class assume monotone-increasing dK_M/dS_e."
979 )
980 raise ValueError(msg)
982 def _k_se(self, s: float) -> float:
983 """``K_M(S_e)`` evaluated at a scalar ``S_e``. Returns 0 at ``S_e = 0``."""
984 if s <= 0.0:
985 return 0.0
986 if s >= 1.0:
987 return self.k_s
988 t = 1.0 - s ** (1.0 / self.m)
989 u = 1.0 - t**self.m
990 return self.k_s * s**self.mualem_l * u * u
992 def _dk_dse(self, s: float) -> float:
993 """Closed-form ``dK_M/dS_e`` at scalar ``S_e``. Inlined at call sites.
995 At ``s → 1`` (saturation), ``dK/dS_e`` diverges because ``t^(m-1) → ∞``
996 for ``m < 1``. The function returns ``+∞`` at and above ``s = 1`` so that
997 ``brentq`` can use ``s = 1`` as a closed upper bracket endpoint.
998 """
999 if s <= 0.0:
1000 # Limit form: K vanishes as S^(L + 2/m), so derivative is 0 at S=0.
1001 return 0.0
1002 s_pow_inv_m = s ** (1.0 / self.m)
1003 t = 1.0 - s_pow_inv_m
1004 if t <= 0.0:
1005 # Numerical underflow or s ≥ 1 — the dK/dS_e singularity at saturation.
1006 return float("inf")
1007 u = 1.0 - t**self.m
1008 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))
1010 def _se_from_c(self, c: float) -> float:
1011 """Invert ``K_M(S_e) = c`` for ``S_e``. Closed form for ``mualem_l = 0``; brentq otherwise.
1013 For the Burdine variant (``L = 0``), ``K_M(S_e) = K_s · [1 − (1 − S_e^{1/m})^m]^2``
1014 is invertible as ``S_e = (1 − (1 − √(K/K_s))^{1/m})^m`` — completely closed
1015 form. For ``L ≠ 0`` (default Mualem ``L = 0.5``), no closed-form inverse
1016 exists; ``scipy.optimize.brentq`` with ``xtol = BRENTQ_XTOL = 1e-14`` is
1017 used. The brentq call is unavoidable in the Mualem case because the
1018 ``K_M(S_e)`` function is transcendental.
1019 """
1020 c_eff = max(float(c), _C_MIN)
1021 if c_eff >= self.k_s:
1022 return 1.0
1023 if self.mualem_l == 0.0:
1024 u = (c_eff / self.k_s) ** 0.5 # U = 1 − (1−S_e^{1/m})^m
1025 one_minus_u = 1.0 - u
1026 one_minus_q_to_inv_m = one_minus_u ** (1.0 / self.m)
1027 q = 1.0 - one_minus_q_to_inv_m
1028 return float(q**self.m)
1029 return float(brentq(lambda s: self._k_se(s) - c_eff, _C_MIN, 1.0, xtol=BRENTQ_XTOL)) # type: ignore[arg-type]
1031 def total_concentration(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
1032 """``C_T = Δθ · S_e(C)``. Returns 0 at C=0 (no clamp)."""
1033 is_array = isinstance(c, np.ndarray)
1034 c_arr = np.maximum(np.asarray(c, dtype=float), 0.0)
1035 flat = c_arr.ravel()
1036 se = np.fromiter(
1037 (self._se_from_c(ci) if ci > 0.0 else 0.0 for ci in flat), dtype=float, count=flat.size
1038 ).reshape(c_arr.shape)
1039 result = self.delta_theta * se
1040 return result if is_array else float(result)
1042 def retardation(self, c: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
1043 """``R = Δθ / (dK_M/dS_e)|_{S_e(C)}``. Uses inlined derivative; clamps C at ``_C_MIN``."""
1044 is_array = isinstance(c, np.ndarray)
1045 c_arr = np.maximum(np.asarray(c, dtype=float), _C_MIN)
1046 flat = c_arr.ravel()
1047 out = np.empty(flat.size, dtype=float)
1048 for i, ci in enumerate(flat):
1049 s = self._se_from_c(ci)
1050 out[i] = self.delta_theta / self._dk_dse(s)
1051 result = out.reshape(c_arr.shape)
1052 return result if is_array else float(result)
1054 def concentration_from_retardation(self, r: float | npt.NDArray[np.float64]) -> float | npt.NDArray[np.float64]:
1055 """Invert ``R(C) = r``. Solve ``dK_M/dS_e(S_e) = Δθ/r`` via brentq, then ``C = K_M(S_e)``."""
1056 is_array = isinstance(r, np.ndarray)
1057 r_arr = np.asarray(r, dtype=float)
1058 flat = r_arr.ravel()
1059 out = np.empty(flat.size, dtype=float)
1060 for i, ri in enumerate(flat):
1061 s = self._se_from_retardation(float(ri))
1062 out[i] = max(self._k_se(s), _C_MIN)
1063 result = out.reshape(r_arr.shape)
1064 return result if is_array else float(result)
1066 def _se_from_retardation(self, r: float) -> float:
1067 """Invert ``dK_M/dS_e(S_e) = Δθ/r`` for ``S_e`` via brentq.
1069 Single root-find for vG-Mualem; shared by ``concentration_from_retardation``
1070 and ``c_and_total_from_retardation`` to avoid duplicate brentq calls.
1071 """
1072 if r <= 0.0:
1073 return _C_MIN
1074 target = self.delta_theta / r
1075 try:
1076 return float(brentq(lambda s, tgt=target: self._dk_dse(s) - tgt, _C_MIN, 1.0, xtol=BRENTQ_XTOL)) # type: ignore[arg-type]
1077 except ValueError:
1078 return _C_MIN
1080 def c_and_total_from_retardation(self, r: float) -> tuple[float, float]:
1081 """Return ``(c, C_T)`` at retardation ``r`` from a SINGLE brentq call.
1083 Overrides the default base-class implementation (which calls
1084 ``concentration_from_retardation`` and ``total_concentration``
1085 separately and ends up doing two independent brentq solves on the same
1086 underlying equation). Halves the iterative-solver cost in the IBP fan
1087 integrators.
1088 """
1089 s = self._se_from_retardation(r)
1090 c = max(self._k_se(s), _C_MIN)
1091 ct = self.delta_theta * s
1092 return c, ct
1095SorptionModel = NonlinearSorption | ConstantRetardation
1096"""Type alias for all sorption models accepted by the front-tracking solver."""
1099def characteristic_speed(c: float, sorption: SorptionModel) -> float:
1100 """Compute characteristic speed dV/dθ = 1/R(C).
1102 In (V, θ) coordinates, every characteristic propagates at a flow-free
1103 speed determined solely by the local concentration and the sorption
1104 isotherm.
1106 Parameters
1107 ----------
1108 c : float
1109 Dissolved concentration [mass/volume].
1110 sorption : SorptionModel
1111 Sorption model.
1113 Returns
1114 -------
1115 speed : float
1116 Characteristic speed dV/dθ.
1118 Examples
1119 --------
1120 >>> sorption = FreundlichSorption(
1121 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
1122 ... )
1123 >>> s = characteristic_speed(c=5.0, sorption=sorption)
1124 >>> s > 0
1125 True
1126 """
1127 r = float(sorption.retardation(c))
1128 return float("inf") if r == 0.0 else 1.0 / r
1131def characteristic_position(
1132 c: float,
1133 sorption: SorptionModel,
1134 theta_start: float,
1135 v_start: float,
1136 theta: float,
1137) -> float | None:
1138 """Compute position of a characteristic at cumulative flow θ.
1140 Characteristics propagate linearly in θ::
1142 V(θ) = v_start + characteristic_speed(C) * (θ - θ_start)
1144 Parameters
1145 ----------
1146 c : float
1147 Concentration carried by characteristic [mass/volume].
1148 sorption : SorptionModel
1149 Sorption model.
1150 theta_start : float
1151 Cumulative flow at which the characteristic starts [m³].
1152 v_start : float
1153 Starting position [m³].
1154 theta : float
1155 Cumulative flow at which to evaluate position [m³].
1157 Returns
1158 -------
1159 position : float or None
1160 Position at θ [m³], or None if θ < θ_start.
1162 Examples
1163 --------
1164 >>> sorption = ConstantRetardation(retardation_factor=2.0)
1165 >>> v = characteristic_position(
1166 ... c=5.0, sorption=sorption, theta_start=0.0, v_start=0.0, theta=1000.0
1167 ... )
1168 >>> bool(np.isclose(v, 500.0)) # v = (1/2) * 1000 = 500
1169 True
1170 """
1171 if theta < theta_start:
1172 return None
1174 return v_start + characteristic_speed(c, sorption) * (theta - theta_start)
1177def compute_first_front_arrival_theta(
1178 cin: npt.NDArray[np.floating],
1179 theta_edges: npt.NDArray[np.floating],
1180 aquifer_pore_volume: float,
1181 sorption: SorptionModel,
1182) -> float:
1183 """Cumulative-flow θ at which ``c_first`` arrives at the outlet (end of spin-up).
1185 "Arrival" means the θ at which the ``c_first`` *level* is fully present at
1186 the outlet, ``θ_emit + V·R(c_first)`` for ``n<1`` and
1187 ``θ_emit + V·C_T(c_first)/c_first`` for ``n>1``/constant retardation.
1189 .. warning::
1191 For ``n<1`` with ``c_min > 0`` (default ``c_min = 1e-12`` in
1192 :class:`FreundlichSorption`), the actual wave emitted is a
1193 :class:`~gwtransport.fronttracking.waves.RarefactionWave` whose head (``c = c_min ≈ 0``) reaches the
1194 outlet at θ ≈ ``V·R(c_min) ≈ V`` — *much* earlier than the value this
1195 function returns (which is the *tail* arrival ``V·R(c_first)``).
1196 The function returns "tail arrival" semantics: the returned θ is a
1197 conservative end-of-spin-up where c ≤ c_first everywhere before it.
1198 Consult the solver event log for the true rarefaction head crossing.
1200 Parameters
1201 ----------
1202 cin : numpy.ndarray
1203 Inlet concentration [mass/volume].
1204 theta_edges : numpy.ndarray
1205 Cumulative-flow edges; length ``len(cin) + 1``.
1206 aquifer_pore_volume : float
1207 Total pore volume [m³]. Must be positive.
1208 sorption : SorptionModel
1209 Sorption model.
1211 Returns
1212 -------
1213 theta_first_arrival : float
1214 Cumulative-flow θ at which ``c_first`` is fully present at the outlet
1215 [m³]. Returns ``np.inf`` only if ``cin`` is identically zero.
1217 Examples
1218 --------
1219 >>> cin = np.array([0.0, 10.0] + [10.0] * 10)
1220 >>> theta_edges = np.arange(0.0, 1300.0, 100.0) # constant flow=100, dt=1
1221 >>> sorption = ConstantRetardation(retardation_factor=2.0)
1222 >>> theta_first = compute_first_front_arrival_theta(
1223 ... cin, theta_edges, 500.0, sorption
1224 ... )
1225 >>> bool(np.isclose(theta_first, 100.0 + 500.0 * 2.0)) # θ_emit + V·R
1226 True
1227 """
1228 nonzero_indices = np.where(cin > 0)[0]
1229 if len(nonzero_indices) == 0:
1230 return float(np.inf)
1232 idx_first = int(nonzero_indices[0])
1233 c_first = float(cin[idx_first])
1235 if isinstance(sorption, FreundlichSorption) and sorption.n < 1.0:
1236 # n<1: the 0→c_first step emits a rarefaction; its tail (c=c_first)
1237 # reaches the outlet after V·R(c_first) units of cumulative flow.
1238 target_volume = aquifer_pore_volume * float(sorption.retardation(c_first))
1239 else:
1240 # n>1 or constant: R-H shock with speed = c / (C_T(c) - C_T(0));
1241 # target volume = V · C_T(c_first) / c_first.
1242 target_volume = aquifer_pore_volume * float(sorption.total_concentration(c_first)) / c_first
1244 return float(theta_edges[idx_first]) + target_volume