Coverage for src/gwtransport/residence_time.py: 47%
62 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-22 18:30 +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:
9- residence_time: Compute the residence time of a retarded compound in the aquifer at indices.
10- residence_time_mean: Compute the mean residence time of a retarded compound in the aquifer
11 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",
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 : pandas.Series, 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 : str, optional
50 Direction of the flow. Either 'extraction' or 'infiltration'. Extraction refers to backward modeling: how many days ago did this extracted water infiltrate. Infiltration refers to forward modeling: how many days will it take for this infiltrated water to be extracted. Default is 'extraction'.
51 return_pandas_series : bool, optional
52 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.
54 Returns
55 -------
56 numpy.ndarray
57 Residence time of the retarded compound in the aquifer [days].
58 """
59 aquifer_pore_volume = np.atleast_1d(aquifer_pore_volume)
61 if flow_tedges is None:
62 msg = "flow_tedges must be provided"
63 raise ValueError(msg)
65 flow_tedges = pd.DatetimeIndex(flow_tedges)
66 if len(flow_tedges) != len(flow) + 1:
67 msg = "flow_tedges must have one more element than flow"
68 raise ValueError(msg)
70 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
71 flow_tdelta = np.diff(flow_tedges_days, prepend=0.0)
72 flow_values = np.concat(([0.0], np.asarray(flow)))
73 flow_cum = (flow_values * flow_tdelta).cumsum() # at flow_tedges and flow_tedges_days. First value is 0.
75 if index is None:
76 # If index is not provided return the residence time that matches with the index of the flow; at the center of the flow bin.
77 index_dates_days_extraction = (flow_tedges_days[:-1] + flow_tedges_days[1:]) / 2
78 flow_cum_at_index = (flow_cum[:-1] + flow_cum[1:]) / 2 # at the center of the flow bin
79 else:
80 index_dates_days_extraction = np.asarray((index - flow_tedges[0]) / np.timedelta64(1, "D"))
81 flow_cum_at_index = linear_interpolate(
82 flow_tedges_days, flow_cum, index_dates_days_extraction, left=np.nan, right=np.nan
83 )
85 if direction == "extraction":
86 # How many days ago was the extraced water infiltrated
87 a = flow_cum_at_index[None, :] - retardation_factor * aquifer_pore_volume[:, None]
88 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
89 data = index_dates_days_extraction - days
90 elif direction == "infiltration":
91 # In how many days the water that is infiltrated now be extracted
92 a = flow_cum_at_index[None, :] + retardation_factor * aquifer_pore_volume[:, None]
93 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
94 data = days - index_dates_days_extraction
95 else:
96 msg = "direction should be 'extraction' or 'infiltration'"
97 raise ValueError(msg)
99 if return_pandas_series:
100 warnings.warn(
101 "return_pandas_series parameter is deprecated and will be removed in a future version. "
102 "The function now returns numpy arrays by default.",
103 DeprecationWarning,
104 stacklevel=2,
105 )
106 if len(aquifer_pore_volume) > 1:
107 msg = "return_pandas_series=True is only supported for a single pore volume"
108 raise ValueError(msg)
109 return pd.Series(data=data[0], index=index, name=f"residence_time_{direction}")
110 return data
113def residence_time_mean(
114 flow, flow_tedges, tedges_out, aquifer_pore_volume, *, direction="extraction", retardation_factor=1.0
115):
116 """
117 Compute the mean residence time of a retarded compound in the aquifer between specified time edges.
119 This function calculates the average residence time of a retarded compound in the aquifer
120 between specified time intervals. It can compute both backward modeling (extraction direction:
121 when was extracted water infiltrated) and forward modeling (infiltration direction: when will
122 infiltrated water be extracted).
124 The function handles time series data by computing the cumulative flow and using linear
125 interpolation and averaging to determine mean residence times between the specified time edges.
127 Parameters
128 ----------
129 flow : array-like
130 Flow rate of water in the aquifer [m3/day]. Should be an array of flow values
131 corresponding to the intervals defined by flow_tedges.
132 flow_tedges : array-like
133 Time edges for the flow data, as datetime64 objects. These define the time
134 intervals for which the flow values are provided.
135 tedges_out : array-like
136 Output time edges as datetime64 objects. These define the intervals for which
137 the mean residence times will be calculated.
138 aquifer_pore_volume : float or array-like
139 Pore volume of the aquifer [m3]. Can be a single value or an array of values
140 for multiple pore volume scenarios.
141 direction : {'extraction', 'infiltration'}, optional
142 Direction of the flow calculation:
143 * 'extraction': Backward modeling - how many days ago was the extracted water infiltrated
144 * 'infiltration': Forward modeling - how many days until the infiltrated water is extracted
145 Default is 'extraction'.
146 retardation_factor : float, optional
147 Retardation factor of the compound in the aquifer [dimensionless].
148 A value greater than 1.0 indicates that the compound moves slower than water.
149 Default is 1.0 (no retardation).
151 Returns
152 -------
153 numpy.ndarray
154 Mean residence time of the retarded compound in the aquifer [days] for each interval
155 defined by tedges_out. The first dimension corresponds to the different pore volumes
156 and the second to the residence times between tedges_out.
158 Notes
159 -----
160 - The function converts datetime objects to days since the start of the time series.
161 - For extraction direction, the function computes how many days ago water was infiltrated.
162 - For infiltration direction, the function computes how many days until water will be extracted.
163 - The function uses linear interpolation for computing residence times at specific points
164 and linear averaging for computing mean values over intervals.
166 Examples
167 --------
168 >>> import pandas as pd
169 >>> import numpy as np
170 >>> # Create sample flow data
171 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
172 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day
173 >>> pore_volume = 200.0 # Aquifer pore volume in m³
174 >>> # Calculate mean residence times
175 >>> mean_times = residence_time_mean(
176 ... flow=flow_values,
177 ... flow_tedges=flow_dates,
178 ... tedges_out=flow_dates,
179 ... aquifer_pore_volume=pore_volume,
180 ... direction="extraction",
181 ... )
182 >>> # With constant flow of 100 m³/day and pore volume of 200 m³,
183 >>> # mean residence time should be approximately 2 days
184 >>> print(mean_times) # Output: [np.nan, np.nan, 2.0, 2.0, 2.0, ..., 2.0]
185 """
186 flow = np.asarray(flow)
187 flow_tedges = pd.DatetimeIndex(flow_tedges)
188 tedges_out = pd.DatetimeIndex(tedges_out)
189 aquifer_pore_volume = np.atleast_1d(aquifer_pore_volume)
191 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
192 tedges_out_days = np.asarray((tedges_out - flow_tedges[0]) / np.timedelta64(1, "D"))
194 # compute cumulative flow at flow_tedges and flow_tedges_days
195 flow_cum = np.diff(flow_tedges_days, prepend=0.0)
196 flow_cum[1:] *= flow
197 flow_cum = flow_cum.cumsum()
199 if direction == "extraction":
200 # How many days ago was the extraced water infiltrated
201 a = flow_cum[None, :] - retardation_factor * aquifer_pore_volume[:, None]
202 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
203 data_edges = flow_tedges_days - days
204 # Process each pore volume (row) separately. Although linear_average supports 2D x_edges,
205 # our use case is different: multiple time series (different y_data) with same edges,
206 # rather than same time series with multiple edge sets.
207 data_avg = np.array([linear_average(flow_tedges_days, y, tedges_out_days)[0] for y in data_edges])
208 elif direction == "infiltration":
209 # In how many days the water that is infiltrated now be extracted
210 a = flow_cum[None, :] + retardation_factor * aquifer_pore_volume[:, None]
211 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
212 data_edges = days - flow_tedges_days
213 # Process each pore volume (row) separately. Although linear_average supports 2D x_edges,
214 # our use case is different: multiple time series (different y_data) with same edges,
215 # rather than same time series with multiple edge sets.
216 data_avg = np.array([linear_average(flow_tedges_days, y, tedges_out_days)[0] for y in data_edges])
217 else:
218 msg = "direction should be 'extraction' or 'infiltration'"
219 raise ValueError(msg)
220 return data_avg