Coverage for src/gwtransport/residence_time.py: 97%
71 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-09-29 20:01 +0000
1"""
2Residence time of a compound in the aquifer.
4This module provides functions to compute the residence time of a compound in the aquifer.
5The residence time is the time it takes for the compound to travel from the infiltration
6point to the extraction point. The compound is retarded in the aquifer with a retardation factor.
8Main functions:
10- residence_time: Compute the residence time of a retarded compound in the aquifer at indices.
11- residence_time_mean: Compute the mean residence time of a retarded compound in the aquifer between specified time edges.
12"""
14import warnings
16import numpy as np
17import numpy.typing as npt
18import pandas as pd
20from gwtransport.utils import linear_average, linear_interpolate
23def residence_time(
24 flow: npt.ArrayLike | None = None,
25 flow_tedges: pd.DatetimeIndex | np.ndarray | None = None,
26 aquifer_pore_volume: npt.ArrayLike | None = None,
27 index: pd.DatetimeIndex | np.ndarray | None = None,
28 retardation_factor: float = 1.0,
29 direction: str = "extraction_to_infiltration",
30 *,
31 return_pandas_series: bool = False,
32):
33 """
34 Compute the residence time of retarded compound in the water in the aquifer.
36 Parameters
37 ----------
38 flow : array-like
39 Flow rate of water in the aquifer [m3/day]. The length of `flow` should match the length of `flow_tedges` minus one.
40 flow_tedges : pandas.DatetimeIndex
41 Time edges for the flow data. Used to compute the cumulative flow.
42 Has a length of one more than `flow`.
43 aquifer_pore_volume : float or array-like of float
44 Pore volume of the aquifer [m3].
45 index : pandas.DatetimeIndex, optional
46 Index at which to compute the residence time. If left to None, the index of `flow` is used.
47 Default is None.
48 retardation_factor : float
49 Retardation factor of the compound in the aquifer [dimensionless].
50 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
51 Direction of the flow calculation:
52 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
53 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
54 Default is 'extraction_to_infiltration'.
55 return_pandas_series : bool, optional
56 If True, return a pandas Series with the residence time at the index provided. Only supported for a single aquifer pore volume. This parameter is deprecated and will be removed in a future version.
58 Returns
59 -------
60 numpy.ndarray
61 Residence time of the retarded compound in the aquifer [days].
62 """
63 aquifer_pore_volume = np.atleast_1d(aquifer_pore_volume)
65 if flow_tedges is None:
66 msg = "flow_tedges must be provided"
67 raise ValueError(msg)
69 if flow is None:
70 msg = "flow must be provided"
71 raise ValueError(msg)
73 flow_tedges = pd.DatetimeIndex(flow_tedges)
74 if len(flow_tedges) != len(flow) + 1:
75 msg = "tedges must have one more element than flow"
76 raise ValueError(msg)
78 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
79 flow_tdelta = np.diff(flow_tedges_days, prepend=0.0)
80 flow_values = np.concatenate(([0.0], np.asarray(flow)))
81 flow_cum = (flow_values * flow_tdelta).cumsum() # at flow_tedges and flow_tedges_days. First value is 0.
83 if index is None:
84 # If index is not provided return the residence time that matches with the index of the flow; at the center of the flow bin.
85 index_dates_days_extraction = (flow_tedges_days[:-1] + flow_tedges_days[1:]) / 2
86 flow_cum_at_index = (flow_cum[:-1] + flow_cum[1:]) / 2 # at the center of the flow bin
87 else:
88 index_dates_days_extraction = np.asarray((index - flow_tedges[0]) / np.timedelta64(1, "D"))
89 flow_cum_at_index = linear_interpolate(
90 flow_tedges_days, flow_cum, index_dates_days_extraction, left=np.nan, right=np.nan
91 )
93 if direction == "extraction_to_infiltration":
94 # How many days ago was the extraced water infiltrated
95 a = flow_cum_at_index[None, :] - retardation_factor * aquifer_pore_volume[:, None]
96 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
97 data = index_dates_days_extraction - days
98 elif direction == "infiltration_to_extraction":
99 # In how many days the water that is infiltrated now be extracted
100 a = flow_cum_at_index[None, :] + retardation_factor * aquifer_pore_volume[:, None]
101 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
102 data = days - index_dates_days_extraction
103 else:
104 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
105 raise ValueError(msg)
107 if return_pandas_series:
108 # If multiple pore volumes were provided, raise the explicit error first so that
109 # callers (and tests) see the ValueError before any deprecation warning. When
110 # running with warnings-as-errors, emitting the warning before the error would
111 # cause the test to fail on the warning instead of the intended ValueError.
112 if len(aquifer_pore_volume) > 1:
113 msg = "return_pandas_series=True is only supported for a single pore volume"
114 raise ValueError(msg)
115 warnings.warn(
116 "return_pandas_series parameter is deprecated and will be removed in a future version. "
117 "The function now returns numpy arrays by default.",
118 FutureWarning,
119 stacklevel=2,
120 )
121 return pd.Series(data=data[0], index=index, name=f"residence_time_{direction}")
122 return data
125def residence_time_mean(
126 flow: npt.ArrayLike,
127 flow_tedges: pd.DatetimeIndex | np.ndarray,
128 tedges_out: pd.DatetimeIndex | np.ndarray,
129 aquifer_pore_volume: npt.ArrayLike,
130 *,
131 direction: str = "extraction_to_infiltration",
132 retardation_factor: float = 1.0,
133):
134 """
135 Compute the mean residence time of a retarded compound in the aquifer between specified time edges.
137 This function calculates the average residence time of a retarded compound in the aquifer
138 between specified time intervals. It can compute both extraction to infiltration modeling (extraction direction:
139 when was extracted water infiltrated) and infiltration to extraction modeling (infiltration direction: when will
140 infiltrated water be extracted).
142 The function handles time series data by computing the cumulative flow and using linear
143 interpolation and averaging to determine mean residence times between the specified time edges.
145 Parameters
146 ----------
147 flow : array-like
148 Flow rate of water in the aquifer [m3/day]. Should be an array of flow values
149 corresponding to the intervals defined by flow_tedges.
150 flow_tedges : array-like
151 Time edges for the flow data, as datetime64 objects. These define the time
152 intervals for which the flow values are provided.
153 tedges_out : array-like
154 Output time edges as datetime64 objects. These define the intervals for which
155 the mean residence times will be calculated.
156 aquifer_pore_volume : float or array-like
157 Pore volume of the aquifer [m3]. Can be a single value or an array of values
158 for multiple pore volume scenarios.
159 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
160 Direction of the flow calculation:
161 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
162 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
163 Default is 'extraction_to_infiltration'.
164 retardation_factor : float, optional
165 Retardation factor of the compound in the aquifer [dimensionless].
166 A value greater than 1.0 indicates that the compound moves slower than water.
167 Default is 1.0 (no retardation).
169 Returns
170 -------
171 numpy.ndarray
172 Mean residence time of the retarded compound in the aquifer [days] for each interval
173 defined by tedges_out. The first dimension corresponds to the different pore volumes
174 and the second to the residence times between tedges_out.
176 Notes
177 -----
178 - The function converts datetime objects to days since the start of the time series.
179 - For extraction_to_infiltration direction, the function computes how many days ago water was infiltrated.
180 - For infiltration_to_extraction direction, the function computes how many days until water will be extracted.
181 - The function uses linear interpolation for computing residence times at specific points
182 and linear averaging for computing mean values over intervals.
184 Examples
185 --------
186 >>> import pandas as pd
187 >>> import numpy as np
188 >>> # Create sample flow data
189 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
190 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day
191 >>> pore_volume = 200.0 # Aquifer pore volume in m³
192 >>> # Calculate mean residence times
193 >>> mean_times = residence_time_mean(
194 ... flow=flow_values,
195 ... flow_tedges=flow_dates,
196 ... tedges_out=flow_dates,
197 ... aquifer_pore_volume=pore_volume,
198 ... direction="extraction_to_infiltration",
199 ... )
200 >>> # With constant flow of 100 m³/day and pore volume of 200 m³,
201 >>> # mean residence time should be approximately 2 days
202 >>> print(mean_times) # Output: [np.nan, np.nan, 2.0, 2.0, 2.0, ..., 2.0]
203 """
204 flow = np.asarray(flow)
205 flow_tedges = pd.DatetimeIndex(flow_tedges)
206 tedges_out = pd.DatetimeIndex(tedges_out)
207 aquifer_pore_volume = np.atleast_1d(aquifer_pore_volume)
209 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
210 tedges_out_days = np.asarray((tedges_out - flow_tedges[0]) / np.timedelta64(1, "D"))
212 # compute cumulative flow at flow_tedges and flow_tedges_days
213 flow_cum = np.diff(flow_tedges_days, prepend=0.0)
214 flow_cum[1:] *= flow
215 flow_cum = flow_cum.cumsum()
217 if direction == "extraction_to_infiltration":
218 # How many days ago was the extraced water infiltrated
219 a = flow_cum[None, :] - retardation_factor * aquifer_pore_volume[:, None]
220 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
221 data_edges = flow_tedges_days - days
222 # Process each pore volume (row) separately. Although linear_average supports 2D x_edges,
223 # our use case is different: multiple time series (different y_data) with same edges,
224 # rather than same time series with multiple edge sets.
225 data_avg = np.array([linear_average(flow_tedges_days, y, tedges_out_days)[0] for y in data_edges])
226 elif direction == "infiltration_to_extraction":
227 # In how many days the water that is infiltrated now be extracted
228 a = flow_cum[None, :] + retardation_factor * aquifer_pore_volume[:, None]
229 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
230 data_edges = days - flow_tedges_days
231 # Process each pore volume (row) separately. Although linear_average supports 2D x_edges,
232 # our use case is different: multiple time series (different y_data) with same edges,
233 # rather than same time series with multiple edge sets.
234 data_avg = np.array([linear_average(flow_tedges_days, y, tedges_out_days)[0] for y in data_edges])
235 else:
236 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'"
237 raise ValueError(msg)
238 return data_avg
241def fraction_explained(
242 rt=None,
243 flow=None,
244 flow_tedges=None,
245 aquifer_pore_volume=None,
246 index=None,
247 retardation_factor=1.0,
248 direction="extraction_to_infiltration",
249):
250 """
251 Compute the fraction of the aquifer that is informed with respect to the retarded flow.
253 Parameters
254 ----------
255 rt : numpy.ndarray, optional
256 Pre-computed residence time array [days]. If not provided, it will be computed.
257 flow : array-like, optional
258 Flow rate of water in the aquifer [m3/day]. The length of `flow` should match the length of `flow_tedges` minus one.
259 flow_tedges : pandas.DatetimeIndex, optional
260 Time edges for the flow data. Used to compute the cumulative flow.
261 Has a length of one more than `flow`. Inbetween neighboring time edges, the flow is assumed constant.
262 aquifer_pore_volume : float or array-like of float, optional
263 Pore volume of the aquifer [m3].
264 index : pandas.DatetimeIndex, optional
265 Index at which to compute the fraction. If left to None, the index of `flow` is used.
266 Default is None.
267 retardation_factor : float or array-like of float, optional
268 Retardation factor of the compound in the aquifer [dimensionless].
269 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional
270 Direction of the flow calculation:
271 * 'extraction_to_infiltration': Extraction to infiltration modeling - how many days ago was the extracted water infiltrated
272 * 'infiltration_to_extraction': Infiltration to extraction modeling - how many days until the infiltrated water is extracted
273 Default is 'extraction_to_infiltration'.
274 return_pandas_series : bool, optional
275 If True, return a pandas Series with the residence time at the index provided. Only supported for a single aquifer pore volume. This parameter is deprecated and will be removed in a future version.
277 Returns
278 -------
279 numpy.ndarray
280 Fraction of the aquifer that is informed with respect to the retarded flow.
281 """
282 if rt is None:
283 rt = residence_time(
284 flow=flow,
285 flow_tedges=flow_tedges,
286 aquifer_pore_volume=aquifer_pore_volume,
287 index=index,
288 retardation_factor=retardation_factor,
289 direction=direction,
290 )
292 n_aquifer_pore_volume = rt.shape[0]
293 return (n_aquifer_pore_volume - np.isnan(rt).sum(axis=0)) / n_aquifer_pore_volume