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