Coverage for src / gwtransport / deposition.py: 75%
138 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-27 06:32 +0000
1"""
2Deposition Analysis for 1D Aquifer Systems.
4This module models compound *source enrichment* in 1D groundwater flow: a deposition rate
5[ng/m²/day] supplies mass from the aquifer matrix into the water along the flow path,
6increasing the extracted concentration. The model is a *source* term (positive deposition
7adds mass to the water); it does NOT model removal processes such as pathogen attachment,
8particle filtration, or chemical precipitation, which would remove mass from the water and
9require an opposite sign convention. Example physical scenarios: dissolution of a sparingly
10soluble mineral coating from the matrix, leaching of a stored solute, mass release from a
11distributed contaminant source on the matrix surface. The module follows advection module
12patterns for consistency in forward (deposition to extraction) and reverse (extraction to
13deposition) calculations.
15Available functions:
17- :func:`deposition_to_extraction` - Compute concentrations from deposition rates (convolution).
18 Given deposition rate time series [ng/m²/day], computes resulting concentration changes in
19 extracted water [ng/m³]. Uses flow-weighted integration over contact area between water and
20 aquifer matrix. Accounts for aquifer geometry (porosity, thickness) and residence time
21 distribution.
23- :func:`extraction_to_deposition` - Compute deposition rates from concentration changes
24 (deconvolution). Given concentration change time series in extracted water [ng/m³], estimates
25 deposition rate history [ng/m²/day] that produced those changes. Uses Tikhonov regularization
26 toward a physically motivated target (transpose-and-normalize of the forward matrix). Handles
27 NaN values in concentration data by excluding corresponding time periods.
29- :func:`extraction_to_deposition_full` - Full-featured inverse solver exposing all options of
30 the nullspace-based solver (:func:`~gwtransport.utils.solve_underdetermined_system`). Allows
31 choosing between different nullspace objectives (``'squared_differences'``,
32 ``'summed_differences'``, or custom callables) and optimization methods.
34- :func:`compute_deposition_weights` - Internal helper function. Compute weight matrix relating
35 deposition rates to concentration changes. Used by both deposition_to_extraction (forward) and
36 extraction_to_deposition (reverse). Calculates contact area between water parcels and aquifer
37 matrix based on streamline geometry and residence times.
39- :func:`spinup_duration` - Compute spinup duration for deposition modeling. Returns the
40 earliest extraction time at which the extracted water was infiltrated at the start of the
41 flow series (equivalently, the time at which cumulative flow first reaches
42 ``retardation_factor * aquifer_pore_volume``). Before this duration the extracted
43 concentration lacks complete deposition history. Useful for determining the valid analysis
44 period and identifying when boundary effects are negligible.
46This file is part of gwtransport which is released under AGPL-3.0 license.
47See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
48"""
50import warnings
51from collections.abc import Callable
53import numpy as np
54import numpy.typing as npt
55import pandas as pd
57from gwtransport.deposition_utils import compute_average_heights
58from gwtransport.residence_time import residence_time
59from gwtransport.utils import compute_reverse_target, linear_interpolate, solve_tikhonov, solve_underdetermined_system
62def compute_deposition_weights(
63 *,
64 flow: npt.ArrayLike,
65 tedges: pd.DatetimeIndex,
66 cout_tedges: pd.DatetimeIndex,
67 aquifer_pore_volume: float,
68 porosity: float,
69 thickness: float,
70 retardation_factor: float = 1.0,
71) -> npt.NDArray[np.floating]:
72 """Compute deposition weights for concentration-deposition convolution.
74 Parameters
75 ----------
76 flow : array-like
77 Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1.
78 tedges : pandas.DatetimeIndex
79 Time bin edges for flow data.
80 cout_tedges : pandas.DatetimeIndex
81 Time bin edges for output concentration data.
82 aquifer_pore_volume : float
83 Aquifer pore volume [m3].
84 porosity : float
85 Aquifer porosity [dimensionless].
86 thickness : float
87 Aquifer thickness [m].
88 retardation_factor : float, optional
89 Compound retardation factor, by default 1.0.
91 Returns
92 -------
93 numpy.ndarray
94 Deposition weights matrix with shape (len(cout_tedges) - 1, len(tedges) - 1).
95 May contain NaN values where residence time cannot be computed.
97 Notes
98 -----
99 The returned weights matrix may contain NaN values in locations where the
100 residence time calculation fails or is undefined. This typically occurs
101 when flow conditions result in invalid or non-physical residence times.
102 """
103 # Convert to days relative to first time edge
104 t0 = tedges[0]
105 tedges_days = ((tedges - t0) / pd.Timedelta(days=1)).values
106 cout_tedges_days = ((cout_tedges - t0) / pd.Timedelta(days=1)).values
108 # Compute residence times and cumulative flow
109 flow_values = np.asarray(flow)
110 cout_rt_at_edges = residence_time(
111 flow=flow_values,
112 flow_tedges=tedges,
113 index=cout_tedges,
114 aquifer_pore_volumes=float(aquifer_pore_volume),
115 retardation_factor=retardation_factor,
116 direction="extraction_to_infiltration",
117 )
118 # residence_time returns shape (n_pore_volumes, n_index); we pass a single
119 # pore volume so select row 0 explicitly.
120 cout_tedges_days_infiltration = cout_tedges_days - cout_rt_at_edges[0, :]
122 flow_cum = np.concatenate(([0.0], np.cumsum(flow_values * np.diff(tedges_days))))
124 # Interpolate volumes at concentration time edges
125 start_vol = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days_infiltration)
126 end_vol = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days)
128 # Compute deposition weights
129 flow_cum_cout = flow_cum[None, :] - start_vol[:, None]
130 volume_array = compute_average_heights(
131 x_edges=tedges_days, y_edges=flow_cum_cout, y_lower=0.0, y_upper=retardation_factor * float(aquifer_pore_volume)
132 )
133 area_array = volume_array / (porosity * thickness)
134 extracted_volume = np.diff(end_vol)
135 # Zero-flow cout bins have extracted_volume == 0 (no water extracted), so
136 # the weight row collapses to 0: no extraction contributes no deposition
137 # signature to cout. Use a sentinel divisor on those rows to avoid the
138 # division-by-zero warning before np.where discards them.
139 numerator = area_array * np.diff(tedges_days)[None, :]
140 mask = extracted_volume > 0
141 safe_denom = np.where(mask, extracted_volume, 1.0)[:, None]
142 return np.where(mask[:, None], numerator / safe_denom, 0.0)
145def deposition_to_extraction(
146 *,
147 dep: npt.ArrayLike,
148 flow: npt.ArrayLike,
149 tedges: pd.DatetimeIndex | np.ndarray,
150 cout_tedges: pd.DatetimeIndex | np.ndarray,
151 aquifer_pore_volume: float,
152 porosity: float,
153 thickness: float,
154 retardation_factor: float = 1.0,
155) -> npt.NDArray[np.floating]:
156 """Compute concentrations from deposition rates (convolution).
158 Parameters
159 ----------
160 dep : array-like
161 Deposition rates [ng/m2/day]. Length must equal len(tedges) - 1.
162 flow : array-like
163 Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1. The model
164 assumes this value is constant over each interval ``[tedges[i], tedges[i+1])``.
165 tedges : pandas.DatetimeIndex
166 Time bin edges for deposition and flow data.
167 cout_tedges : pandas.DatetimeIndex
168 Time bin edges for output concentration data.
169 aquifer_pore_volume : float
170 Aquifer pore volume [m3].
171 porosity : float
172 Aquifer porosity [dimensionless].
173 thickness : float
174 Aquifer thickness [m].
175 retardation_factor : float, optional
176 Compound retardation factor, by default 1.0.
178 Returns
179 -------
180 numpy.ndarray
181 Concentration changes [ng/m3] with length len(cout_tedges) - 1.
183 Raises
184 ------
185 ValueError
186 If tedges does not have one more element than dep or flow, if input
187 arrays contain NaN values, or if physical parameters are out of
188 valid range (porosity not in (0, 1), non-positive thickness or
189 aquifer pore volume).
191 See Also
192 --------
193 extraction_to_deposition : Inverse operation (deconvolution)
194 gwtransport.advection.infiltration_to_extraction : For concentration transport without deposition
195 :ref:`concept-transport-equation` : Flow-weighted averaging approach
197 Examples
198 --------
199 >>> import pandas as pd
200 >>> import numpy as np
201 >>> from gwtransport.deposition import deposition_to_extraction
202 >>> dates = pd.date_range("2020-01-01", "2020-01-10", freq="D")
203 >>> tedges = pd.date_range("2019-12-31 12:00", "2020-01-10 12:00", freq="D")
204 >>> cout_tedges = pd.date_range("2020-01-03 12:00", "2020-01-12 12:00", freq="D")
205 >>> dep = np.ones(len(dates))
206 >>> flow = np.full(len(dates), 100.0)
207 >>> cout = deposition_to_extraction(
208 ... dep=dep,
209 ... flow=flow,
210 ... tedges=tedges,
211 ... cout_tedges=cout_tedges,
212 ... aquifer_pore_volume=500.0,
213 ... porosity=0.3,
214 ... thickness=10.0,
215 ... )
216 """
217 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
218 dep_values, flow_values = np.asarray(dep), np.asarray(flow)
220 # Validate input dimensions and values
221 if len(tedges) != len(dep_values) + 1:
222 msg = "tedges must have one more element than dep"
223 raise ValueError(msg)
224 if len(tedges) != len(flow_values) + 1:
225 msg = "tedges must have one more element than flow"
226 raise ValueError(msg)
227 if np.any(np.isnan(dep_values)) or np.any(np.isnan(flow_values)):
228 msg = "Input arrays cannot contain NaN values"
229 raise ValueError(msg)
230 if np.any(flow_values < 0):
231 msg = "flow must be non-negative (negative flow not supported)"
232 raise ValueError(msg)
234 # Validate physical parameters
235 if not 0 < porosity < 1:
236 msg = f"Porosity must be in (0, 1), got {porosity}"
237 raise ValueError(msg)
238 if thickness <= 0:
239 msg = f"Thickness must be positive, got {thickness}"
240 raise ValueError(msg)
241 if aquifer_pore_volume <= 0:
242 msg = f"Aquifer pore volume must be positive, got {aquifer_pore_volume}"
243 raise ValueError(msg)
245 # Compute deposition weights
246 deposition_weights = compute_deposition_weights(
247 flow=flow_values,
248 tedges=tedges,
249 cout_tedges=cout_tedges,
250 aquifer_pore_volume=aquifer_pore_volume,
251 porosity=porosity,
252 thickness=thickness,
253 retardation_factor=retardation_factor,
254 )
256 return deposition_weights.dot(dep_values)
259def extraction_to_deposition(
260 *,
261 cout: npt.ArrayLike,
262 flow: npt.ArrayLike,
263 tedges: pd.DatetimeIndex | np.ndarray,
264 cout_tedges: pd.DatetimeIndex | np.ndarray,
265 aquifer_pore_volume: float,
266 porosity: float,
267 thickness: float,
268 retardation_factor: float = 1.0,
269 regularization_strength: float = 1e-10,
270) -> npt.NDArray[np.floating]:
271 """Compute deposition rates from concentration changes (deconvolution).
273 Inverts the forward model by solving ``W @ dep = cout`` where ``W`` is
274 the weight matrix from :func:`compute_deposition_weights`. Uses Tikhonov
275 regularization to smoothly blend data fitting with a physically motivated
276 target (transpose-and-normalize of the forward matrix).
278 Well-determined modes (large singular values relative to ``sqrt(λ)``) are
279 dominated by the data; poorly-determined modes are pulled toward the
280 target.
282 Parameters
283 ----------
284 cout : array-like
285 Concentration changes in extracted water [ng/m3]. Length must equal
286 len(cout_tedges) - 1. May contain NaN values, which will be excluded
287 from the computation along with corresponding rows in the weight matrix.
288 The model assumes this value is constant over each interval
289 ``[cout_tedges[i], cout_tedges[i+1])``.
290 flow : array-like
291 Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1.
292 Must not contain NaN values. The model assumes this value is constant
293 over each interval ``[tedges[i], tedges[i+1])``.
294 tedges : pandas.DatetimeIndex
295 Time bin edges for deposition and flow data. Length must equal
296 len(flow) + 1.
297 cout_tedges : pandas.DatetimeIndex
298 Time bin edges for output concentration data. Length must equal
299 len(cout) + 1.
300 aquifer_pore_volume : float
301 Aquifer pore volume [m3].
302 porosity : float
303 Aquifer porosity [dimensionless].
304 thickness : float
305 Aquifer thickness [m].
306 retardation_factor : float, optional
307 Compound retardation factor, by default 1.0. Values > 1.0 indicate
308 slower transport due to sorption/interaction.
309 regularization_strength : float, optional
310 Tikhonov regularization parameter λ. Controls the tradeoff between
311 fitting the data (``||W dep - cout||²``) and staying close to the
312 regularization target (``λ ||dep - dep_target||²``). The target is
313 the transpose-and-normalize of the forward matrix applied to cout.
315 Larger values trust the target more (smoother, more biased); smaller
316 values trust the data more (noisier, less biased). Default is 1e-10.
318 Returns
319 -------
320 numpy.ndarray
321 Mean deposition rates [ng/m2/day] between tedges. Length equals
322 len(tedges) - 1.
324 Raises
325 ------
326 ValueError
327 If input dimensions are incompatible or if flow contains NaN values.
329 Warns
330 -----
331 UserWarning
332 When the weight matrix is rank-deficient. This occurs with constant
333 flow when the residence time is an integer multiple of the time step
334 width. The deposition weight matrix then acts as a uniform moving
335 average whose transfer function has exact zeros, making certain
336 deposition patterns invisible in the concentration signal. To fix,
337 adjust ``aquifer_pore_volume`` slightly (e.g., multiply by 1.001).
339 See Also
340 --------
341 deposition_to_extraction : Forward operation (convolution)
342 extraction_to_deposition_full : Full solver with nullspace options
343 gwtransport.advection.extraction_to_infiltration : For concentration transport without deposition
344 gwtransport.utils.solve_tikhonov : Solver used for inversion
345 :ref:`concept-transport-equation` : Flow-weighted averaging approach
347 Notes
348 -----
349 The forward model is ``W @ dep = cout``, where the weight matrix ``W``
350 encodes the physical relationship between deposition rates and
351 concentrations. Unlike advection (where rows sum to ~1), deposition rows
352 sum to ``residence_time / (porosity * thickness)``. The system is
353 row-normalized before solving so that each observation contributes equally
354 and ``compute_reverse_target`` gives the correct scale for the
355 regularization target. Rows where the residence time cannot be computed
356 (spin-up period) contain NaN and are excluded automatically. NaN values
357 in ``cout`` are also excluded.
359 Examples
360 --------
361 >>> import pandas as pd
362 >>> import numpy as np
363 >>> from gwtransport.deposition import extraction_to_deposition
364 >>>
365 >>> dates = pd.date_range("2020-01-01", "2020-01-10", freq="D")
366 >>> tedges = pd.date_range("2019-12-31 12:00", "2020-01-10 12:00", freq="D")
367 >>> cout_tedges = pd.date_range("2020-01-03 12:00", "2020-01-12 12:00", freq="D")
368 >>>
369 >>> flow = np.full(len(dates), 100.0) # m3/day
370 >>> cout = np.ones(len(cout_tedges) - 1) * 10.0 # ng/m3
371 >>>
372 >>> dep = extraction_to_deposition(
373 ... cout=cout,
374 ... flow=flow,
375 ... tedges=tedges,
376 ... cout_tedges=cout_tedges,
377 ... aquifer_pore_volume=500.0,
378 ... porosity=0.3,
379 ... thickness=10.0,
380 ... )
381 >>> print(f"Deposition rates shape: {dep.shape}")
382 Deposition rates shape: (10,)
383 >>> print(f"Mean deposition rate: {np.nanmean(dep):.2f} ng/m2/day")
384 Mean deposition rate: 6.00 ng/m2/day
385 """
386 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
387 cout_values, flow_values = np.asarray(cout), np.asarray(flow)
389 # Validate input dimensions and values
390 if len(cout_tedges) != len(cout_values) + 1:
391 msg = "cout_tedges must have one more element than cout"
392 raise ValueError(msg)
393 if len(tedges) != len(flow_values) + 1:
394 msg = "tedges must have one more element than flow"
395 raise ValueError(msg)
396 if np.any(np.isnan(flow_values)):
397 msg = "flow array cannot contain NaN values"
398 raise ValueError(msg)
399 if np.any(flow_values < 0):
400 msg = "flow must be non-negative (negative flow not supported)"
401 raise ValueError(msg)
403 # Validate physical parameters
404 if not 0 < porosity < 1:
405 msg = f"Porosity must be in (0, 1), got {porosity}"
406 raise ValueError(msg)
407 if thickness <= 0:
408 msg = f"Thickness must be positive, got {thickness}"
409 raise ValueError(msg)
410 if aquifer_pore_volume <= 0:
411 msg = f"Aquifer pore volume must be positive, got {aquifer_pore_volume}"
412 raise ValueError(msg)
414 # Compute deposition weights
415 deposition_weights = compute_deposition_weights(
416 flow=flow_values,
417 tedges=tedges,
418 cout_tedges=cout_tedges,
419 aquifer_pore_volume=aquifer_pore_volume,
420 porosity=porosity,
421 thickness=thickness,
422 retardation_factor=retardation_factor,
423 )
425 n_dep = len(tedges) - 1
427 # Normalize weight matrix rows to sum to 1. The deposition weight matrix
428 # W has row sums equal to residence_time/(porosity*thickness), not 1 like
429 # advection. Normalizing makes each observation equally important and
430 # gives compute_reverse_target the correct scale for the target.
431 # NaN rows (spin-up), all-zero rows (zero-flow cout bins; carry no info),
432 # and NaN values in cout are excluded by solve_tikhonov.
433 valid_rows = ~np.isnan(deposition_weights).any(axis=1) & (np.abs(deposition_weights).sum(axis=1) > 0)
434 valid_weights = deposition_weights[valid_rows]
435 row_sums = valid_weights.sum(axis=1, keepdims=True)
436 col_active = np.sum(np.abs(valid_weights), axis=0) > 0
438 if not np.any(col_active):
439 return np.full(n_dep, np.nan)
441 # Build normalized system: W_norm @ dep = cout_norm
442 w_norm = valid_weights / row_sums
443 cout_norm = cout_values[valid_rows] / row_sums.ravel()
445 # Check for rank deficiency. With constant flow and integer RT/dt, the
446 # deposition weight matrix acts as a uniform moving average whose transfer
447 # function has exact zeros, making certain deposition patterns invisible.
448 n_active = int(col_active.sum())
449 rank = np.linalg.matrix_rank(w_norm[:, col_active])
450 if rank < n_active:
451 warnings.warn(
452 f"Weight matrix is rank-deficient (rank {rank} < {n_active} active "
453 f"columns). This occurs with constant flow when the residence time "
454 f"is an integer multiple of the time step width, creating deposition "
455 f"patterns that are invisible in the concentration signal. The "
456 f"underdetermined modes will be pulled toward the regularization "
457 f"target instead of being determined by data. To achieve full rank, "
458 f"adjust aquifer_pore_volume slightly (e.g., multiply by 1.001).",
459 stacklevel=2,
460 )
462 # Reconstruct full arrays with NaN for invalid rows
463 full_w_norm = np.full_like(deposition_weights, np.nan)
464 full_w_norm[valid_rows] = w_norm
465 full_cout_norm = np.full(len(cout_values), np.nan)
466 full_cout_norm[valid_rows] = cout_norm
468 x_target = compute_reverse_target(coeff_matrix=w_norm, rhs_vector=cout_norm)
470 dep_solved = solve_tikhonov(
471 coefficient_matrix=full_w_norm,
472 rhs_vector=full_cout_norm,
473 x_target=x_target,
474 regularization_strength=regularization_strength,
475 )
477 out = np.full(n_dep, np.nan)
478 out[col_active] = dep_solved[col_active]
479 return out
482def extraction_to_deposition_full(
483 *,
484 cout: npt.ArrayLike,
485 flow: npt.ArrayLike,
486 tedges: pd.DatetimeIndex | np.ndarray,
487 cout_tedges: pd.DatetimeIndex | np.ndarray,
488 aquifer_pore_volume: float,
489 porosity: float,
490 thickness: float,
491 retardation_factor: float = 1.0,
492 nullspace_objective: str | Callable = "squared_differences",
493 optimization_method: str = "BFGS",
494 rcond: float | None = None,
495 x_target: npt.NDArray[np.floating] | None = None,
496) -> npt.NDArray[np.floating]:
497 """Compute deposition rates from concentration changes using nullspace solver.
499 Full-featured inverse solver exposing all options of
500 :func:`~gwtransport.utils.solve_underdetermined_system`. For most use
501 cases, prefer :func:`extraction_to_deposition` which uses Tikhonov
502 regularization.
504 Parameters
505 ----------
506 cout : array-like
507 Concentration changes in extracted water [ng/m3]. Length must equal
508 len(cout_tedges) - 1. May contain NaN values, which will be excluded
509 from the computation along with corresponding rows in the weight matrix.
510 flow : array-like
511 Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1.
512 Must not contain NaN values.
513 tedges : pandas.DatetimeIndex
514 Time bin edges for deposition and flow data. Length must equal
515 len(flow) + 1.
516 cout_tedges : pandas.DatetimeIndex
517 Time bin edges for output concentration data. Length must equal
518 len(cout) + 1.
519 aquifer_pore_volume : float
520 Aquifer pore volume [m3].
521 porosity : float
522 Aquifer porosity [dimensionless].
523 thickness : float
524 Aquifer thickness [m].
525 retardation_factor : float, optional
526 Compound retardation factor, by default 1.0.
527 nullspace_objective : str or callable, optional
528 Objective function to minimize in the nullspace. Options:
530 * ``"squared_differences"`` : Minimize sum of squared differences
531 between adjacent deposition rates (default, smooth solutions).
532 * ``"summed_differences"`` : Minimize sum of absolute differences
533 (sparse/piecewise constant solutions).
534 * callable : Custom objective ``f(coeffs, x_ls, nullspace_basis)``.
536 optimization_method : str, optional
537 Scipy optimization method. Default is ``"BFGS"``.
538 rcond : float or None, optional
539 Cutoff for small singular values in the least-squares step.
540 Default is None (uses numpy default).
541 x_target : numpy.ndarray or None, optional
542 Optional target solution for the nullspace optimization.
543 Default is None.
545 Returns
546 -------
547 numpy.ndarray
548 Mean deposition rates [ng/m2/day] between tedges. Length equals
549 len(tedges) - 1.
551 Raises
552 ------
553 ValueError
554 If cout_tedges does not have one more element than cout, if tedges
555 does not have one more element than flow, if flow contains NaN
556 values, or if physical parameters are out of valid range (porosity
557 not in (0, 1), non-positive thickness or aquifer pore volume).
559 See Also
560 --------
561 extraction_to_deposition : Recommended solver using Tikhonov regularization.
562 gwtransport.utils.solve_underdetermined_system : Underlying solver.
563 """
564 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
565 cout_values, flow_values = np.asarray(cout), np.asarray(flow)
567 # Validate input dimensions and values
568 if len(cout_tedges) != len(cout_values) + 1:
569 msg = "cout_tedges must have one more element than cout"
570 raise ValueError(msg)
571 if len(tedges) != len(flow_values) + 1:
572 msg = "tedges must have one more element than flow"
573 raise ValueError(msg)
574 if np.any(np.isnan(flow_values)):
575 msg = "flow array cannot contain NaN values"
576 raise ValueError(msg)
577 if np.any(flow_values < 0):
578 msg = "flow must be non-negative (negative flow not supported)"
579 raise ValueError(msg)
581 # Validate physical parameters
582 if not 0 < porosity < 1:
583 msg = f"Porosity must be in (0, 1), got {porosity}"
584 raise ValueError(msg)
585 if thickness <= 0:
586 msg = f"Thickness must be positive, got {thickness}"
587 raise ValueError(msg)
588 if aquifer_pore_volume <= 0:
589 msg = f"Aquifer pore volume must be positive, got {aquifer_pore_volume}"
590 raise ValueError(msg)
592 # Compute deposition weights
593 deposition_weights = compute_deposition_weights(
594 flow=flow_values,
595 tedges=tedges,
596 cout_tedges=cout_tedges,
597 aquifer_pore_volume=aquifer_pore_volume,
598 porosity=porosity,
599 thickness=thickness,
600 retardation_factor=retardation_factor,
601 )
603 return solve_underdetermined_system(
604 coefficient_matrix=deposition_weights,
605 rhs_vector=cout_values,
606 nullspace_objective=nullspace_objective,
607 optimization_method=optimization_method,
608 rcond=rcond,
609 x_target=x_target,
610 )
613def spinup_duration(
614 *,
615 flow: np.ndarray,
616 flow_tedges: pd.DatetimeIndex,
617 aquifer_pore_volume: float,
618 retardation_factor: float,
619) -> float:
620 """
621 Compute the spinup duration for deposition modeling.
623 The spinup duration is the smallest extraction time ``t*`` (relative to
624 ``flow_tedges[0]``) at which the extracted water was infiltrated exactly
625 at ``flow_tedges[0]``: equivalently, the time at which the cumulative
626 flow first reaches ``retardation_factor * aquifer_pore_volume``. For
627 extraction times earlier than ``t*`` the extracted concentration lacks
628 complete deposition history. Under constant flow this equals
629 ``aquifer_pore_volume * retardation_factor / flow``.
631 Parameters
632 ----------
633 flow : numpy.ndarray
634 Flow rate of water in the aquifer [m3/day].
635 flow_tedges : pandas.DatetimeIndex
636 Time edges for the flow data.
637 aquifer_pore_volume : float
638 Pore volume of the aquifer [m3].
639 retardation_factor : float
640 Retardation factor of the compound in the aquifer [dimensionless].
642 Returns
643 -------
644 float
645 Spinup duration in days.
647 Raises
648 ------
649 ValueError
650 If the cumulative flow over the entire ``flow_tedges`` window does
651 not reach ``retardation_factor * aquifer_pore_volume``, indicating
652 the flow timeseries is too short to fully characterise the aquifer.
653 """
654 # Spin-up is the residence time of water *currently being extracted*: how
655 # far back in history we must know deposition to fully characterise the
656 # extracted concentration. This uses the ``extraction_to_infiltration``
657 # direction. Under variable flow this differs from
658 # ``infiltration_to_extraction`` (which would describe how long ahead
659 # water infiltrated at the first time step will take to be extracted, a
660 # forward-in-time question that is not what spin-up means).
661 #
662 # The smallest extraction time t* at which the extracted water was
663 # infiltrated exactly at flow_tedges[0] satisfies
664 # ``flow_cum(t*) = R * V_pore``; the spin-up duration is then
665 # ``t* - 0 = t*``. Inverting the cumulative flow gives this value
666 # exactly (no quantisation to flow_tedges spacing). Under constant flow
667 # this matches V*R/Q.
668 flow_arr = np.asarray(flow)
669 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
670 flow_cum = np.concatenate(([0.0], np.cumsum(flow_arr * np.diff(flow_tedges_days))))
671 target_cum = retardation_factor * float(aquifer_pore_volume)
672 if not flow_cum[-1] >= target_cum:
673 msg = "Residence time at the first time step is NaN. This indicates that the aquifer is not fully informed: flow timeseries too short."
674 raise ValueError(msg)
675 rt_value = float(linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=target_cum))
676 if np.isnan(rt_value):
677 msg = "Residence time at the first time step is NaN. This indicates that the aquifer is not fully informed: flow timeseries too short."
678 raise ValueError(msg)
679 return rt_value