Coverage for src / gwtransport / residence_time.py: 91%
91 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-27 06:32 +0000
1"""
2Residence Time Calculations for Retarded Compound Transport.
4This module provides functions to compute residence times for compounds traveling through
5aquifer systems, accounting for flow variability, pore volume, and retardation due to
6physical or chemical interactions with the aquifer matrix. Residence time represents the
7duration a compound spends traveling from infiltration to extraction points, depending on
8flow rate (higher flow yields shorter residence time), pore volume (larger volume yields
9longer residence time), and retardation factor (interaction with matrix yields longer
10residence time).
12Available functions:
14- :func:`residence_time` - Compute residence times at specific time indices. Supports both
15 forward (infiltration to extraction) and reverse (extraction to infiltration) directions.
16 Handles single or multiple pore volumes (2D output for multiple volumes). Returns residence
17 times in days using cumulative flow integration for accurate time-varying flow handling.
19- :func:`residence_time_mean` - Compute mean residence times over time intervals. Calculates
20 average residence time between specified time edges using linear averaging to properly weight
21 time-varying residence times. Supports same directional options as residence_time. Particularly
22 useful for time-binned analysis.
24- :func:`fraction_explained` - Compute fraction of aquifer pore volumes with valid residence
25 times. Indicates how many pore volumes have sufficient flow history to compute residence time.
26 Returns values in [0, 1] where 1.0 means all volumes are fully informed. Useful for assessing
27 spin-up periods and data coverage. NaN residence times indicate insufficient flow history.
29This file is part of gwtransport which is released under AGPL-3.0 license.
30See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
31"""
33import numpy as np
34import numpy.typing as npt
35import pandas as pd
37from gwtransport.utils import linear_average, linear_interpolate
40def residence_time(
41 *,
42 flow: npt.ArrayLike,
43 flow_tedges: pd.DatetimeIndex | np.ndarray,
44 aquifer_pore_volumes: npt.ArrayLike,
45 index: pd.DatetimeIndex | np.ndarray | None = None,
46 direction: str = "extraction_to_infiltration",
47 retardation_factor: float = 1.0,
48) -> npt.NDArray[np.floating]:
49 """
50 Compute the residence time of retarded compound in the water in the aquifer.
52 Parameters
53 ----------
54 flow : array-like
55 Flow rate of water in the aquifer [m3/day]. The length of `flow` should match the length of `flow_tedges` minus one.
56 flow_tedges : pandas.DatetimeIndex
57 Time edges for the flow data. Used to compute the cumulative flow.
58 Has a length of one more than `flow`.
59 aquifer_pore_volumes : float or array-like of float
60 Pore volume(s) of the aquifer [m3]. Can be a single value or an array
61 of pore volumes representing different flow paths.
62 index : pandas.DatetimeIndex, optional
63 Index at which to compute the residence time. If left to None, flow_tedges is used.
64 Default is None.
65 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
66 Direction of the flow calculation:
67 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
68 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
69 Default is 'extraction_to_infiltration'.
70 retardation_factor : float, optional
71 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
73 Returns
74 -------
75 numpy.ndarray
76 Residence time of the retarded compound in the aquifer [days].
78 Raises
79 ------
80 ValueError
81 If ``flow_tedges`` does not have exactly one more element than ``flow``.
82 If ``direction`` is not ``'extraction_to_infiltration'`` or
83 ``'infiltration_to_extraction'``.
85 See Also
86 --------
87 residence_time_mean : Compute mean residence time over time intervals
88 gwtransport.advection.gamma_infiltration_to_extraction : Use residence times for transport
89 gwtransport.logremoval.residence_time_to_log_removal : Convert residence time to log removal
90 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
91 :ref:`concept-retardation-factor` : Slower movement due to sorption
92 """
93 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes)
94 flow_tedges = pd.DatetimeIndex(flow_tedges)
96 # Convert to arrays for type safety
97 flow = np.asarray(flow)
98 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes)
100 if len(flow_tedges) != len(flow) + 1:
101 msg = "tedges must have one more element than flow"
102 raise ValueError(msg)
104 # Check for negative flow values - physically invalid
105 if np.any(flow < 0):
106 # Return NaN array with correct shape
107 n_output = len(flow_tedges) - 1 if index is None else len(index)
108 n_pore_volumes = len(aquifer_pore_volumes)
109 return np.full((n_pore_volumes, n_output), np.nan)
111 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
112 flow_tdelta = np.diff(flow_tedges_days)
113 flow_cum = np.concatenate((
114 [0.0],
115 (flow * flow_tdelta).cumsum(),
116 )) # at flow_tedges and flow_tedges_days. First value is 0.
118 if index is None:
119 # If index is not provided return the residence time that matches with the index of the flow; at the center of the flow bin.
120 index_dates_days_extraction = (flow_tedges_days[:-1] + flow_tedges_days[1:]) / 2
121 flow_cum_at_index = (flow_cum[:-1] + flow_cum[1:]) / 2 # at the center of the flow bin
122 else:
123 index_dates_days_extraction = np.asarray((index - flow_tedges[0]) / np.timedelta64(1, "D"))
124 flow_cum_at_index = linear_interpolate(
125 x_ref=flow_tedges_days, y_ref=flow_cum, x_query=index_dates_days_extraction, left=np.nan, right=np.nan
126 )
128 if direction == "extraction_to_infiltration":
129 # How many days ago was the extraced water infiltrated
130 a = flow_cum_at_index[None, :] - retardation_factor * aquifer_pore_volumes[:, None]
131 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
132 data = index_dates_days_extraction - days
133 elif direction == "infiltration_to_extraction":
134 # In how many days the water that is infiltrated now be extracted
135 a = flow_cum_at_index[None, :] + retardation_factor * aquifer_pore_volumes[:, None]
136 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
137 data = days - index_dates_days_extraction
138 else:
139 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
140 raise ValueError(msg)
142 return data
145def residence_time_mean(
146 *,
147 flow: npt.ArrayLike,
148 flow_tedges: pd.DatetimeIndex | np.ndarray,
149 tedges_out: pd.DatetimeIndex | np.ndarray,
150 aquifer_pore_volumes: npt.ArrayLike,
151 direction: str = "extraction_to_infiltration",
152 retardation_factor: float = 1.0,
153) -> npt.NDArray[np.floating]:
154 """
155 Compute the mean residence time of a retarded compound in the aquifer between specified time edges.
157 This function calculates the average residence time of a retarded compound in the aquifer
158 between specified time intervals. It can compute both extraction to infiltration modeling (extraction direction:
159 when was extracted water infiltrated) and infiltration to extraction modeling (infiltration direction: when will
160 infiltrated water be extracted).
162 The function handles time series data by computing the cumulative flow and using linear
163 interpolation and averaging to determine mean residence times between the specified time edges.
165 Parameters
166 ----------
167 flow : array-like
168 Flow rate of water in the aquifer [m3/day]. Should be an array of flow values
169 corresponding to the intervals defined by flow_tedges.
170 flow_tedges : array-like
171 Time edges for the flow data, as datetime64 objects. These define the time
172 intervals for which the flow values are provided.
173 tedges_out : array-like
174 Output time edges as datetime64 objects. These define the intervals for which
175 the mean residence times will be calculated.
176 aquifer_pore_volumes : float or array-like
177 Pore volume(s) of the aquifer [m3]. Can be a single value or an array of values
178 for multiple pore volume scenarios.
179 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
180 Direction of the flow calculation:
181 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
182 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
183 Default is 'extraction_to_infiltration'.
184 retardation_factor : float, optional
185 Retardation factor of the compound in the aquifer [dimensionless].
186 A value greater than 1.0 indicates that the compound moves slower than water.
187 Default is 1.0 (no retardation).
189 Returns
190 -------
191 numpy.ndarray
192 Mean residence time of the retarded compound in the aquifer [days] for each interval
193 defined by tedges_out. The first dimension corresponds to the different pore volumes
194 and the second to the residence times between tedges_out.
196 Raises
197 ------
198 ValueError
199 If ``direction`` is not ``'extraction_to_infiltration'`` or
200 ``'infiltration_to_extraction'``.
202 See Also
203 --------
204 residence_time : Compute residence time at specific time indices
205 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
207 Notes
208 -----
209 - The function converts datetime objects to days since the start of the time series.
210 - For extraction_to_infiltration direction, the function computes how many days ago water was infiltrated.
211 - For infiltration_to_extraction direction, the function computes how many days until water will be extracted.
212 - The function uses linear interpolation for computing residence times at specific points
213 and linear averaging for computing mean values over intervals.
215 Examples
216 --------
217 >>> import pandas as pd
218 >>> import numpy as np
219 >>> from gwtransport.residence_time import residence_time_mean
220 >>> # Create sample flow data
221 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
222 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day
223 >>> pore_volume = 200.0 # Aquifer pore volume in m³
224 >>> # Calculate mean residence times
225 >>> mean_times = residence_time_mean(
226 ... flow=flow_values,
227 ... flow_tedges=flow_dates,
228 ... tedges_out=flow_dates,
229 ... aquifer_pore_volumes=pore_volume,
230 ... direction="extraction_to_infiltration",
231 ... )
232 >>> # With constant flow of 100 m³/day and pore volume of 200 m³,
233 >>> # mean residence time should be approximately 2 days
234 >>> print(mean_times) # doctest: +NORMALIZE_WHITESPACE
235 [[nan nan 2. 2. 2. 2. 2. 2. 2.]]
236 """
237 flow = np.asarray(flow)
238 flow_tedges = pd.DatetimeIndex(flow_tedges)
239 tedges_out = pd.DatetimeIndex(tedges_out)
240 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes)
242 # Check for negative flow values - physically invalid
243 if np.any(flow < 0):
244 # Return NaN array with correct shape
245 n_pore_volumes = len(aquifer_pore_volumes)
246 n_output_bins = len(tedges_out) - 1
247 return np.full((n_pore_volumes, n_output_bins), np.nan)
249 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
250 tedges_out_days = np.asarray((tedges_out - flow_tedges[0]) / np.timedelta64(1, "D"))
252 # compute cumulative flow at flow_tedges
253 flow_cum = np.concatenate(([0.0], np.cumsum(flow * np.diff(flow_tedges_days))))
255 if direction == "extraction_to_infiltration":
256 # How many days ago was the extraced water infiltrated
257 a = flow_cum[None, :] - retardation_factor * aquifer_pore_volumes[:, None]
258 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
259 data_edges = flow_tedges_days - days
260 elif direction == "infiltration_to_extraction":
261 # In how many days the water that is infiltrated now be extracted
262 a = flow_cum[None, :] + retardation_factor * aquifer_pore_volumes[:, None]
263 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
264 data_edges = days - flow_tedges_days
265 else:
266 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
267 raise ValueError(msg)
269 # Vectorized linear average across all pore volumes in one call. ``data_edges`` is 2D
270 # with shape (n_pore_volumes, n_flow_edges); each row shares the same x_data and the
271 # same x_edges. linear_average's 2D-y_data path handles per-row NaN segments cleanly.
272 return linear_average(x_data=flow_tedges_days, y_data=data_edges, x_edges=tedges_out_days)
275def fraction_explained(
276 *,
277 rt: npt.NDArray[np.floating] | None = None,
278 flow: npt.ArrayLike | None = None,
279 flow_tedges: pd.DatetimeIndex | np.ndarray | None = None,
280 aquifer_pore_volumes: npt.ArrayLike | None = None,
281 index: pd.DatetimeIndex | np.ndarray | None = None,
282 direction: str = "extraction_to_infiltration",
283 retardation_factor: float = 1.0,
284) -> npt.NDArray[np.floating]:
285 """
286 Compute the fraction of the aquifer that is informed with respect to the retarded flow.
288 Parameters
289 ----------
290 rt : numpy.ndarray, optional
291 Pre-computed residence time array [days]. If not provided, it will be computed.
292 flow : array-like, optional
293 Flow rate of water in the aquifer [m3/day]. The length of `flow` should match the length of `flow_tedges` minus one.
294 flow_tedges : pandas.DatetimeIndex, optional
295 Time edges for the flow data. Used to compute the cumulative flow.
296 Has a length of one more than `flow`. Inbetween neighboring time edges, the flow is assumed constant.
297 aquifer_pore_volumes : float or array-like of float, optional
298 Pore volume(s) of the aquifer [m3]. Can be a single value or an array
299 of pore volumes representing different flow paths.
300 index : pandas.DatetimeIndex, optional
301 Index at which to compute the fraction. If left to None, the index of `flow` is used.
302 Default is None.
303 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
304 Direction of the flow calculation:
305 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
306 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
307 Default is 'extraction_to_infiltration'.
308 retardation_factor : float, optional
309 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
311 Returns
312 -------
313 numpy.ndarray
314 Fraction of the aquifer that is informed with respect to the retarded flow.
316 Raises
317 ------
318 ValueError
319 If ``rt`` is not provided and any of ``flow``, ``flow_tedges``, or
320 ``aquifer_pore_volumes`` are missing. If ``rt`` is provided but is not 2D.
321 """
322 if rt is None:
323 # Validate that required parameters are provided for computing rt
324 if flow is None:
325 msg = "Either rt or flow must be provided"
326 raise ValueError(msg)
327 if flow_tedges is None:
328 msg = "Either rt or flow_tedges must be provided"
329 raise ValueError(msg)
330 if aquifer_pore_volumes is None:
331 msg = "Either rt or aquifer_pore_volumes must be provided"
332 raise ValueError(msg)
334 rt = residence_time(
335 flow=flow,
336 flow_tedges=flow_tedges,
337 aquifer_pore_volumes=aquifer_pore_volumes,
338 index=index,
339 direction=direction,
340 retardation_factor=retardation_factor,
341 )
343 expected_ndim = 2
344 if rt.ndim != expected_ndim:
345 msg = f"rt must be 2D with shape (n_pore_volumes, n_times), got {rt.ndim}D"
346 raise ValueError(msg)
348 n_aquifer_pore_volume = rt.shape[0]
349 return (n_aquifer_pore_volume - np.isnan(rt).sum(axis=0)) / n_aquifer_pore_volume
352def freundlich_retardation(
353 *,
354 concentration: npt.ArrayLike,
355 freundlich_k: float,
356 freundlich_n: float,
357 bulk_density: float,
358 porosity: float,
359) -> npt.NDArray[np.floating]:
360 """
361 Compute concentration-dependent retardation factors using Freundlich isotherm.
363 The Freundlich isotherm relates sorbed concentration S to aqueous concentration C:
364 S = rho_f * C^n
366 The retardation factor is computed as:
367 R = 1 + (rho_b/θ) * dS/dC = 1 + (rho_b/θ) * rho_f * n * C^(n-1)
369 Parameters
370 ----------
371 concentration : array-like
372 Concentration of compound in water [mass/volume].
373 Length should match flow (i.e., len(flow_tedges) - 1).
374 freundlich_k : float
375 Freundlich coefficient [(m³/kg)^(1/n)].
376 freundlich_n : float
377 Freundlich sorption exponent [dimensionless].
378 bulk_density : float
379 Bulk density of aquifer material [mass/volume].
380 porosity : float
381 Porosity of aquifer [dimensionless, 0-1].
383 Returns
384 -------
385 numpy.ndarray
386 Retardation factors for each flow interval.
387 Length equals len(concentration) for use as retardation_factor in residence_time.
389 Raises
390 ------
391 ValueError
392 If ``porosity`` is not in ``(0, 1)``, if ``bulk_density`` is not positive, if
393 ``freundlich_k`` is negative, or if any ``concentration`` is non-positive while
394 ``freundlich_n < 1`` (the retardation factor diverges as ``C -> 0``).
396 See Also
397 --------
398 residence_time : Compute residence times with retardation
399 gwtransport.advection.infiltration_to_extraction_front_tracking : Transport with nonlinear sorption
400 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and concentration-dependent retardation
402 Examples
403 --------
404 >>> concentration = np.array([0.1, 0.2, 0.3]) # same length as flow
405 >>> R = freundlich_retardation(
406 ... concentration=concentration,
407 ... freundlich_k=0.5,
408 ... freundlich_n=0.9,
409 ... bulk_density=1600, # kg/m3
410 ... porosity=0.35,
411 ... )
412 >>> # Use R in residence_time as retardation_factor
413 """
414 concentration = np.asarray(concentration)
416 # Validate physical parameters
417 if not 0 < porosity < 1:
418 msg = f"Porosity must be in (0, 1), got {porosity}"
419 raise ValueError(msg)
420 if bulk_density <= 0:
421 msg = f"Bulk density must be positive, got {bulk_density}"
422 raise ValueError(msg)
423 if freundlich_k < 0:
424 msg = f"Freundlich K must be non-negative, got {freundlich_k}"
425 raise ValueError(msg)
427 # For n < 1 the Freundlich retardation factor 1 + (rho_b/theta) * k_f * n * C^(n-1)
428 # diverges as C -> 0. Silently clamping concentration would produce a very large but
429 # finite value that depends on an arbitrary regularization constant; instead, refuse
430 # the call so the user can decide how to handle non-positive concentrations.
431 if freundlich_n < 1.0 and np.any(concentration <= 0):
432 msg = "concentration must be strictly positive when freundlich_n < 1 (retardation diverges as C -> 0)"
433 raise ValueError(msg)
435 return 1.0 + (bulk_density / porosity) * freundlich_k * freundlich_n * np.power(concentration, freundlich_n - 1)