Coverage for src/gwtransport/deposition.py: 0%
118 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 20:54 +0000
1"""
2Deposition Analysis for 1D Aquifer Systems.
4Areal deposition supplies mass to the groundwater, mixed instantaneously over the height of the
5aquifer. The aquifer has a constant thickness with a finite pore volume; water with zero
6concentration infiltrates at one end and is extracted at the other, whether the flow is radial or
7orthogonal. Transport is 1D advection with linear sorption; there is no microdispersion, molecular
8diffusion, or macrodispersion. Forward and backward modeling are supported.
10The model is a *source* term (positive deposition adds mass to the water); it does NOT model removal
11processes such as pathogen attachment, particle filtration, or chemical precipitation, which would
12remove mass from the water and require the opposite sign convention.
14Available functions:
16- :func:`deposition_to_extraction` - Forward model: from areal deposition rates [g/m²/day] and
17 flow, compute the concentration [g/m³] of the extracted water on ``cout_tedges``. Each output
18 bin is the flow-weighted average over the water extracted in that bin; a parcel's concentration
19 gain is proportional to the residence-time contact with the areal flux, mixed over
20 ``porosity * thickness``. Bins whose residence time is not yet resolved (spin-up) return NaN,
21 while bins that extract zero volume -- including bins lying outside the flow record -- return 0.0.
23- :func:`extraction_to_deposition` - Inverse model (deconvolution): from the concentration of the
24 extracted water [g/m³], estimate the mean deposition rate [g/m²/day] in each ``tedges`` bin.
25 Solves the banded forward system with Tikhonov regularization toward a physically motivated
26 target (transpose-and-normalize of the forward operator), with ``regularization_strength``
27 trading data fit against that target. Rows without a resolved residence time, zero-flow output
28 bins and NaN entries in ``cout`` are excluded from the solve. This is the recommended inverse.
30- :func:`extraction_to_deposition_full` - Same inversion routed through the dense nullspace solver
31 :func:`gwtransport.utils.solve_underdetermined_system`, exposing its options: the nullspace
32 objective (``"squared_differences"`` for smooth solutions, ``"summed_differences"`` for
33 piecewise-constant ones, or a callable), the scipy ``optimization_method``, and the
34 least-squares ``rcond`` cutoff. Reach for it when the choice of nullspace component matters;
35 otherwise prefer :func:`extraction_to_deposition`, which is cheaper and needs no dense operator.
37- :func:`compute_deposition_weights` - Build the operator relating deposition rates to extracted
38 concentrations in a compact banded layout, returning the band values, the first cin column of
39 each row, and the row masks for valid and spin-up output bins. Row ``k`` sums to
40 ``residence_time_k / (retardation_factor * porosity * thickness)`` rather than to one, because a
41 deposition rate maps to a concentration through the residence time. Exposed for custom inverse
42 solvers; the forward and both inverse entry points build on it.
44- :func:`spinup_duration` - Time in days, measured from ``tedges[0]``, at which the cumulative flow
45 first reaches ``retardation_factor * aquifer_pore_volume``: the earliest extraction time whose
46 water carries a complete deposition history, and hence the start of the valid analysis period.
47 Under constant flow this is ``retardation_factor * aquifer_pore_volume / flow``. Raises
48 ``ValueError`` when the flow record is too short to reach that volume.
50This file is part of gwtransport which is released under AGPL-3.0 license.
51See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
52"""
54from collections.abc import Callable
56import numpy as np
57import numpy.typing as npt
58import pandas as pd
60from gwtransport._time import tedges_to_days
61from gwtransport._validation import (
62 _validate_no_nan,
63 _validate_non_negative_array,
64 _validate_positive_scalar,
65 _validate_retardation_factor,
66 _validate_tedges_parity,
67)
68from gwtransport.advection_utils import _densify_weights, _resolve_spinup_inputs
69from gwtransport.deposition_utils import _clipped_linear_integral
70from gwtransport.utils import (
71 _make_strictly_monotone,
72 cumulative_flow_volume,
73 linear_interpolate,
74 solve_inverse_transport_banded,
75 solve_underdetermined_system,
76)
79def _validate_deposition_inputs(
80 *,
81 tedges: pd.DatetimeIndex,
82 flow_values: np.ndarray,
83 aquifer_pore_volume: float,
84 porosity: float,
85 thickness: float,
86 retardation_factor: float = 1.0,
87 spinup: str | None = "constant",
88 cout_tedges: pd.DatetimeIndex | None = None,
89 cout_values: np.ndarray | None = None,
90 dep_values: np.ndarray | None = None,
91) -> None:
92 """Validate inputs common to deposition forward / reverse / full entry points.
94 Activates checks per the kwargs that are not None:
96 - ``dep_values`` provided => ``tedges``-parity vs ``dep`` + a combined dep+flow
97 NaN-check (forward path) under the "Input arrays cannot contain NaN values"
98 message.
99 - ``cout_values`` + ``cout_tedges`` provided => ``cout_tedges``-parity check
100 (inverse paths). ``cout_values`` itself is intentionally NOT NaN-checked
101 -- NaN in ``cout`` is allowed and excluded downstream by the inverse solve.
102 - ``flow_values`` + ``tedges`` always => parity check + non-negative;
103 additionally, in the inverse path (``dep_values is None``), a flow-only
104 NaN-check fires with the "flow array cannot contain NaN values" message.
105 - Physical params (``porosity``, ``thickness``, ``aquifer_pore_volume``,
106 ``retardation_factor``) always validated.
108 Raises
109 ------
110 ValueError
111 If any of the activated checks fails. The specific message names which
112 invariant was violated (see body for the verbatim strings).
113 NotImplementedError
114 If ``spinup`` is anything other than ``None`` or ``"constant"``. The
115 fraction-threshold mode is not implemented for deposition (matching the
116 diffusion family); floats, ints, and bools are all rejected.
117 """
118 if spinup is not None and spinup != "constant":
119 # Accept only None and "constant"; reject everything else (floats, ints, bools,
120 # typo'd strings). The deposition callers discard the fraction-threshold value from
121 # _resolve_spinup_inputs, so any other request would be silently ignored.
122 msg = (
123 "deposition's spinup parameter only supports None or 'constant'; "
124 f"other values are not yet implemented (got {spinup!r})"
125 )
126 raise NotImplementedError(msg)
127 if dep_values is not None:
128 _validate_tedges_parity(tedges, dep_values, tedges_name="tedges", values_name="dep")
129 _validate_tedges_parity(tedges, flow_values, tedges_name="tedges", values_name="flow")
130 if cout_values is not None and cout_tedges is not None:
131 _validate_tedges_parity(cout_tedges, cout_values, tedges_name="cout_tedges", values_name="cout")
132 if dep_values is not None:
133 # Compound NaN-check covers both ``dep`` and ``flow`` under one message; mapping to
134 # two separate ``_validate_no_nan`` calls would change wording and which array name
135 # surfaces first.
136 if np.any(np.isnan(dep_values)) or np.any(np.isnan(flow_values)):
137 msg = "Input arrays cannot contain NaN values"
138 raise ValueError(msg)
139 else:
140 _validate_no_nan(flow_values, name="flow", message="flow array cannot contain NaN values")
141 _validate_non_negative_array(
142 flow_values, name="flow", message="flow must be non-negative (negative flow not supported)"
143 )
144 if not 0 < porosity < 1:
145 msg = f"Porosity must be in (0, 1), got {porosity}"
146 raise ValueError(msg)
147 _validate_positive_scalar(thickness, name="thickness", message=f"Thickness must be positive, got {thickness}")
148 _validate_positive_scalar(
149 aquifer_pore_volume,
150 name="aquifer_pore_volume",
151 message=f"Aquifer pore volume must be positive, got {aquifer_pore_volume}",
152 )
153 _validate_retardation_factor(retardation_factor)
156def compute_deposition_weights(
157 *,
158 flow: npt.ArrayLike,
159 tedges: pd.DatetimeIndex,
160 cout_tedges: pd.DatetimeIndex,
161 aquifer_pore_volume: float,
162 porosity: float,
163 thickness: float,
164 retardation_factor: float = 1.0,
165) -> tuple[
166 npt.NDArray[np.floating],
167 npt.NDArray[np.intp],
168 npt.NDArray[np.bool_],
169 npt.NDArray[np.bool_],
170]:
171 """Build the deposition weight operator in a compact banded layout.
173 Row ``k`` of the dense ``(n_cout, n_cin)`` operator is ``band_vals[k]``
174 placed at columns ``[col_start[k], col_start[k] + full_band)``. The operator
175 is genuinely banded -- row ``k`` is nonzero only on the cin bins whose
176 cumulative through-flow volume lies in the residence-time window
177 ``[min(start_vol_k, start_vol_{k+1}), max(start_vol_k, start_vol_{k+1}) +
178 R * aquifer_pore_volume]`` -- so each band has at most ``full_band`` slots,
179 bounded by ``R * aquifer_pore_volume`` in volume (independent of record
180 length ``n_cin``). The window is located by :func:`numpy.searchsorted` on the
181 cumulative flow volume ``flow_cum``; the per-cell math reuses
182 ``gwtransport.deposition_utils._clipped_linear_integral`` restricted to
183 the band columns, so each row sums to
184 ``r_k = residence_time_k / (retardation_factor * porosity * thickness)``.
185 Reconstruct the dense
186 ``(n_cout, n_cin)`` matrix with
187 ``gwtransport.advection_utils._densify_weights`` when a dense operator
188 is required (the nullspace inverse).
190 Parameters
191 ----------
192 flow : array-like
193 Flow rates in aquifer [m³/day]. Length must equal ``len(tedges) - 1``.
194 tedges : pandas.DatetimeIndex
195 Time bin edges for flow data.
196 cout_tedges : pandas.DatetimeIndex
197 Time bin edges for output concentration data.
198 aquifer_pore_volume : float
199 Aquifer pore volume [m³].
200 porosity : float
201 Aquifer porosity [dimensionless].
202 thickness : float
203 Aquifer thickness [m].
204 retardation_factor : float, optional
205 Compound retardation factor, by default 1.0.
207 Returns
208 -------
209 band_vals : numpy.ndarray
210 Banded weights of shape ``(n_cout, full_band)``. Slot ``band_vals[k, b]``
211 is the weight on cin bin ``col_start[k] + b``. Row ``k`` sums to
212 ``r_k = residence_time_k / (retardation_factor * porosity * thickness)``;
213 invalid rows (NaN residence time, zero-flow cout bins) are zero.
214 col_start : numpy.ndarray of int
215 First cin bin index of each cout row's band, shape ``(n_cout,)``.
216 row_valid : numpy.ndarray of bool
217 True for cout bins whose residence-time window is fully defined and
218 carries flow (the finite, nonzero rows), shape ``(n_cout,)``.
219 spinup_row : numpy.ndarray of bool
220 True for cout bins whose residence time is undefined (spin-up period),
221 shape ``(n_cout,)``. These rows carry an all-zero band; the forward
222 path returns NaN for these bins (distinct from zero-flow cout bins,
223 which return 0).
225 See Also
226 --------
227 ``gwtransport.advection_utils._densify_weights`` : Reconstruct the dense matrix.
228 """
229 t0 = tedges[0]
230 tedges_days = tedges_to_days(tedges, ref=t0)
231 cout_tedges_days = tedges_to_days(cout_tedges, ref=t0)
233 flow_values = np.asarray(flow, dtype=float)
234 flow_cum = cumulative_flow_volume(flow_values, np.diff(tedges_days))
235 end_vol = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days)
236 r_apv = retardation_factor * float(aquifer_pore_volume)
238 # Infiltration-side cumulative volume of each cout edge is its extraction-side volume minus the
239 # retarded pore volume -- the direct cumulative-volume identity, avoiding the residence-time
240 # round-trip. NaN where the cout edge is outside the flow record or the look-back precedes the
241 # record start (the spin-up NaN the round-trip produced).
242 in_record = (cout_tedges_days >= tedges_days[0]) & (cout_tedges_days <= tedges_days[-1])
243 start_vol = end_vol - r_apv
244 start_vol = np.where(in_record & (start_vol >= flow_cum[0]), start_vol, np.nan)
246 n_cin = len(tedges) - 1
247 n_cout = len(cout_tedges) - 1
248 extracted_volume = np.diff(end_vol)
249 dt = np.diff(tedges_days)
251 # Row k's clipped trapezoid spans cout edges k (top: start_vol[k]) and k+1
252 # (bottom: start_vol[k+1]). It is nonzero only on cin bins j whose cumulative
253 # volume window [flow_cum[j], flow_cum[j+1]] overlaps the clip window
254 # [lower_k, upper_k] in flow_cum space, located by searchsorted.
255 sv_top, sv_bot = start_vol[:-1], start_vol[1:]
256 nan_row = np.isnan(sv_top) | np.isnan(sv_bot)
257 lower = np.minimum(sv_top, sv_bot)
258 upper = np.maximum(sv_top, sv_bot) + r_apv
259 j_first = np.clip(np.searchsorted(flow_cum, np.where(nan_row, flow_cum[0], lower), side="right") - 1, 0, n_cin - 1)
260 j_last = np.clip(np.searchsorted(flow_cum, np.where(nan_row, flow_cum[0], upper), side="left") - 1, 0, n_cin - 1)
261 width = np.where(nan_row, 0, j_last - j_first + 1)
262 full_band = int(max(1, width.max(initial=0)))
264 col_start = np.where(nan_row, 0, j_first).astype(np.intp)
265 row_valid = ~nan_row & (extracted_volume > 0)
266 # A row goes NaN in the forward only when its residence time is undefined AND
267 # it extracts water (left-edge spin-up). Out-of-range rows have undefined
268 # residence but zero extracted volume; they stay at the zeros sentinel (0).
269 spinup_row = nan_row & (extracted_volume > 0)
271 # Gather each row's band of cin edges (full_band + 1 edges) and bin widths,
272 # then evaluate the clipped-trapezoid integral on those columns only. Out-of-
273 # range / right-pad slots gather the last edge (zero-width contribution).
274 band_edge = col_start[:, None] + np.arange(full_band + 1)[None, :]
275 band_edge_clipped = np.clip(band_edge, 0, n_cin)
276 y_at_edge = flow_cum[band_edge_clipped] # (n_cout, full_band + 1)
277 y_top = y_at_edge - sv_top[:, None]
278 y_bot = y_at_edge - sv_bot[:, None]
279 widths = dt[np.clip(band_edge[:, :-1], 0, n_cin - 1)]
280 # Zero the width of right-pad slots (band slot index >= width) so they add nothing.
281 slot = np.arange(full_band)[None, :]
282 widths = np.where(slot < width[:, None], widths, 0.0)
284 # _clipped_linear_integral returns the volume*time measure (m³*day) of the
285 # parcel's residence-window overlap with each cin bin; dividing by
286 # (porosity*thickness) converts that overlap to plan-area * time, the areal
287 # footprint (and duration) over which the areal deposition flux is mixed down
288 # through the aquifer thickness. The bin width is already folded into the
289 # integral, so do NOT multiply by widths again. The R-scaled window (r_apv)
290 # stretches this measure by R, but the steady-state areal mass balance
291 # Q*cout = dep*A_plan is R-independent (retardation delays breakthrough, it
292 # does not raise the outlet concentration), so divide out that same R here.
293 top_integral = _clipped_linear_integral(y_top[:, :-1], y_top[:, 1:], widths, 0.0, r_apv)
294 bottom_integral = _clipped_linear_integral(y_bot[:, :-1], y_bot[:, 1:], widths, 0.0, r_apv)
295 contact_volume_time = np.maximum(top_integral - bottom_integral, 0.0)
296 numerator = contact_volume_time / (porosity * thickness * retardation_factor)
298 band_vals = np.zeros((n_cout, full_band))
299 # row_valid implies extracted_volume > 0, so the masked divide never sees a zero.
300 band_vals[row_valid] = numerator[row_valid] / extracted_volume[row_valid, None]
301 return band_vals, col_start, row_valid, spinup_row
304def deposition_to_extraction(
305 *,
306 dep: npt.ArrayLike,
307 flow: npt.ArrayLike,
308 tedges: pd.DatetimeIndex | np.ndarray,
309 cout_tedges: pd.DatetimeIndex | np.ndarray,
310 aquifer_pore_volume: float,
311 porosity: float,
312 thickness: float,
313 retardation_factor: float = 1.0,
314 spinup: str | None = "constant",
315) -> npt.NDArray[np.floating]:
316 """Compute concentrations from deposition rates (convolution).
318 Parameters
319 ----------
320 dep : array-like
321 Deposition rates [g/m²/day]. Length must equal len(tedges) - 1.
322 flow : array-like
323 Flow rates in aquifer [m³/day]. Length must equal len(tedges) - 1. The model
324 assumes this value is constant over each interval ``[tedges[i], tedges[i+1])``.
325 tedges : pandas.DatetimeIndex
326 Time bin edges for deposition and flow data.
327 cout_tedges : pandas.DatetimeIndex
328 Time bin edges for output concentration data.
329 aquifer_pore_volume : float
330 Aquifer pore volume [m³].
331 porosity : float
332 Aquifer porosity [dimensionless].
333 thickness : float
334 Aquifer thickness [m].
335 retardation_factor : float, optional
336 Compound retardation factor, by default 1.0.
337 spinup : {"constant"} | None, optional
338 Spin-up policy applied before computing deposition weights.
339 Default ``"constant"`` shifts ``tedges[0]`` backward by at least
340 ``retardation_factor * aquifer_pore_volume / flow[0]`` (rounded up to
341 whole bins, plus one extra bin) and treats ``dep`` and ``flow`` as
342 constant at their first observed values over the prepended interval.
343 When the warm-start is undefined -- ``flow[0] <= 0``, zero/NaN pore
344 volume, a non-positive first bin width, or a pad exceeding the memory
345 cap -- ``"constant"`` silently reverts to the strict-validity
346 behavior. ``None`` keeps the existing strict-validity behavior (NaN
347 cout rows during spin-up). A float raises ``NotImplementedError`` --
348 the fraction-threshold mode is not implemented for deposition
349 (matching the diffusion family).
351 Returns
352 -------
353 numpy.ndarray
354 Concentration changes [g/m³] with length len(cout_tedges) - 1.
356 Zero-extraction-flow cout bins (no water leaves the aquifer over the
357 bin) return ``0.0``, not NaN. This deliberately differs from advection,
358 which returns NaN for its undefined zero-flow output: the deposition
359 source term is defined even with no water (an areal flux still supplies
360 mass), and a bin that extracts zero volume carries zero mass, so ``0.0``
361 is the physically correct value rather than an undefined result. NaN is
362 reserved for spin-up bins whose residence time is not yet resolved.
364 Cout bins lying entirely outside the flow record (before ``tedges[0]``
365 or after ``tedges[-1]``) also return ``0.0``: their out-of-record edges
366 clamp to the record boundaries, so they extract zero volume and fall
367 under the same zero-mass convention.
369 Raises
370 ------
371 ValueError
372 If tedges does not have one more element than dep or flow, if input
373 arrays contain NaN values, if ``flow`` contains negative values, if
374 ``retardation_factor`` is less than 1.0, or if physical parameters are
375 out of valid range (porosity not in (0, 1), non-positive thickness or
376 aquifer pore volume).
377 NotImplementedError
378 If ``spinup`` is anything other than ``None`` or ``"constant"`` (the
379 fraction-threshold mode is not implemented for deposition).
381 See Also
382 --------
383 extraction_to_deposition : Inverse operation (deconvolution)
384 spinup_duration : Earliest extraction time with a fully resolved deposition history
385 gwtransport.advection.infiltration_to_extraction : For concentration transport without deposition
386 :ref:`concept-transport-equation` : Flow-weighted averaging approach
388 Notes
389 -----
390 This is a *source* term -- positive ``dep`` raises ``cout``. Sink
391 processes (pathogen attachment, first-order decay, particle filtration)
392 require the opposite sign convention and are not modelled here.
394 Examples
395 --------
396 >>> import pandas as pd
397 >>> import numpy as np
398 >>> from gwtransport.deposition import deposition_to_extraction
399 >>> dates = pd.date_range("2020-01-01", "2020-01-10", freq="D")
400 >>> tedges = pd.date_range("2019-12-31 12:00", "2020-01-10 12:00", freq="D")
401 >>> cout_tedges = pd.date_range("2020-01-03 12:00", "2020-01-12 12:00", freq="D")
402 >>> dep = np.ones(len(dates))
403 >>> flow = np.full(len(dates), 100.0)
404 >>> cout = deposition_to_extraction(
405 ... dep=dep,
406 ... flow=flow,
407 ... tedges=tedges,
408 ... cout_tedges=cout_tedges,
409 ... aquifer_pore_volume=500.0,
410 ... porosity=0.3,
411 ... thickness=10.0,
412 ... )
413 >>> print(f"First finite cout: {cout[np.isfinite(cout)][0]:.4f} g/m³")
414 First finite cout: 1.6667 g/m³
415 """
416 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
417 dep_values, flow_values = np.asarray(dep), np.asarray(flow)
419 _validate_deposition_inputs(
420 tedges=tedges,
421 flow_values=flow_values,
422 aquifer_pore_volume=aquifer_pore_volume,
423 porosity=porosity,
424 thickness=thickness,
425 retardation_factor=retardation_factor,
426 spinup=spinup,
427 dep_values=dep_values,
428 )
430 # Apply spinup policy: optionally prepend warm-start bins to tedges/flow/dep.
431 weight_tedges, weight_flow, weight_dep, _, _ = _resolve_spinup_inputs(
432 spinup,
433 tedges=tedges,
434 flow=flow_values,
435 aquifer_pore_volumes=np.array([aquifer_pore_volume]),
436 retardation_factor=retardation_factor,
437 cin=dep_values,
438 )
439 assert weight_dep is not None # noqa: S101 -- narrowed: cin was passed in
441 # Build the banded forward operator and apply it as a banded einsum instead of
442 # a dense W.dot(dep). Spin-up rows (NaN residence time) carry an all-zero band
443 # and must return NaN. Zero-flow cout bins (extracted_volume == 0) carry a zero
444 # band and return 0.
445 band_vals, col_start, _, spinup_row = compute_deposition_weights(
446 flow=weight_flow,
447 tedges=weight_tedges,
448 cout_tedges=cout_tedges,
449 aquifer_pore_volume=aquifer_pore_volume,
450 porosity=porosity,
451 thickness=thickness,
452 retardation_factor=retardation_factor,
453 )
454 n_cin = len(weight_tedges) - 1
455 cols = np.clip(col_start[:, None] + np.arange(band_vals.shape[1]), 0, n_cin - 1)
456 cout = np.einsum("kb,kb->k", band_vals, weight_dep[cols])
457 cout[spinup_row] = np.nan
458 return cout
461def extraction_to_deposition(
462 *,
463 cout: npt.ArrayLike,
464 flow: npt.ArrayLike,
465 tedges: pd.DatetimeIndex | np.ndarray,
466 cout_tedges: pd.DatetimeIndex | np.ndarray,
467 aquifer_pore_volume: float,
468 porosity: float,
469 thickness: float,
470 retardation_factor: float = 1.0,
471 regularization_strength: float = 1e-10,
472 spinup: str | None = "constant",
473) -> npt.NDArray[np.floating]:
474 """Compute deposition rates from concentration changes (deconvolution).
476 Inverts the forward model by solving ``W @ dep = cout`` where ``W`` is
477 the weight matrix from :func:`compute_deposition_weights`. Uses Tikhonov
478 regularization to smoothly blend data fitting with a physically motivated
479 target (transpose-and-normalize of the forward matrix).
481 Well-determined modes (large singular values relative to ``sqrt(λ)``) are
482 dominated by the data; poorly-determined modes are pulled toward the
483 target.
485 Parameters
486 ----------
487 cout : array-like
488 Concentration changes in extracted water [g/m³]. Length must equal
489 len(cout_tedges) - 1. May contain NaN values, which will be excluded
490 from the computation along with corresponding rows in the weight matrix.
491 The model assumes this value is constant over each interval
492 ``[cout_tedges[i], cout_tedges[i+1])``.
493 flow : array-like
494 Flow rates in aquifer [m³/day]. Length must equal len(tedges) - 1.
495 Must not contain NaN values. The model assumes this value is constant
496 over each interval ``[tedges[i], tedges[i+1])``.
497 tedges : pandas.DatetimeIndex
498 Time bin edges for deposition and flow data. Length must equal
499 len(flow) + 1.
500 cout_tedges : pandas.DatetimeIndex
501 Time bin edges for output concentration data. Length must equal
502 len(cout) + 1.
503 aquifer_pore_volume : float
504 Aquifer pore volume [m³].
505 porosity : float
506 Aquifer porosity [dimensionless].
507 thickness : float
508 Aquifer thickness [m].
509 retardation_factor : float, optional
510 Compound retardation factor, by default 1.0. Values > 1.0 indicate
511 slower transport due to sorption/interaction.
512 regularization_strength : float, optional
513 Tikhonov regularization parameter λ. Controls the tradeoff between
514 fitting the data (``||W dep - cout||²``) and staying close to the
515 regularization target (``λ ||dep - dep_target||²``). The target is
516 the transpose-and-normalize of the forward matrix applied to cout.
518 Larger values trust the target more (smoother, more biased); smaller
519 values trust the data more (noisier, less biased). Default is 1e-10.
520 spinup : {"constant"} | None, optional
521 Spin-up policy applied before building the forward weight matrix.
522 Default ``"constant"`` shifts ``tedges[0]`` backward by at least
523 ``retardation_factor * aquifer_pore_volume / flow[0]`` (rounded up to
524 whole bins, plus one extra bin) and treats flow as constant at its
525 first value over the prepended interval; the recovered deposition
526 vector is sliced back to the original ``tedges`` length so the public
527 output shape is unchanged. When the warm-start is undefined
528 (``flow[0] <= 0``, zero/NaN pore volume, non-positive first bin width,
529 or a pad exceeding the memory cap), ``"constant"`` silently reverts to
530 strict-validity behavior. ``None`` keeps strict-validity behavior. A
531 float raises ``NotImplementedError`` -- the fraction-threshold mode is
532 not implemented for deposition (matching the diffusion family).
534 Returns
535 -------
536 numpy.ndarray
537 Mean deposition rates [g/m²/day] between tedges. Length equals
538 len(tedges) - 1.
540 Raises
541 ------
542 ValueError
543 If input dimensions are incompatible, if flow contains NaN or negative
544 values, if ``retardation_factor`` is less than 1.0, or if physical
545 parameters are out of valid range (porosity not in (0, 1),
546 non-positive thickness or aquifer pore volume).
547 NotImplementedError
548 If ``spinup`` is anything other than ``None`` or ``"constant"`` (the
549 fraction-threshold mode is not implemented for deposition).
551 See Also
552 --------
553 deposition_to_extraction : Forward operation (convolution)
554 extraction_to_deposition_full : Full solver with nullspace options
555 spinup_duration : Earliest extraction time with a fully resolved deposition history
556 gwtransport.advection.extraction_to_infiltration : For concentration transport without deposition
557 gwtransport.utils.solve_inverse_transport_banded : Banded Tikhonov solver used for inversion
558 :ref:`concept-transport-equation` : Flow-weighted averaging approach
560 Notes
561 -----
562 This is a *source* term -- positive ``dep`` raises ``cout``. Sink
563 processes (pathogen attachment, first-order decay, particle filtration)
564 require the opposite sign convention and are not modelled here.
566 The forward model is ``W @ dep = cout``, where the weight matrix ``W``
567 encodes the physical relationship between deposition rates and
568 concentrations. ``W`` is genuinely banded -- row ``i`` is nonzero only on
569 the cin bins inside its residence-time window -- and is *built* in a
570 compact banded layout (peak memory ``O(n_cin * band)``); the Tikhonov
571 *solve* (:func:`gwtransport.utils.solve_inverse_transport_banded`)
572 transiently materializes the dense matrix and its Gram product
573 (``O(n_cout * n_cin + n_cin**2)``). Unlike advection (where rows sum to ~1), deposition
574 rows sum to ``r_i = residence_time_i / (retardation_factor * porosity *
575 thickness)``. Rows are
576 rescaled by ``r_i`` before solving: when ``W`` has full column rank and
577 ``cout`` lies in its column space this preserves the exact ``dep`` (with
578 the default square ``spinup="constant"`` system the warm-start padding
579 makes ``W`` structurally rank-deficient, so even a forward-generated
580 ``cout`` is recovered only up to the nullspace, regardless of
581 ``regularization_strength``), while for
582 overdetermined systems with noise it is equivalent to weighted least
583 squares with weights ``1 / r_i^2`` (shorter residence times get more
584 weight; under constant flow all ``r_i`` are equal and this reduces to
585 OLS). The rescaling puts the regularization target (transpose-and-normalize
586 of ``W`` applied to ``cout``) on the same scale as ``dep``, which controls
587 the regularization scale. Rows where the residence time cannot be computed
588 (spin-up period) and zero-flow cout bins are excluded automatically; NaN
589 values in ``cout`` are also excluded. The banded Tikhonov solve stays
590 well-defined via ``regularization_strength`` even when ``W`` is
591 rank-deficient (constant flow with integer ``RT/dt`` makes it a uniform
592 moving average with exact transfer-function zeros), so no rank-deficiency
593 warning is emitted.
595 Examples
596 --------
597 >>> import pandas as pd
598 >>> import numpy as np
599 >>> from gwtransport.deposition import extraction_to_deposition
600 >>>
601 >>> dates = pd.date_range("2020-01-01", "2020-01-10", freq="D")
602 >>> tedges = pd.date_range("2019-12-31 12:00", "2020-01-10 12:00", freq="D")
603 >>> cout_tedges = pd.date_range("2020-01-03 12:00", "2020-01-12 12:00", freq="D")
604 >>>
605 >>> flow = np.full(len(dates), 100.0) # m³/day
606 >>> cout = np.ones(len(cout_tedges) - 1) * 10.0 # g/m³
607 >>>
608 >>> dep = extraction_to_deposition(
609 ... cout=cout,
610 ... flow=flow,
611 ... tedges=tedges,
612 ... cout_tedges=cout_tedges,
613 ... aquifer_pore_volume=500.0,
614 ... porosity=0.3,
615 ... thickness=10.0,
616 ... )
617 >>> print(f"Deposition rates shape: {dep.shape}")
618 Deposition rates shape: (10,)
619 >>> print(f"Mean deposition rate: {np.nanmean(dep):.2f} g/m²/day")
620 Mean deposition rate: 6.00 g/m²/day
621 """
622 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
623 cout_values, flow_values = np.asarray(cout), np.asarray(flow)
625 _validate_deposition_inputs(
626 tedges=tedges,
627 flow_values=flow_values,
628 aquifer_pore_volume=aquifer_pore_volume,
629 porosity=porosity,
630 thickness=thickness,
631 retardation_factor=retardation_factor,
632 spinup=spinup,
633 cout_tedges=cout_tedges,
634 cout_values=cout_values,
635 )
637 # Apply spinup policy: optionally prepend warm-start bins to tedges/flow.
638 weight_tedges, weight_flow, _, _, n_pad = _resolve_spinup_inputs(
639 spinup,
640 tedges=tedges,
641 flow=flow_values,
642 aquifer_pore_volumes=np.array([aquifer_pore_volume]),
643 retardation_factor=retardation_factor,
644 )
646 # Build the banded forward operator (rows sum to r_k = RT_k/(porosity*thickness)).
647 band_vals, col_start, row_valid, _ = compute_deposition_weights(
648 flow=weight_flow,
649 tedges=weight_tedges,
650 cout_tedges=cout_tedges,
651 aquifer_pore_volume=aquifer_pore_volume,
652 porosity=porosity,
653 thickness=thickness,
654 retardation_factor=retardation_factor,
655 )
656 n_cin_padded = len(weight_tedges) - 1
658 # Per-row rescaling: normalize valid rows to sum 1 (w_norm = W_valid / r_k) and
659 # feed the banded solver observed = cout / r_k -- the SAME 1/r_k scaling, REQUIRED
660 # since deposition rows sum to r_k != 1 (unlike advection). Excluded rows -- NaN
661 # residence time / zero-flow cout bins (~row_valid) and NaN values in cout -- are
662 # zeroed in the band and given observed = 0 so they drop out of the normal
663 # equations (a zero band contributes nothing).
664 row_sums = band_vals.sum(axis=1)
665 keep = row_valid & ~np.isnan(cout_values)
666 safe_sum = np.where(keep, row_sums, 1.0)
667 band_norm = np.where(keep[:, None], band_vals / safe_sum[:, None], 0.0)
668 observed = np.where(keep, cout_values / safe_sum, 0.0)
670 # The banded Tikhonov solve is well-defined via regularization even when the
671 # operator is rank-deficient (constant flow with integer RT/dt makes it a
672 # uniform moving average with exact transfer-function zeros).
673 dep_padded = solve_inverse_transport_banded(
674 band_vals=band_norm,
675 col_start=col_start,
676 observed=observed,
677 n_output=n_cin_padded,
678 regularization_strength=regularization_strength,
679 )
680 # Drop warm-start prefix so output aligns with the user-provided tedges.
681 return dep_padded[n_pad:]
684def extraction_to_deposition_full(
685 *,
686 cout: npt.ArrayLike,
687 flow: npt.ArrayLike,
688 tedges: pd.DatetimeIndex | np.ndarray,
689 cout_tedges: pd.DatetimeIndex | np.ndarray,
690 aquifer_pore_volume: float,
691 porosity: float,
692 thickness: float,
693 retardation_factor: float = 1.0,
694 nullspace_objective: str | Callable = "squared_differences",
695 optimization_method: str = "BFGS",
696 rcond: float | None = None,
697 spinup: str | None = "constant",
698) -> npt.NDArray[np.floating]:
699 """Compute deposition rates from concentration changes using nullspace solver.
701 Full-featured inverse solver exposing all options of
702 :func:`~gwtransport.utils.solve_underdetermined_system`. For most use
703 cases, prefer :func:`extraction_to_deposition` which uses Tikhonov
704 regularization.
706 Parameters
707 ----------
708 cout : array-like
709 Concentration changes in extracted water [g/m³]. Length must equal
710 len(cout_tedges) - 1. May contain NaN values, which will be excluded
711 from the computation along with corresponding rows in the weight matrix.
712 flow : array-like
713 Flow rates in aquifer [m³/day]. Length must equal len(tedges) - 1.
714 Must not contain NaN values.
715 tedges : pandas.DatetimeIndex
716 Time bin edges for deposition and flow data. Length must equal
717 len(flow) + 1.
718 cout_tedges : pandas.DatetimeIndex
719 Time bin edges for output concentration data. Length must equal
720 len(cout) + 1.
721 aquifer_pore_volume : float
722 Aquifer pore volume [m³].
723 porosity : float
724 Aquifer porosity [dimensionless].
725 thickness : float
726 Aquifer thickness [m].
727 retardation_factor : float, optional
728 Compound retardation factor, by default 1.0.
729 nullspace_objective : str or callable, optional
730 Objective function to minimize in the nullspace. Options:
732 * ``"squared_differences"`` : Minimize sum of squared differences
733 between adjacent deposition rates (default, smooth solutions).
734 * ``"summed_differences"`` : Minimize sum of absolute differences
735 (sparse/piecewise constant solutions).
736 * callable : Custom objective ``f(coeffs, x_ls, nullspace_basis)``.
738 optimization_method : str, optional
739 Scipy optimization method. Default is ``"BFGS"``.
740 rcond : float or None, optional
741 Cutoff for small singular values in the least-squares step.
742 Default is None (uses numpy default).
743 spinup : {"constant"} | None, optional
744 Spin-up policy applied before building the forward weight matrix.
745 Default ``"constant"`` shifts ``tedges[0]`` backward by at least
746 ``retardation_factor * aquifer_pore_volume / flow[0]`` (rounded up to
747 whole bins, plus one extra bin), silently reverting to strict-validity
748 when the warm-start is undefined; the recovered deposition is sliced
749 back to the original ``tedges`` length. ``None`` keeps strict-validity
750 behavior. A float raises ``NotImplementedError`` -- the
751 fraction-threshold mode is not implemented for deposition (matching
752 the diffusion family). See :func:`extraction_to_deposition` for full
753 semantics.
755 Returns
756 -------
757 numpy.ndarray
758 Mean deposition rates [g/m²/day] between tedges. Length equals
759 len(tedges) - 1.
761 Raises
762 ------
763 ValueError
764 If cout_tedges does not have one more element than cout, if tedges
765 does not have one more element than flow, if flow contains NaN or
766 negative values, if ``retardation_factor`` is less than 1.0, or if
767 physical parameters are out of valid range (porosity not in (0, 1),
768 non-positive thickness or aquifer pore volume).
769 NotImplementedError
770 If ``spinup`` is anything other than ``None`` or ``"constant"`` (the
771 fraction-threshold mode is not implemented for deposition).
773 See Also
774 --------
775 extraction_to_deposition : Recommended solver using Tikhonov regularization.
776 spinup_duration : Earliest extraction time with a fully resolved deposition history.
777 gwtransport.utils.solve_underdetermined_system : Underlying solver.
779 Notes
780 -----
781 This is a *source* term -- positive ``dep`` raises ``cout``. Sink
782 processes (pathogen attachment, first-order decay, particle filtration)
783 require the opposite sign convention and are not modelled here.
784 """
785 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
786 cout_values, flow_values = np.asarray(cout), np.asarray(flow)
788 _validate_deposition_inputs(
789 tedges=tedges,
790 flow_values=flow_values,
791 aquifer_pore_volume=aquifer_pore_volume,
792 porosity=porosity,
793 thickness=thickness,
794 retardation_factor=retardation_factor,
795 spinup=spinup,
796 cout_tedges=cout_tedges,
797 cout_values=cout_values,
798 )
800 # Apply spinup policy: optionally prepend warm-start bins to tedges/flow.
801 weight_tedges, weight_flow, _, _, n_pad = _resolve_spinup_inputs(
802 spinup,
803 tedges=tedges,
804 flow=flow_values,
805 aquifer_pore_volumes=np.array([aquifer_pore_volume]),
806 retardation_factor=retardation_factor,
807 )
809 # The nullspace solver (lstsq + null_space SVD) genuinely needs a dense matrix,
810 # so build the band and densify it. Spin-up rows carry no complete deposition
811 # history and are therefore set entirely to NaN.
812 band_vals, col_start, _, spinup_row = compute_deposition_weights(
813 flow=weight_flow,
814 tedges=weight_tedges,
815 cout_tedges=cout_tedges,
816 aquifer_pore_volume=aquifer_pore_volume,
817 porosity=porosity,
818 thickness=thickness,
819 retardation_factor=retardation_factor,
820 )
821 n_cin_padded = len(weight_tedges) - 1
822 deposition_weights = _densify_weights(band_vals, col_start, n_cin_padded)
823 deposition_weights[spinup_row] = np.nan
825 dep_padded = solve_underdetermined_system(
826 coefficient_matrix=deposition_weights,
827 rhs_vector=cout_values,
828 nullspace_objective=nullspace_objective,
829 optimization_method=optimization_method,
830 rcond=rcond,
831 )
832 # Drop warm-start prefix so output aligns with the user-provided tedges.
833 return dep_padded[n_pad:]
836def spinup_duration(
837 *,
838 flow: npt.ArrayLike,
839 tedges: pd.DatetimeIndex,
840 aquifer_pore_volume: float,
841 retardation_factor: float = 1.0,
842) -> float:
843 """
844 Compute the spinup duration for deposition modeling.
846 The spinup duration is the smallest extraction time ``t*`` (relative to
847 ``tedges[0]``) at which the extracted water was infiltrated exactly at
848 ``tedges[0]``: equivalently, the time at which the cumulative flow first
849 reaches ``retardation_factor * aquifer_pore_volume``. For extraction times
850 earlier than ``t*`` the extracted concentration lacks complete deposition
851 history. Under constant flow this equals
852 ``aquifer_pore_volume * retardation_factor / flow``.
854 Parameters
855 ----------
856 flow : array-like
857 Flow rate of water in the aquifer [m³/day].
858 tedges : pandas.DatetimeIndex
859 Time edges for the flow data.
860 aquifer_pore_volume : float
861 Pore volume of the aquifer [m³].
862 retardation_factor : float, optional
863 Retardation factor of the compound in the aquifer [dimensionless], by
864 default 1.0.
866 Returns
867 -------
868 float
869 Spinup duration in days.
871 Raises
872 ------
873 ValueError
874 If the cumulative flow over the entire ``tedges`` window does not
875 reach ``retardation_factor * aquifer_pore_volume``, indicating the
876 flow timeseries is too short to characterise the spin-up duration.
878 See Also
879 --------
880 deposition_to_extraction : Forward solver that uses the spin-up duration to resolve NaN cout rows.
881 extraction_to_deposition : Inverse solver.
882 """
883 # Spin-up is the residence time of water *currently being extracted*: how
884 # far back in history we must know deposition to fully characterise the
885 # extracted concentration. For this boundary quantity -- the extraction
886 # time of the water infiltrated exactly at tedges[0] -- the two transport
887 # directions coincide: the infiltration->extraction and extraction->
888 # infiltration maps are inverses, so both solve V(t*) = R * V_p below.
889 #
890 # The smallest extraction time t* at which the extracted water was
891 # infiltrated exactly at tedges[0] satisfies
892 # ``flow_cum(t*) = R * V_pore``; the spin-up duration is then
893 # ``t* - 0 = t*``. Inverting the cumulative flow gives this value
894 # exactly (no quantisation to tedges spacing). Under constant flow
895 # this matches V*R/Q.
896 flow_arr = np.asarray(flow)
897 tedges_days = tedges_to_days(tedges)
898 dt_days = np.diff(tedges_days)
899 target_cum = retardation_factor * float(aquifer_pore_volume)
900 # Feasibility guard on the *un-bumped* cumulative total: the request is infeasible iff
901 # R*V_pore exceeds the true total infiltrated volume. (The monotone bump below would
902 # otherwise lift a trailing Q=0 plateau above target_cum and admit an infeasible request.)
903 flow_cum_raw = cumulative_flow_volume(flow_arr, dt_days)
904 if not flow_cum_raw[-1] >= target_cum:
905 msg = (
906 f"Cumulative flow over the entire tedges window ({flow_cum_raw[-1]:.6g} m³) does not reach "
907 f"retardation_factor * aquifer_pore_volume ({target_cum:.6g} m³); the flow timeseries is too "
908 "short to characterise the spin-up duration."
909 )
910 raise ValueError(msg)
911 # Plateaus in flow_cum from Q = 0 bins make V → t inversion multi-valued; bump duplicates
912 # by the smallest representable amount so np.interp resolves consistently at plateau levels.
913 # Reuse the raw cumsum (bit-identical to cumulative_flow_volume(..., strictly_monotone=True),
914 # which applies the same _make_strictly_monotone to the same array).
915 flow_cum = _make_strictly_monotone(flow_cum_raw)
916 return float(linear_interpolate(x_ref=flow_cum, y_ref=tedges_days, x_query=target_cum))