Coverage for src / gwtransport / residence_time.py: 91%
91 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-04 18:51 +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 See Also
79 --------
80 residence_time_mean : Compute mean residence time over time intervals
81 gwtransport.advection.gamma_infiltration_to_extraction : Use residence times for transport
82 gwtransport.logremoval.residence_time_to_log_removal : Convert residence time to log removal
83 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
84 :ref:`concept-retardation-factor` : Slower movement due to sorption
85 """
86 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes)
87 flow_tedges = pd.DatetimeIndex(flow_tedges)
89 # Convert to arrays for type safety
90 flow = np.asarray(flow)
91 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes)
93 if len(flow_tedges) != len(flow) + 1:
94 msg = "tedges must have one more element than flow"
95 raise ValueError(msg)
97 # Check for negative flow values - physically invalid
98 if np.any(flow < 0):
99 # Return NaN array with correct shape
100 n_output = len(flow_tedges) - 1 if index is None else len(index)
101 n_pore_volumes = len(aquifer_pore_volumes)
102 return np.full((n_pore_volumes, n_output), np.nan)
104 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
105 flow_tdelta = np.diff(flow_tedges_days)
106 flow_cum = np.concatenate((
107 [0.0],
108 (flow * flow_tdelta).cumsum(),
109 )) # at flow_tedges and flow_tedges_days. First value is 0.
111 if index is None:
112 # If index is not provided return the residence time that matches with the index of the flow; at the center of the flow bin.
113 index_dates_days_extraction = (flow_tedges_days[:-1] + flow_tedges_days[1:]) / 2
114 flow_cum_at_index = (flow_cum[:-1] + flow_cum[1:]) / 2 # at the center of the flow bin
115 else:
116 index_dates_days_extraction = np.asarray((index - flow_tedges[0]) / np.timedelta64(1, "D"))
117 flow_cum_at_index = linear_interpolate(
118 x_ref=flow_tedges_days, y_ref=flow_cum, x_query=index_dates_days_extraction, left=np.nan, right=np.nan
119 )
121 if direction == "extraction_to_infiltration":
122 # How many days ago was the extraced water infiltrated
123 a = flow_cum_at_index[None, :] - retardation_factor * aquifer_pore_volumes[:, None]
124 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
125 data = index_dates_days_extraction - days
126 elif direction == "infiltration_to_extraction":
127 # In how many days the water that is infiltrated now be extracted
128 a = flow_cum_at_index[None, :] + retardation_factor * aquifer_pore_volumes[:, None]
129 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
130 data = days - index_dates_days_extraction
131 else:
132 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
133 raise ValueError(msg)
135 return data
138def residence_time_mean(
139 *,
140 flow: npt.ArrayLike,
141 flow_tedges: pd.DatetimeIndex | np.ndarray,
142 tedges_out: pd.DatetimeIndex | np.ndarray,
143 aquifer_pore_volumes: npt.ArrayLike,
144 direction: str = "extraction_to_infiltration",
145 retardation_factor: float = 1.0,
146) -> npt.NDArray[np.floating]:
147 """
148 Compute the mean residence time of a retarded compound in the aquifer between specified time edges.
150 This function calculates the average residence time of a retarded compound in the aquifer
151 between specified time intervals. It can compute both extraction to infiltration modeling (extraction direction:
152 when was extracted water infiltrated) and infiltration to extraction modeling (infiltration direction: when will
153 infiltrated water be extracted).
155 The function handles time series data by computing the cumulative flow and using linear
156 interpolation and averaging to determine mean residence times between the specified time edges.
158 Parameters
159 ----------
160 flow : array-like
161 Flow rate of water in the aquifer [m3/day]. Should be an array of flow values
162 corresponding to the intervals defined by flow_tedges.
163 flow_tedges : array-like
164 Time edges for the flow data, as datetime64 objects. These define the time
165 intervals for which the flow values are provided.
166 tedges_out : array-like
167 Output time edges as datetime64 objects. These define the intervals for which
168 the mean residence times will be calculated.
169 aquifer_pore_volumes : float or array-like
170 Pore volume(s) of the aquifer [m3]. Can be a single value or an array of values
171 for multiple pore volume scenarios.
172 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
173 Direction of the flow calculation:
174 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
175 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
176 Default is 'extraction_to_infiltration'.
177 retardation_factor : float, optional
178 Retardation factor of the compound in the aquifer [dimensionless].
179 A value greater than 1.0 indicates that the compound moves slower than water.
180 Default is 1.0 (no retardation).
182 Returns
183 -------
184 numpy.ndarray
185 Mean residence time of the retarded compound in the aquifer [days] for each interval
186 defined by tedges_out. The first dimension corresponds to the different pore volumes
187 and the second to the residence times between tedges_out.
189 Notes
190 -----
191 - The function converts datetime objects to days since the start of the time series.
192 - For extraction_to_infiltration direction, the function computes how many days ago water was infiltrated.
193 - For infiltration_to_extraction direction, the function computes how many days until water will be extracted.
194 - The function uses linear interpolation for computing residence times at specific points
195 and linear averaging for computing mean values over intervals.
197 Examples
198 --------
199 >>> import pandas as pd
200 >>> import numpy as np
201 >>> from gwtransport.residence_time import residence_time_mean
202 >>> # Create sample flow data
203 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
204 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day
205 >>> pore_volume = 200.0 # Aquifer pore volume in m³
206 >>> # Calculate mean residence times
207 >>> mean_times = residence_time_mean(
208 ... flow=flow_values,
209 ... flow_tedges=flow_dates,
210 ... tedges_out=flow_dates,
211 ... aquifer_pore_volumes=pore_volume,
212 ... direction="extraction_to_infiltration",
213 ... )
214 >>> # With constant flow of 100 m³/day and pore volume of 200 m³,
215 >>> # mean residence time should be approximately 2 days
216 >>> print(mean_times) # doctest: +NORMALIZE_WHITESPACE
217 [[nan nan 2. 2. 2. 2. 2. 2. 2.]]
219 See Also
220 --------
221 residence_time : Compute residence time at specific time indices
222 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction
223 """
224 flow = np.asarray(flow)
225 flow_tedges = pd.DatetimeIndex(flow_tedges)
226 tedges_out = pd.DatetimeIndex(tedges_out)
227 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes)
229 # Check for negative flow values - physically invalid
230 if np.any(flow < 0):
231 # Return NaN array with correct shape
232 n_pore_volumes = len(aquifer_pore_volumes)
233 n_output_bins = len(tedges_out) - 1
234 return np.full((n_pore_volumes, n_output_bins), np.nan)
236 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
237 tedges_out_days = np.asarray((tedges_out - flow_tedges[0]) / np.timedelta64(1, "D"))
239 # compute cumulative flow at flow_tedges
240 flow_cum = np.concatenate(([0.0], np.cumsum(flow * np.diff(flow_tedges_days))))
242 if direction == "extraction_to_infiltration":
243 # How many days ago was the extraced water infiltrated
244 a = flow_cum[None, :] - retardation_factor * aquifer_pore_volumes[:, None]
245 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
246 data_edges = flow_tedges_days - days
247 # Process each pore volume (row) separately. Although linear_average supports 2D x_edges,
248 # our use case is different: multiple time series (different y_data) with same edges,
249 # rather than same time series with multiple edge sets.
250 data_avg = np.array([
251 linear_average(x_data=flow_tedges_days, y_data=y, x_edges=tedges_out_days)[0] for y in data_edges
252 ])
253 elif direction == "infiltration_to_extraction":
254 # In how many days the water that is infiltrated now be extracted
255 a = flow_cum[None, :] + retardation_factor * aquifer_pore_volumes[:, None]
256 days = linear_interpolate(x_ref=flow_cum, y_ref=flow_tedges_days, x_query=a, left=np.nan, right=np.nan)
257 data_edges = days - flow_tedges_days
258 # Process each pore volume (row) separately. Although linear_average supports 2D x_edges,
259 # our use case is different: multiple time series (different y_data) with same edges,
260 # rather than same time series with multiple edge sets.
261 data_avg = np.array([
262 linear_average(x_data=flow_tedges_days, y_data=y, x_edges=tedges_out_days)[0] for y in data_edges
263 ])
264 else:
265 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
266 raise ValueError(msg)
267 return data_avg
270def fraction_explained(
271 *,
272 rt: npt.NDArray[np.floating] | None = None,
273 flow: npt.ArrayLike | None = None,
274 flow_tedges: pd.DatetimeIndex | np.ndarray | None = None,
275 aquifer_pore_volumes: npt.ArrayLike | None = None,
276 index: pd.DatetimeIndex | np.ndarray | None = None,
277 direction: str = "extraction_to_infiltration",
278 retardation_factor: float = 1.0,
279) -> npt.NDArray[np.floating]:
280 """
281 Compute the fraction of the aquifer that is informed with respect to the retarded flow.
283 Parameters
284 ----------
285 rt : numpy.ndarray, optional
286 Pre-computed residence time array [days]. If not provided, it will be computed.
287 flow : array-like, optional
288 Flow rate of water in the aquifer [m3/day]. The length of `flow` should match the length of `flow_tedges` minus one.
289 flow_tedges : pandas.DatetimeIndex, optional
290 Time edges for the flow data. Used to compute the cumulative flow.
291 Has a length of one more than `flow`. Inbetween neighboring time edges, the flow is assumed constant.
292 aquifer_pore_volumes : float or array-like of float, optional
293 Pore volume(s) of the aquifer [m3]. Can be a single value or an array
294 of pore volumes representing different flow paths.
295 index : pandas.DatetimeIndex, optional
296 Index at which to compute the fraction. If left to None, the index of `flow` is used.
297 Default is None.
298 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
299 Direction of the flow calculation:
300 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
301 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
302 Default is 'extraction_to_infiltration'.
303 retardation_factor : float, optional
304 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
306 Returns
307 -------
308 numpy.ndarray
309 Fraction of the aquifer that is informed with respect to the retarded flow.
310 """
311 if rt is None:
312 # Validate that required parameters are provided for computing rt
313 if flow is None:
314 msg = "Either rt or flow must be provided"
315 raise ValueError(msg)
316 if flow_tedges is None:
317 msg = "Either rt or flow_tedges must be provided"
318 raise ValueError(msg)
319 if aquifer_pore_volumes is None:
320 msg = "Either rt or aquifer_pore_volumes must be provided"
321 raise ValueError(msg)
323 rt = residence_time(
324 flow=flow,
325 flow_tedges=flow_tedges,
326 aquifer_pore_volumes=aquifer_pore_volumes,
327 index=index,
328 direction=direction,
329 retardation_factor=retardation_factor,
330 )
332 _expected_ndim = 2
333 if rt.ndim != _expected_ndim:
334 msg = f"rt must be 2D with shape (n_pore_volumes, n_times), got {rt.ndim}D"
335 raise ValueError(msg)
337 n_aquifer_pore_volume = rt.shape[0]
338 return (n_aquifer_pore_volume - np.isnan(rt).sum(axis=0)) / n_aquifer_pore_volume
341def freundlich_retardation(
342 *,
343 concentration: npt.ArrayLike,
344 freundlich_k: float,
345 freundlich_n: float,
346 bulk_density: float,
347 porosity: float,
348) -> npt.NDArray[np.floating]:
349 """
350 Compute concentration-dependent retardation factors using Freundlich isotherm.
352 The Freundlich isotherm relates sorbed concentration S to aqueous concentration C:
353 S = rho_f * C^n
355 The retardation factor is computed as:
356 R = 1 + (rho_b/θ) * dS/dC = 1 + (rho_b/θ) * rho_f * n * C^(n-1)
358 Parameters
359 ----------
360 concentration : array-like
361 Concentration of compound in water [mass/volume].
362 Length should match flow (i.e., len(flow_tedges) - 1).
363 freundlich_k : float
364 Freundlich coefficient [(m³/kg)^(1/n)].
365 freundlich_n : float
366 Freundlich sorption exponent [dimensionless].
367 bulk_density : float
368 Bulk density of aquifer material [mass/volume].
369 porosity : float
370 Porosity of aquifer [dimensionless, 0-1].
372 Returns
373 -------
374 numpy.ndarray
375 Retardation factors for each flow interval.
376 Length equals len(concentration) for use as retardation_factor in residence_time.
378 Examples
379 --------
380 >>> concentration = np.array([0.1, 0.2, 0.3]) # same length as flow
381 >>> R = freundlich_retardation(
382 ... concentration=concentration,
383 ... freundlich_k=0.5,
384 ... freundlich_n=0.9,
385 ... bulk_density=1600, # kg/m3
386 ... porosity=0.35,
387 ... )
388 >>> # Use R in residence_time as retardation_factor
390 See Also
391 --------
392 residence_time : Compute residence times with retardation
393 gwtransport.advection.infiltration_to_extraction_front_tracking : Transport with nonlinear sorption
394 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and concentration-dependent retardation
395 """
396 concentration = np.asarray(concentration)
398 # Validate physical parameters
399 if not 0 < porosity < 1:
400 msg = f"Porosity must be in (0, 1), got {porosity}"
401 raise ValueError(msg)
402 if bulk_density <= 0:
403 msg = f"Bulk density must be positive, got {bulk_density}"
404 raise ValueError(msg)
405 if freundlich_k < 0:
406 msg = f"Freundlich K must be non-negative, got {freundlich_k}"
407 raise ValueError(msg)
409 concentration_safe = np.maximum(concentration, 1e-12) # Avoid zero concentration issues
410 return 1.0 + (bulk_density / porosity) * freundlich_k * freundlich_n * np.power(
411 concentration_safe, freundlich_n - 1
412 )