Coverage for src/gwtransport/percolation.py: 100%
72 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1"""
2Percolation through thick unsaturated zones via the Kinematic Wave method.
4This module provides one public function:
6- :func:`root_zone_to_water_table_kinematic_wave` — exact front-tracking
7 solver for gravity-driven percolation between the bottom of the root
8 zone and the water table, following the Kinematic-Wave method described
9 in Olsthoorn (2026, *Stromingen* 32(1)). Supports Brooks-Corey and
10 van Genuchten-Mualem constitutive curves and a time-varying
11 multiplicative scaling of K(θ) (e.g. for temperature-corrected
12 viscosity).
14**Forward-only.** Inverse mapping ``water_table_to_root_zone`` is not
15provided. The KW unsaturated-zone problem is fundamentally one-way under
16gravity: multiple ``q_root_zone(t)`` series produce indistinguishable
17``q_water_table(t)`` after the column's intrinsic low-pass response,
18making the inverse ill-posed. Users wanting an inverse should formulate
19it as a regularised inverse problem outside this package.
21**Cumulative pore-volume coordinate.** The position axis is *cumulative
22pore volume per unit cross-sectional area* (units of length), not
23geometric depth. For a soil of constant porosity ``n_p ≡ θ_s`` and
24water-table depth ``z_wt``, the conversion is ``V_out = θ_s · z_wt``.
25The docstring of :func:`root_zone_to_water_table_kinematic_wave`
26spells out the recovery rule and the layered-porosity generalisation.
28The full Kinematic-Wave derivation and the constitutive-curve references
29are documented on :func:`root_zone_to_water_table_kinematic_wave`.
31This file is part of gwtransport which is released under AGPL-3.0 license.
32See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
33"""
35import numpy as np
36import numpy.typing as npt
37import pandas as pd
39from gwtransport._time import tedges_to_days
40from gwtransport.advection import _flow_weighted_front_tracking_output
41from gwtransport.fronttracking.math import (
42 BrooksCoreyConductivity,
43 VanGenuchtenMualemConductivity,
44)
45from gwtransport.fronttracking.solver import FrontTracker
46from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave
49def root_zone_to_water_table_kinematic_wave(
50 *,
51 q_root_zone: npt.ArrayLike,
52 tedges: pd.DatetimeIndex,
53 q_water_table_tedges: pd.DatetimeIndex,
54 cumulative_pore_volumes_outlet: npt.ArrayLike,
55 theta_r: float,
56 theta_s: float,
57 k_s: float,
58 brooks_corey_lambda: float | None = None,
59 van_genuchten_n: float | None = None,
60 mualem_l: float = 0.5,
61 k_scaling: npt.ArrayLike | None = None,
62 max_iterations: int = 10000,
63) -> tuple[npt.NDArray[np.floating], list[dict]]:
64 r"""Percolation flux at the water table by exact Kinematic-Wave front tracking.
66 Solves the nonlinear scalar conservation law
68 .. math::
69 \\frac{\\partial \\theta_m}{\\partial t} +
70 \\frac{\\partial K(\\theta_m)}{\\partial z} = 0
72 exactly via :class:`gwtransport.fronttracking.solver.FrontTracker`,
73 using either a Brooks-Corey or a van Genuchten-Mualem constitutive
74 curve. Implements the Kinematic-Wave method (see [3]_ for the general
75 theory) described in Olsthoorn (2026) [1]_. The capillary term
76 ``∂ψ/∂z`` is dropped (gravity drainage only); real fronts are slightly
77 smoothed by capillarity, so if smoothing matters use the Munsflow-style
78 approach in :mod:`gwtransport.diffusion` instead.
80 Parameters
81 ----------
82 q_root_zone : array-like
83 Root-zone leakage entering the unsaturated zone at the top
84 boundary [length/time, e.g. m/day]. Piecewise constant over each
85 ``[tedges[i], tedges[i+1])`` bin. Non-negative.
86 Length = ``len(tedges) - 1``. At any bin, ``q_root_zone <= f·K_s``
87 must hold (with ``f = k_scaling`` or 1) for the inlet inversion
88 to be well-defined; the validator raises ``ValueError`` otherwise.
89 tedges : pandas.DatetimeIndex
90 Time bin edges of the input series. Length ``n + 1`` for ``n`` bins.
91 q_water_table_tedges : pandas.DatetimeIndex
92 Output time bin edges. Free monotone index when ``k_scaling`` is
93 None; **must equal** ``tedges`` when ``k_scaling`` is set (the
94 back-transform ``q_wt = f · cout`` is exact only on the input grid).
95 Must lie within the input window ``[tedges[0], tedges[-1]]`` (the flow
96 series defines the system only there); querying beyond it raises.
97 cumulative_pore_volumes_outlet : array-like
98 Cumulative pore volume per unit cross-sectional area at the water
99 table [length]. For a soil of constant porosity (``n_p ≡ θ_s``)
100 and water-table depth ``z_wt``, this is ``θ_s · z_wt``. For
101 layered porosity, ``∫₀^{z_wt} n_p(z') dz'``. The geometric depth
102 is recovered as ``z_wt = V_out / θ_s`` (uniform case). Array-like
103 to support a distribution of column lengths in parallel
104 (analogous to :func:`gwtransport.advection.gamma_infiltration_to_extraction`); each
105 entry must be positive.
106 theta_r : float
107 Residual volumetric moisture content [-]. Must satisfy
108 ``0 <= theta_r < theta_s``.
109 theta_s : float
110 Saturated volumetric moisture content [-]. Equal to the porosity
111 for typical soils. Must satisfy ``theta_r < theta_s < 1``.
112 k_s : float
113 Saturated hydraulic conductivity [length/time]. Positive.
114 brooks_corey_lambda : float or None, optional
115 Brooks-Corey pore-size distribution index [-]. Set to use the
116 Brooks-Corey branch. Mutually exclusive with ``van_genuchten_n``.
117 Tabulated soil values are available in the Staringreeks [2]_.
118 van_genuchten_n : float or None, optional
119 Van Genuchten shape parameter ``n_vG > 1``. Set to use the
120 van Genuchten-Mualem branch (numerical inversion via brentq).
121 Mutually exclusive with ``brooks_corey_lambda``.
122 mualem_l : float, optional
123 Mualem pore-connectivity parameter ``L``. Default 0.5
124 (standard Mualem). Honored only when ``van_genuchten_n`` is set.
125 k_scaling : array-like or None, optional
126 Dimensionless time-only multiplicative factor ``f(t)`` applied
127 to the entire ``K(θ)`` curve:
128 ``K(θ, t) = f(t) · K_reference(θ)``. Length ``n``. Default None
129 means ``f ≡ 1``. All entries must be strictly positive.
131 The cumulative-flow trick in the underlying front-tracking solver
132 absorbs ``f(t)`` exactly: wave dynamics in cumulative effective
133 time remain flow-free. Typical usage is a temperature-corrected
134 viscosity ``f(t) = μ_ref / μ(T(t))``; ``μ`` varies ~60% between
135 5 °C and 25 °C, so seasonal swings of 30-50% in effective ``K_s``
136 are realistic for shallow soils.
137 max_iterations : int, optional
138 Maximum number of solver events. Default 10000.
140 Returns
141 -------
142 q_water_table : ndarray
143 Bin-averaged percolation flux at the water table [same units as
144 ``q_root_zone``], length ``len(q_water_table_tedges) - 1``,
145 averaged across the columns in ``cumulative_pore_volumes_outlet``.
146 structures : list of dict
147 Per-column simulation structures (same schema as
148 :func:`gwtransport.advection.infiltration_to_extraction_nonlinear_sorption`,
149 with ``aquifer_pore_volume`` renamed to
150 ``cumulative_pore_volume_outlet``):
152 - ``waves`` — all wave objects.
153 - ``events`` — event history; each record has ``"theta"`` (cumulative
154 effective time) and ``"type"`` keys. Translate ``theta`` to wall-clock
155 time via ``tracker_state.t_at_theta(event["theta"])``.
156 - ``theta_first_arrival`` — cumulative effective time at which
157 the first nonzero arrival reaches the outlet.
158 - ``n_events``, ``n_shocks``, ``n_rarefactions``,
159 ``n_characteristics`` — counts.
160 - ``theta_current`` — final cumulative effective time.
161 - ``sorption`` — the sorption object.
162 - ``tracker_state`` — complete :class:`~gwtransport.fronttracking.solver.FrontTrackerState` for the
163 column (use ``state.t_at_theta`` to translate ``θ → t``).
164 - ``cumulative_pore_volume_outlet`` — the V_out for this column.
166 Raises
167 ------
168 ValueError
169 If inputs are inconsistent (wrong lengths, NaN, negative ``q_root_zone``
170 or ``k_scaling``, non-finite or non-positive
171 ``cumulative_pore_volumes_outlet`` or ``k_s``), if neither or both
172 sorption-parameter groups are supplied,
173 if ``q_root_zone > f(t) * k_s`` at any bin (saturation/ponding limit),
174 or if ``q_water_table_tedges`` does not equal ``tedges`` while
175 ``k_scaling`` is provided.
177 Warns
178 -----
179 UserWarning
180 If output θ-bins extend beyond the inlet θ-window (i.e. the drying tail
181 of ``q_root_zone`` reaches zero and the column has not yet equilibrated
182 by the last output bin). Bin averages in that region are clamped to zero.
184 See Also
185 --------
186 gwtransport.advection.infiltration_to_extraction_nonlinear_sorption :
187 Solute transport with nonlinear sorption (analogous front-tracking
188 algorithm in the saturated-zone domain).
189 gwtransport.diffusion :
190 Munsflow-style linearised advection-diffusion (complementary;
191 smoothed fronts).
192 gwtransport.fronttracking.math.BrooksCoreyConductivity :
193 Brooks-Corey constitutive class.
194 gwtransport.fronttracking.math.VanGenuchtenMualemConductivity :
195 van Genuchten-Mualem constitutive class.
196 :ref:`concept-kinematic-wave` : Background on the Kinematic-Wave method for
197 unsaturated-zone percolation.
199 Notes
200 -----
201 **Cumulative pore-volume coordinate.** The internal V axis is
202 ``V(z) = int_0^z n_p(z') dz'`` (units of length). For a uniform soil
203 with ``n_p = theta_s``, ``V = theta_s * z``; depth is recovered as
204 ``z = V / theta_s``. The solver-side identification
205 ``flow = theta_s * f(t)`` (with ``f`` the optional K-scaling) follows
206 from the chain rule ``d/dz = theta_s * d/dV``.
208 **Inlet boundary inversion.** The solver works in a reference frame
209 where ``K = K_ref(theta_m)``; the time-varying scaling is moved to the
210 boundary as ``cin_solver(t) = q_root_zone(t) / f(t)`` and recovered
211 at the outlet as ``q_water_table(t) = f(t) * cout(t)``. The
212 requirement ``cin_solver <= k_s`` (i.e. ``q_root_zone <= f * k_s``)
213 is the saturation/ponding admissibility check enforced by the
214 validator.
216 **The KW approximation.** Capillary stresses are neglected; flow
217 is gravity-only. Wetting fronts are sharp shocks satisfying
218 Rankine-Hugoniot ``V_f = (K_1 - K_2)/(theta_1 - theta_2)``. Drying tails are
219 self-similar rarefaction fans. Real fronts are slightly capillary-
220 smoothed; if that smoothing matters, use Munsflow-style
221 advection-diffusion (the article's Munsflow method, mapped to
222 :mod:`gwtransport.diffusion` in this package).
224 **Initial condition.** The column starts at ``theta_m = theta_r`` (i.e.
225 ``K = 0``) everywhere. To start from field capacity or a long-term
226 equilibrium, prepend a constant-q spin-up to the input series.
228 **Exact mass conservation.** Both Brooks-Corey and van Genuchten-Mualem
229 fan integrals use a closed-form integration-by-parts antiderivative
230 derived from the universal identity ``R = dC_T/dC``: for the spatial
231 fan integral ``G(u) = C_T(c) * u - kappa * c``, and for the temporal
232 fan integral ``F(theta) = c * (theta - theta_origin) - Delta_v * C_T(c)``.
233 For Brooks-Corey both ``c`` and ``C_T`` at the endpoints are closed
234 form; for van Genuchten-Mualem they require a single ``brentq`` call
235 per endpoint (transcendental ``K(theta)``). The Burdine variant
236 (``mualem_l = 0``) admits a closed-form inverse and is fully
237 free of root-finding.
239 References
240 ----------
241 .. [1] Olsthoorn, T.N. (2026). Percolation through thick unsaturated
242 zones — Munsflow vs. the Kinematic Wave. *Stromingen* 32(1).
243 .. [2] Heinen, M., Bakker, G., Wösten, J.M.H. (2020). *Waterretentie
244 en Doorlatendheidskarakteristieken van boven- en ondergronden in
245 Nederland: de Staringreeks. Update 2018.* Wageningen Environmental
246 Research, Report 2978.
247 .. [3] Charbeneau, R.J. (2000). *Groundwater Hydraulics and Pollutant
248 Transport.* Prentice Hall.
250 Examples
251 --------
252 .. disable_try_examples
254 Reproduce a 10-year step-response for the article's soil O05
255 (coarse sand, Brooks-Corey)::
257 import numpy as np
258 import pandas as pd
259 from gwtransport.percolation import (
260 root_zone_to_water_table_kinematic_wave,
261 )
263 tedges = pd.date_range("1995-01-01", "2005-01-01", freq="D")
264 q_root = np.full(len(tedges) - 1, 1e-3) # 1 mm/day
266 q_wt, structures = root_zone_to_water_table_kinematic_wave(
267 q_root_zone=q_root,
268 tedges=tedges,
269 q_water_table_tedges=tedges,
270 cumulative_pore_volumes_outlet=np.array([0.337 * 20.0]),
271 theta_r=0.01,
272 theta_s=0.337,
273 k_s=0.174,
274 brooks_corey_lambda=0.25,
275 )
277 With time-varying water viscosity::
279 days = ((tedges[:-1] - tedges[0]) / pd.Timedelta(days=1)).values
280 T = 10.0 + 5.0 * np.sin(2 * np.pi * days / 365.25) # °C
281 mu_ref, dmu_dT = 1.31, -0.027 # mPa·s, linear around 10 °C
282 mu = mu_ref + dmu_dT * (T - 10.0)
283 k_scaling = mu_ref / mu
285 q_wt_visc, _ = root_zone_to_water_table_kinematic_wave(
286 q_root_zone=q_root,
287 tedges=tedges,
288 q_water_table_tedges=tedges,
289 cumulative_pore_volumes_outlet=np.array([0.337 * 20.0]),
290 theta_r=0.01,
291 theta_s=0.337,
292 k_s=0.174,
293 brooks_corey_lambda=0.25,
294 k_scaling=k_scaling,
295 )
296 """
297 q_root_zone_arr = np.asarray(q_root_zone, dtype=float)
298 # Promote a scalar / 0-d outlet pore volume to a single-column 1-d array so the
299 # downstream per-column iteration and len() are well-defined (a lone scalar is a
300 # valid single column, consistent with the "array-like" contract in the docstring).
301 cumulative_pore_volumes_outlet_arr = np.atleast_1d(np.asarray(cumulative_pore_volumes_outlet, dtype=float))
302 tedges = pd.DatetimeIndex(tedges)
303 q_water_table_tedges = pd.DatetimeIndex(q_water_table_tedges)
305 n_bins = len(q_root_zone_arr)
306 if len(tedges) != n_bins + 1:
307 msg = f"tedges must have length len(q_root_zone) + 1, got {len(tedges)} vs {n_bins + 1}"
308 raise ValueError(msg)
309 if q_water_table_tedges[0] < tedges[0] or q_water_table_tedges[-1] > tedges[-1]:
310 msg = (
311 f"q_water_table_tedges must lie within the input window [{tedges[0]}, {tedges[-1]}], got "
312 f"[{q_water_table_tedges[0]}, {q_water_table_tedges[-1]}]. The flow series defines the "
313 "system only over the input window; querying beyond it is ill-posed (extend q_root_zone instead)."
314 )
315 raise ValueError(msg)
316 if np.any(q_root_zone_arr < 0):
317 msg = "q_root_zone must be non-negative"
318 raise ValueError(msg)
319 if np.any(np.isnan(q_root_zone_arr)):
320 msg = "q_root_zone must not contain NaN"
321 raise ValueError(msg)
322 if cumulative_pore_volumes_outlet_arr.size == 0 or not np.all(
323 np.isfinite(cumulative_pore_volumes_outlet_arr) & (cumulative_pore_volumes_outlet_arr > 0)
324 ):
325 msg = "cumulative_pore_volumes_outlet must be non-empty with all entries positive and finite"
326 raise ValueError(msg)
327 if not (0.0 <= theta_r < theta_s < 1.0):
328 msg = f"theta_r, theta_s must satisfy 0 <= theta_r < theta_s < 1, got theta_r={theta_r}, theta_s={theta_s}"
329 raise ValueError(msg)
330 if not (np.isfinite(k_s) and k_s > 0):
331 msg = f"k_s must be positive and finite, got {k_s}"
332 raise ValueError(msg)
333 if (brooks_corey_lambda is None) == (van_genuchten_n is None):
334 msg = "Exactly one of brooks_corey_lambda or van_genuchten_n must be provided"
335 raise ValueError(msg)
337 if k_scaling is None:
338 f = np.ones(n_bins, dtype=float)
339 else:
340 f = np.asarray(k_scaling, dtype=float)
341 if f.shape != (n_bins,):
342 msg = f"k_scaling must have shape ({n_bins},), got {f.shape}"
343 raise ValueError(msg)
344 if np.any(np.isnan(f)) or np.any(f <= 0):
345 msg = "k_scaling must be strictly positive and contain no NaN"
346 raise ValueError(msg)
347 if len(q_water_table_tedges) != len(tedges) or not (q_water_table_tedges == tedges).all():
348 msg = (
349 "q_water_table_tedges must equal tedges when k_scaling is provided "
350 "(the back-transform q_wt = f * cout is exact only on the input grid)"
351 )
352 raise ValueError(msg)
354 # Saturation/ponding admissibility: K_ref(θ_m at inlet) = q_root/f must be <= k_s.
355 cin_solver = q_root_zone_arr / f
356 if float(cin_solver.max()) > k_s:
357 bin_idx = int(cin_solver.argmax())
358 msg = (
359 f"Inlet saturation/ponding limit exceeded at bin {bin_idx}: "
360 f"q_root_zone/k_scaling = {cin_solver[bin_idx]:.6g} > k_s = {k_s:.6g}. "
361 "Reduce q_root_zone or increase k_scaling (warmer water → lower viscosity → higher k_s effective)."
362 )
363 raise ValueError(msg)
365 if brooks_corey_lambda is not None:
366 sorption: BrooksCoreyConductivity | VanGenuchtenMualemConductivity = BrooksCoreyConductivity(
367 theta_r=theta_r, theta_s=theta_s, k_s=k_s, brooks_corey_lambda=brooks_corey_lambda
368 )
369 else:
370 assert van_genuchten_n is not None # noqa: S101 # narrowed by validation above
371 sorption = VanGenuchtenMualemConductivity(
372 theta_r=theta_r, theta_s=theta_s, k_s=k_s, van_genuchten_n=van_genuchten_n, mualem_l=mualem_l
373 )
375 # Solver-frame arrays: flow_solver = θ_s · f(t) (porosity ≡ θ_s for unsaturated KW);
376 # cin_solver = q_root/f was already formed for the admissibility check above.
377 flow_solver = theta_s * f
379 flow_tedges_days = tedges_to_days(tedges)
380 cout_tedges_days = tedges_to_days(q_water_table_tedges, ref=tedges[0])
381 n_out = len(q_water_table_tedges) - 1
383 q_wt_all = np.zeros((len(cumulative_pore_volumes_outlet_arr), n_out))
384 structures: list[dict] = []
386 for i, v_out in enumerate(cumulative_pore_volumes_outlet_arr):
387 tracker = FrontTracker(
388 cin=cin_solver,
389 flow=flow_solver,
390 tedges=tedges,
391 aquifer_pore_volume=float(v_out),
392 sorption=sorption,
393 )
394 tracker.run(max_iterations=max_iterations)
396 cout_ref = _flow_weighted_front_tracking_output(
397 cout_tedges_days=cout_tedges_days,
398 flow_tedges_days=flow_tedges_days,
399 flow=flow_solver,
400 v_outlet=float(v_out),
401 waves=tracker.state.waves,
402 sorption=sorption,
403 theta_edges=tracker.state.theta_edges,
404 cin=cin_solver,
405 )
407 # Back-transform to physical flux: q_wt = f · cout_ref. When k_scaling is set the
408 # validator requires q_water_table_tedges == tedges, so f aligns with cout_ref. When
409 # k_scaling is None, f ≡ 1 and the output grid may be coarser than the input, so
410 # cout_ref is used directly (multiplying by the input-length f would mis-broadcast).
411 q_wt_all[i, :] = cout_ref if k_scaling is None else f * cout_ref
413 structures.append({
414 "waves": tracker.state.waves,
415 "events": tracker.state.events,
416 "theta_first_arrival": tracker.theta_first_arrival,
417 "n_events": len(tracker.state.events),
418 "n_shocks": sum(1 for w in tracker.state.waves if isinstance(w, ShockWave)),
419 "n_rarefactions": sum(1 for w in tracker.state.waves if isinstance(w, RarefactionWave)),
420 "n_characteristics": sum(1 for w in tracker.state.waves if isinstance(w, CharacteristicWave)),
421 "theta_current": tracker.state.theta_current,
422 "sorption": sorption,
423 "tracker_state": tracker.state,
424 "cumulative_pore_volume_outlet": float(v_out),
425 })
427 return np.mean(q_wt_all, axis=0), structures