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