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