Coverage for src/gwtransport/residence_time.py: 100%
54 statements
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-06 15: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:
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 numpy as np
15import pandas as pd
17from gwtransport.utils import compute_time_edges, linear_average, linear_interpolate
20def residence_time(
21 flow=None,
22 flow_tedges=None,
23 flow_tstart=None,
24 flow_tend=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 timestamps of the flow data should be aligned with the time edges provided in `flow_tedges`.
39 If left to None, the function will raise an error. The length of `flow` should match the length of `flow_tedges` minus one.
40 flow_tedges : pandas.DatetimeIndex, optional
41 Time edges for the flow data. If provided, it is used to compute the cumulative flow.
42 If left to None, the index of `flow` is used. Has a length of one more than `flow`. Default is None.
43 flow_tstart : pandas.DatetimeIndex, optional
44 Timestamps aligned to the start of the flow measurement intervals. Preferably use flow_tedges,
45 but if not available this approach can be used for convenience. Has the same length as `flow`.
46 flow_tend : pandas.DatetimeIndex, optional
47 Timestamps aligned to the end of the flow measurement intervals. Preferably use flow_tedges,
48 but if not available this approach can be used for convenience. Has the same length as `flow`.
49 aquifer_pore_volume : float or array-like of float
50 Pore volume of the aquifer [m3].
51 index : pandas.DatetimeIndex, optional
52 Index at which to compute the residence time. If left to None, the index of `flow` is used.
53 If Default is None.
54 retardation_factor : float
55 Retardation factor of the compound in the aquifer [dimensionless].
56 direction : str, optional
57 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'.
58 return_pandas_series : bool, optional
59 If True, return a pandas Series with the residence time at the index provided. Only supported for a single aquifer pore volume.
61 Returns
62 -------
63 array
64 Residence time of the retarded compound in the aquifer [days].
65 """
66 aquifer_pore_volume = np.atleast_1d(aquifer_pore_volume)
68 flow_tedges = compute_time_edges(flow_tedges, flow_tstart, flow_tend, len(flow))
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 if len(aquifer_pore_volume) > 1:
101 msg = "return_pandas_series=True is only supported for a single pore volume"
102 raise ValueError(msg)
103 return pd.Series(data=data[0], index=index, name=f"residence_time_{direction}")
104 return data
107def residence_time_mean(
108 flow, flow_tedges, tedges_out, aquifer_pore_volume, *, direction="extraction", retardation_factor=1.0
109):
110 """
111 Compute the mean residence time of a retarded compound in the aquifer between specified time edges.
113 This function calculates the average residence time of a retarded compound in the aquifer
114 between specified time intervals. It can compute both backward modeling (extraction direction:
115 when was extracted water infiltrated) and forward modeling (infiltration direction: when will
116 infiltrated water be extracted).
118 The function handles time series data by computing the cumulative flow and using linear
119 interpolation and averaging to determine mean residence times between the specified time edges.
121 Parameters
122 ----------
123 flow : array-like
124 Flow rate of water in the aquifer [m3/day]. Should be an array of flow values
125 corresponding to the intervals defined by flow_tedges.
126 flow_tedges : array-like
127 Time edges for the flow data, as datetime64 objects. These define the time
128 intervals for which the flow values are provided.
129 tedges_out : array-like
130 Output time edges as datetime64 objects. These define the intervals for which
131 the mean residence times will be calculated.
132 aquifer_pore_volume : float or array-like
133 Pore volume of the aquifer [m3]. Can be a single value or an array of values
134 for multiple pore volume scenarios.
135 direction : {'extraction', 'infiltration'}, optional
136 Direction of the flow calculation:
137 * 'extraction': Backward modeling - how many days ago was the extracted water infiltrated
138 * 'infiltration': Forward modeling - how many days until the infiltrated water is extracted
139 Default is 'extraction'.
140 retardation_factor : float, optional
141 Retardation factor of the compound in the aquifer [dimensionless].
142 A value greater than 1.0 indicates that the compound moves slower than water.
143 Default is 1.0 (no retardation).
145 Returns
146 -------
147 numpy.ndarray
148 Mean residence time of the retarded compound in the aquifer [days] for each interval
149 defined by tedges_out. The first dimension corresponds to the different pore volumes
150 and the second to the residence times between tedges_out.
152 Notes
153 -----
154 - The function converts datetime objects to days since the start of the time series.
155 - For extraction direction, the function computes how many days ago water was infiltrated.
156 - For infiltration direction, the function computes how many days until water will be extracted.
157 - The function uses linear interpolation for computing residence times at specific points
158 and linear averaging for computing mean values over intervals.
160 Examples
161 --------
162 >>> import pandas as pd
163 >>> import numpy as np
164 >>> # Create sample flow data
165 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
166 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day
167 >>> pore_volume = 200.0 # Aquifer pore volume in m³
168 >>> # Calculate mean residence times
169 >>> mean_times = residence_time_mean(
170 ... flow=flow_values,
171 ... flow_tedges=flow_dates,
172 ... tedges_out=flow_dates,
173 ... aquifer_pore_volume=pore_volume,
174 ... direction="extraction",
175 ... )
176 >>> # With constant flow of 100 m³/day and pore volume of 200 m³,
177 >>> # mean residence time should be approximately 2 days
178 >>> print(mean_times) # Output: [np.nan, np.nan, 2.0, 2.0, 2.0, ..., 2.0]
179 """
180 flow = np.asarray(flow)
181 flow_tedges = pd.DatetimeIndex(flow_tedges)
182 tedges_out = pd.DatetimeIndex(tedges_out)
183 aquifer_pore_volume = np.atleast_1d(aquifer_pore_volume)
185 flow_tedges_days = np.asarray((flow_tedges - flow_tedges[0]) / np.timedelta64(1, "D"))
186 tedges_out_days = np.asarray((tedges_out - flow_tedges[0]) / np.timedelta64(1, "D"))
188 # compute cumulative flow at flow_tedges and flow_tedges_days
189 flow_cum = np.diff(flow_tedges_days, prepend=0.0)
190 flow_cum[1:] *= flow
191 flow_cum = flow_cum.cumsum()
193 if direction == "extraction":
194 # How many days ago was the extraced water infiltrated
195 a = flow_cum[None, :] - retardation_factor * aquifer_pore_volume[:, None]
196 days = linear_interpolate(flow_cum, flow_tedges_days, a, left=np.nan, right=np.nan)
197 data_edges = flow_tedges_days - days
198 data_avg = np.array([linear_average(flow_tedges_days, y, tedges_out_days) for y in data_edges])
199 elif direction == "infiltration":
200 # In how many days the water that is infiltrated now be extracted
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 = days - flow_tedges_days
204 data_avg = np.array([linear_average(flow_tedges_days, y, tedges_out_days) for y in data_edges])
205 else:
206 msg = "direction should be 'extraction' or 'infiltration'"
207 raise ValueError(msg)
208 return data_avg