Coverage for src/gwtransport/advection.py: 55%
51 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-22 18:30 +0000
1"""
2Advection Analysis for 1D Aquifer Systems.
4This module provides functions to analyze compound transport by advection
5in aquifer systems. It includes tools for computing concentrations of the extracted water
6based on the concentration of the infiltrating water, extraction data and aquifer properties.
8The model assumes requires the groundwaterflow to be reduced to a 1D system. On one side,
9water with a certain concentration infiltrates ('cin'), the water flows through the aquifer and
10the compound of interest flows through the aquifer with a retarded velocity. The water is
11extracted ('cout').
13Main functions:
14- forward: Compute the concentration of the extracted water by shifting cin with its residence time. This corresponds to a convolution operation.
15- gamma_forward: Similar to forward, but for a gamma distribution of aquifer pore volumes.
16- distribution_forward: Similar to forward, but for an arbitrairy distribution of aquifer pore volumes.
17"""
19import warnings
21import numpy as np
22import pandas as pd
24from gwtransport import gamma
25from gwtransport.residence_time import residence_time
26from gwtransport.utils import compute_time_edges, interp_series, linear_interpolate
29def forward(cin_series, flow_series, aquifer_pore_volume, retardation_factor=1.0, cout_index="cin"):
30 """
31 Compute the concentration of the extracted water by shifting cin with its residence time.
33 The compound is retarded in the aquifer with a retardation factor. The residence
34 time is computed based on the flow rate of the water in the aquifer and the pore volume
35 of the aquifer.
37 This function represents a forward operation (equivalent to convolution).
39 Parameters
40 ----------
41 cin_series : pandas.Series
42 Concentration of the compound in the extracted water [ng/m3]. The cin_series should be the average concentration of a time bin. The index should be a pandas.DatetimeIndex
43 and is labeled at the end of the time bin.
44 flow_series : pandas.Series
45 Flow rate of water in the aquifer [m3/day]. The flow_series should be the average flow of a time bin. The index should be a pandas.DatetimeIndex
46 and is labeled at the end of the time bin.
47 aquifer_pore_volume : float
48 Pore volume of the aquifer [m3].
49 cout_index : str, optional
50 The index of the output series. Can be 'cin', 'flow', or 'cout'. Default is 'cin'.
51 - 'cin': The output series will have the same index as `cin_series`.
52 - 'flow': The output series will have the same index as `flow_series`.
53 - 'cout': The output series will have the same index as `cin_series + residence_time`.
55 Returns
56 -------
57 numpy.ndarray
58 Concentration of the compound in the extracted water [ng/m3].
59 """
60 # Create flow tedges from the flow series index (assuming it's at the end of bins)
61 flow_tedges = compute_time_edges(tedges=None, tstart=None, tend=flow_series.index, number_of_bins=len(flow_series))
62 rt_series = residence_time(
63 flow=flow_series,
64 flow_tedges=flow_tedges,
65 index=cin_series.index,
66 aquifer_pore_volume=aquifer_pore_volume,
67 retardation_factor=retardation_factor,
68 direction="infiltration",
69 return_pandas_series=True,
70 )
72 rt = pd.to_timedelta(rt_series.values, unit="D", errors="coerce")
73 index = cin_series.index + rt
75 cout = pd.Series(data=cin_series.values, index=index, name="cout")
77 if cout_index == "cin":
78 return interp_series(cout, cin_series.index)
79 if cout_index == "flow":
80 # If cout_index is 'flow', we need to resample cout to the flow index
81 return interp_series(cout, flow_series.index)
82 if cout_index == "cout":
83 # If cout_index is 'cout', we return the cout as is
84 return cout.values
86 msg = f"Invalid cout_index: {cout_index}. Must be 'cin', 'flow', or 'cout'."
87 raise ValueError(msg)
90def backward(cout, flow, aquifer_pore_volume, retardation_factor=1.0, resample_dates=None):
91 """
92 Compute the concentration of the infiltrating water by shifting cout with its residence time.
94 This function represents a backward operation (equivalent to deconvolution).
96 Parameters
97 ----------
98 cout : pandas.Series
99 Concentration of the compound in the extracted water [ng/m3].
100 flow : pandas.Series
101 Flow rate of water in the aquifer [m3/day].
102 aquifer_pore_volume : float
103 Pore volume of the aquifer [m3].
105 Returns
106 -------
107 numpy.ndarray
108 Concentration of the compound in the infiltrating water [ng/m3].
109 """
110 msg = "Backward advection (deconvolution) is not implemented yet"
111 raise NotImplementedError(msg)
114def gamma_forward(
115 *,
116 cin,
117 cin_tedges,
118 cout_tedges,
119 flow,
120 flow_tedges,
121 alpha=None,
122 beta=None,
123 mean=None,
124 std=None,
125 n_bins=100,
126 retardation_factor=1.0,
127):
128 """
129 Compute the concentration of the extracted water by shifting cin with its residence time.
131 The compound is retarded in the aquifer with a retardation factor. The residence
132 time is computed based on the flow rate of the water in the aquifer and the pore volume
133 of the aquifer. The aquifer pore volume is approximated by a gamma distribution, with
134 parameters alpha and beta.
136 This function represents a forward operation (equivalent to convolution).
138 Provide either alpha and beta or mean and std.
140 Parameters
141 ----------
142 cin : pandas.Series
143 Concentration of the compound in infiltrating water or temperature of infiltrating
144 water.
145 cin_tedges : pandas.DatetimeIndex
146 Time edges for the concentration data. Used to compute the cumulative concentration.
147 Has a length of one more than `cin`.
148 cout_tedges : pandas.DatetimeIndex
149 Time edges for the output data. Used to compute the cumulative concentration.
150 Has a length of one more than `flow`.
151 flow : pandas.Series
152 Flow rate of water in the aquifer [m3/day].
153 flow_tedges : pandas.DatetimeIndex
154 Time edges for the flow data. Used to compute the cumulative flow.
155 Has a length of one more than `flow`.
156 alpha : float, optional
157 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0)
158 beta : float, optional
159 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0)
160 mean : float, optional
161 Mean of the gamma distribution.
162 std : float, optional
163 Standard deviation of the gamma distribution.
164 n_bins : int
165 Number of bins to discretize the gamma distribution.
166 retardation_factor : float
167 Retardation factor of the compound in the aquifer.
169 Returns
170 -------
171 numpy.ndarray
172 Concentration of the compound in the extracted water [ng/m3] or temperature.
173 """
174 bins = gamma.bins(alpha=alpha, beta=beta, mean=mean, std=std, n_bins=n_bins)
175 return distribution_forward(
176 cin=cin,
177 cin_tedges=cin_tedges,
178 cout_tedges=cout_tedges,
179 flow=flow,
180 flow_tedges=flow_tedges,
181 aquifer_pore_volumes=bins["expected_value"],
182 retardation_factor=retardation_factor,
183 )
186def gamma_backward(cout, flow, alpha, beta, n_bins=100, retardation_factor=1.0):
187 """
188 Compute the concentration of the infiltrating water by shifting cout with its residence time.
190 This function represents a backward operation (equivalent to deconvolution).
192 Parameters
193 ----------
194 cout : pandas.Series
195 Concentration of the compound in the extracted water [ng/m3].
196 flow : pandas.Series
197 Flow rate of water in the aquifer [m3/day].
198 alpha : float
199 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0)
200 beta : float
201 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0)
202 n_bins : int
203 Number of bins to discretize the gamma distribution.
204 retardation_factor : float
205 Retardation factor of the compound in the aquifer.
207 Returns
208 -------
209 NotImplementedError
210 This function is not yet implemented.
211 """
212 msg = "Backward advection gamma (deconvolution) is not implemented yet"
213 raise NotImplementedError(msg)
216def distribution_forward(
217 *,
218 cin,
219 cin_tedges,
220 cout_tedges,
221 flow,
222 flow_tedges,
223 aquifer_pore_volumes,
224 retardation_factor=1.0,
225):
226 """
227 Similar to forward_advection, but with a distribution of aquifer pore volumes.
229 Parameters
230 ----------
231 cin : pandas.Series
232 Concentration of the compound in infiltrating water or temperature of infiltrating
233 water.
234 cin_tedges : pandas.DatetimeIndex
235 Time edges for the concentration data. Used to compute the cumulative concentration.
236 Has a length of one more than `cin`.
237 cout_tedges : pandas.DatetimeIndex
238 Time edges for the output data. Used to compute the cumulative concentration.
239 Has a length of one more than `flow`.
240 flow : pandas.Series
241 Flow rate of water in the aquifer [m3/day].
242 flow_tedges : pandas.DatetimeIndex
243 Time edges for the flow data. Used to compute the cumulative flow.
244 Has a length of one more than `flow`.
245 aquifer_pore_volumes : array-like
246 Array of aquifer pore volumes [m3] representing the distribution.
247 retardation_factor : float
248 Retardation factor of the compound in the aquifer.
250 Returns
251 -------
252 numpy.ndarray
253 Concentration of the compound in the extracted water or temperature. Same units as cin.
254 """
255 cin_tedges = pd.DatetimeIndex(cin_tedges)
256 cout_tedges = pd.DatetimeIndex(cout_tedges)
257 if len(cin_tedges) != len(cin) + 1:
258 msg = "cin_tedges must have one more element than cin"
259 raise ValueError(msg)
260 if len(cout_tedges) != len(flow) + 1:
261 msg = "cout_tedges must have one more element than flow"
262 raise ValueError(msg)
264 cout_time = cout_tedges[:-1] + (cout_tedges[1:] - cout_tedges[:-1]) / 2
266 # Use residence time at cout_time for all pore volumes
267 rt_array = residence_time(
268 flow=flow,
269 flow_tedges=flow_tedges,
270 index=cout_time,
271 aquifer_pore_volume=aquifer_pore_volumes,
272 retardation_factor=retardation_factor,
273 direction="extraction",
274 ).astype("timedelta64[D]")
275 day_of_infiltration_array = cout_time.values[None] - rt_array
277 cin_sum = np.concat(([0.0], cin.cumsum())) # Add a zero at the beginning for cumulative sum
278 cin_sum_interpolated = linear_interpolate(cin_tedges, cin_sum, day_of_infiltration_array)
279 n_measurements = linear_interpolate(cin_tedges, np.arange(cin_tedges.size), day_of_infiltration_array)
280 cout_arr = np.diff(cin_sum_interpolated, axis=0) / np.diff(n_measurements, axis=0)
282 with warnings.catch_warnings():
283 warnings.filterwarnings(action="ignore", message="Mean of empty slice")
284 return np.nanmean(cout_arr, axis=0)
287def distribution_backward(cout, flow, aquifer_pore_volume_edges, retardation_factor=1.0):
288 """
289 Compute the concentration of the infiltrating water from the extracted water concentration considering a distribution of aquifer pore volumes.
291 This function represents a backward operation (equivalent to deconvolution).
293 Parameters
294 ----------
295 cout : pandas.Series
296 Concentration of the compound in the extracted water [ng/m3].
297 flow : pandas.Series
298 Flow rate of water in the aquifer [m3/day].
299 aquifer_pore_volume_edges : array-like
300 Edges of the bins that define the distribution of the aquifer pore volume.
301 Of size nbins + 1 [m3].
302 retardation_factor : float
303 Retardation factor of the compound in the aquifer.
305 Returns
306 -------
307 NotImplementedError
308 This function is not yet implemented.
309 """
310 msg = "Backward advection distribution (deconvolution) is not implemented yet"
311 raise NotImplementedError(msg)