Coverage for src/gwtransport/advection.py: 83%

1""" 

2Advection Analysis for 1D Aquifer Systems. 

3 

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. 

7 

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'). 

12 

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""" 

18 

19import warnings 

20 

21import numpy as np 

22import pandas as pd 

23 

24from gwtransport import gamma 

25from gwtransport.residence_time import residence_time 

26from gwtransport.utils import compute_time_edges, interp_series, linear_interpolate 

27 

28 

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. 

32 

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. 

36 

37 This function represents a forward operation (equivalent to convolution). 

38 

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`. 

54 

55 Returns 

56 ------- 

57 pandas.Series 

58 Concentration of the compound in the extracted water [ng/m3]. 

59 """ 

60 rt_series = residence_time( 

61 flow=flow_series, 

62 flow_tedges=None, 

63 flow_tstart=None, 

64 flow_tend=flow_series.index, 

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 ) 

71 

72 rt = pd.to_timedelta(rt_series.values, unit="D", errors="coerce") 

73 index = cin_series.index + rt 

74 

75 cout = pd.Series(data=cin_series.values, index=index, name="cout") 

76 

77 if cout_index == "cin": 

78 return pd.Series(interp_series(cout, cin_series.index), index=cin_series.index, name="cout") 

79 if cout_index == "flow": 

80 # If cout_index is 'flow', we need to resample cout to the flow index 

81 return pd.Series(interp_series(cout, flow_series.index), index=flow_series.index, name="cout") 

82 if cout_index == "cout": 

83 # If cout_index is 'cout', we return the cout as is 

84 return cout 

85 

86 msg = f"Invalid cout_index: {cout_index}. Must be 'cin', 'flow', or 'cout'." 

87 raise ValueError(msg) 

88 

89 

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. 

93 

94 This function represents a backward operation (equivalent to deconvolution). 

95 

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]. 

104 

105 Returns 

106 ------- 

107 pandas.Series 

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) 

112 

113 

114def gamma_forward( 

115 *, 

116 cin, 

117 cin_tedges=None, 

118 cin_tstart=None, 

119 cin_tend=None, 

120 cout_tedges=None, 

121 cout_tstart=None, 

122 cout_tend=None, 

123 flow, 

124 flow_tedges=None, 

125 flow_tstart=None, 

126 flow_tend=None, 

127 alpha=None, 

128 beta=None, 

129 mean=None, 

130 std=None, 

131 n_bins=100, 

132 retardation_factor=1.0, 

133): 

134 """ 

135 Compute the concentration of the extracted water by shifting cin with its residence time. 

136 

137 The compound is retarded in the aquifer with a retardation factor. The residence 

138 time is computed based on the flow rate of the water in the aquifer and the pore volume 

139 of the aquifer. The aquifer pore volume is approximated by a gamma distribution, with 

140 parameters alpha and beta. 

141 

142 This function represents a forward operation (equivalent to convolution). 

143 

144 Provide: 

145 - either alpha and beta or mean and std. 

146 - either cin_tedges or cin_tstart or cin_tend. 

147 - either cout_tedges or cout_tstart or cout_tend. 

148 - either flow_tedges or flow_tstart or flow_tend. 

149 

150 Parameters 

151 ---------- 

152 cin : pandas.Series 

153 Concentration of the compound in infiltrating water or temperature of infiltrating 

154 water. 

155 cin_tedges : pandas.DatetimeIndex, optional 

156 Time edges for the concentration data. If provided, it is used to compute the cumulative concentration. 

157 cin_tstart : pandas.DatetimeIndex, optional 

158 Timestamps aligned to the start of the concentration measurement intervals. Preferably use cin_tedges, 

159 but if not available this approach can be used for convenience. 

160 cin_tend : pandas.DatetimeIndex, optional 

161 Timestamps aligned to the end of the concentration measurement intervals. Preferably use cin_tedges, 

162 but if not available this approach can be used for convenience. 

163 cout_tedges : pandas.DatetimeIndex, optional 

164 Time edges for the output data. If provided, it is used to compute the cumulative concentration. 

165 cout_tstart : pandas.DatetimeIndex, optional 

166 Timestamps aligned to the start of the output measurement intervals. Preferably use cout_tedges, 

167 but if not available this approach can be used for convenience. 

168 cout_tend : pandas.DatetimeIndex, optional 

169 Timestamps aligned to the end of the output measurement intervals. Preferably use cout_tedges, 

170 but if not available this approach can be used for convenience. 

171 flow : pandas.Series 

172 Flow rate of water in the aquifer [m3/day]. 

173 flow_tedges : pandas.DatetimeIndex, optional 

174 Time edges for the flow data. If provided, it is used to compute the cumulative flow. 

175 If left to None, the index of `flow` is used. Has a length of one more than `flow`. Default is None. 

176 flow_tstart : pandas.DatetimeIndex, optional 

177 Timestamps aligned to the start of the flow measurement intervals. Preferably use flow_tedges, 

178 but if not available this approach can be used for convenience. Has the same length as `flow`. 

179 flow_tend : pandas.DatetimeIndex, optional 

180 Timestamps aligned to the end of the flow measurement intervals. Preferably use flow_tedges, 

181 but if not available this approach can be used for convenience. Has the same length as `flow`. 

182 alpha : float, optional 

183 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0) 

184 beta : float, optional 

185 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0) 

186 mean : float, optional 

187 Mean of the gamma distribution. 

188 std : float, optional 

189 Standard deviation of the gamma distribution. 

190 n_bins : int 

191 Number of bins to discretize the gamma distribution. 

192 retardation_factor : float 

193 Retardation factor of the compound in the aquifer. 

194 

195 Returns 

196 ------- 

197 pandas.Series 

198 Concentration of the compound in the extracted water [ng/m3] or temperature. 

199 """ 

200 bins = gamma.bins(alpha=alpha, beta=beta, mean=mean, std=std, n_bins=n_bins) 

201 return distribution_forward( 

202 cin=cin, 

203 cin_tedges=cin_tedges, 

204 cin_tstart=cin_tstart, 

205 cin_tend=cin_tend, 

206 cout_tedges=cout_tedges, 

207 cout_tstart=cout_tstart, 

208 cout_tend=cout_tend, 

209 flow=flow, 

210 flow_tedges=flow_tedges, 

211 flow_tstart=flow_tstart, 

212 flow_tend=flow_tend, 

213 aquifer_pore_volume_edges=bins["edges"], 

214 retardation_factor=retardation_factor, 

215 ) 

216 

217 

218def gamma_backward(cout, flow, alpha, beta, n_bins=100, retardation_factor=1.0): 

219 """ 

220 Compute the concentration of the infiltrating water by shifting cout with its residence time. 

221 

222 This function represents a backward operation (equivalent to deconvolution). 

223 

224 Parameters 

225 ---------- 

226 cout : pandas.Series 

227 Concentration of the compound in the extracted water [ng/m3]. 

228 flow : pandas.Series 

229 Flow rate of water in the aquifer [m3/day]. 

230 alpha : float 

231 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0) 

232 beta : float 

233 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0) 

234 n_bins : int 

235 Number of bins to discretize the gamma distribution. 

236 retardation_factor : float 

237 Retardation factor of the compound in the aquifer. 

238 

239 Returns 

240 ------- 

241 NotImplementedError 

242 This function is not yet implemented. 

243 """ 

244 msg = "Backward advection gamma (deconvolution) is not implemented yet" 

245 raise NotImplementedError(msg) 

246 

247 

248def distribution_forward( 

249 *, 

250 cin, 

251 cin_tedges=None, 

252 cin_tstart=None, 

253 cin_tend=None, 

254 cout_tedges=None, 

255 cout_tstart=None, 

256 cout_tend=None, 

257 flow, 

258 flow_tedges=None, 

259 flow_tstart=None, 

260 flow_tend=None, 

261 aquifer_pore_volume_edges=None, 

262 retardation_factor=1.0, 

263): 

264 """ 

265 Similar to forward_advection, but with a distribution of aquifer pore volumes. 

266 

267 Provide: 

268 - either cin_tedges or cin_tstart or cin_tend. 

269 - either cout_tedges or cout_tstart or cout_tend. 

270 - either flow_tedges or flow_tstart or flow_tend. 

271 

272 Parameters 

273 ---------- 

274 cin : pandas.Series 

275 Concentration of the compound in infiltrating water or temperature of infiltrating 

276 water. 

277 cin_tedges : pandas.DatetimeIndex, optional 

278 Time edges for the concentration data. If provided, it is used to compute the cumulative concentration. 

279 cin_tstart : pandas.DatetimeIndex, optional 

280 Timestamps aligned to the start of the concentration measurement intervals. Preferably use cin_tedges, 

281 but if not available this approach can be used for convenience. 

282 cin_tend : pandas.DatetimeIndex, optional 

283 Timestamps aligned to the end of the concentration measurement intervals. Preferably use cin_tedges, 

284 but if not available this approach can be used for convenience. 

285 cout_tedges : pandas.DatetimeIndex, optional 

286 Time edges for the output data. If provided, it is used to compute the cumulative concentration. 

287 cout_tstart : pandas.DatetimeIndex, optional 

288 Timestamps aligned to the start of the output measurement intervals. Preferably use cout_tedges, 

289 but if not available this approach can be used for convenience. 

290 cout_tend : pandas.DatetimeIndex, optional 

291 Timestamps aligned to the end of the output measurement intervals. Preferably use cout_tedges, 

292 but if not available this approach can be used for convenience. 

293 flow : pandas.Series 

294 Flow rate of water in the aquifer [m3/day]. 

295 flow_tedges : pandas.DatetimeIndex, optional 

296 Time edges for the flow data. If provided, it is used to compute the cumulative flow. 

297 If left to None, the index of `flow` is used. Has a length of one more than `flow`. Default is None. 

298 flow_tstart : pandas.DatetimeIndex, optional 

299 Timestamps aligned to the start of the flow measurement intervals. Preferably use flow_tedges, 

300 but if not available this approach can be used for convenience. Has the same length as `flow`. 

301 flow_tend : pandas.DatetimeIndex, optional 

302 Timestamps aligned to the end of the flow measurement intervals. Preferably use flow_tedges, 

303 but if not available this approach can be used for convenience. Has the same length as `flow`. 

304 aquifer_pore_volume_edges : array-like 

305 Edges of the bins that define the distribution of the aquifer pore volume. 

306 Of size nbins + 1 [m3]. 

307 retardation_factor : float 

308 Retardation factor of the compound in the aquifer. 

309 

310 Returns 

311 ------- 

312 pandas.Series 

313 Concentration of the compound in the extracted water or temperature. Same units as cin. 

314 """ 

315 cin_tedges = compute_time_edges(tedges=cin_tedges, tstart=cin_tstart, tend=cin_tend, number_of_bins=len(cin)) 

316 cout_tedges = compute_time_edges(tedges=cout_tedges, tstart=cout_tstart, tend=cout_tend, number_of_bins=len(flow)) 

317 

318 if cout_tstart is not None: 

319 cout_time = cout_tstart 

320 elif cout_tend is not None: 

321 cout_time = cout_tend 

322 elif cout_tedges is not None: 

323 cout_time = cout_tedges[:-1] + (cout_tedges[1:] - cout_tedges[:-1]) / 2 

324 else: 

325 msg = "Either cout_tedges, cout_tstart, or cout_tend must be provided." 

326 raise ValueError(msg) 

327 

328 # Use residence time at cout_tedges 

329 rt_edges = residence_time( 

330 flow=flow, 

331 flow_tedges=flow_tedges, 

332 flow_tstart=flow_tstart, 

333 flow_tend=flow_tend, 

334 index=cout_time, 

335 aquifer_pore_volume=aquifer_pore_volume_edges, 

336 retardation_factor=retardation_factor, 

337 direction="extraction", 

338 ).astype("timedelta64[D]") 

339 day_of_infiltration_edges = cout_time.values[None] - rt_edges 

340 

341 cin_sum = np.concat(([0.0], cin.cumsum())) # Add a zero at the beginning for cumulative sum 

342 cin_sum_edges = linear_interpolate(cin_tedges, cin_sum, day_of_infiltration_edges) 

343 n_measurements = linear_interpolate(cin_tedges, np.arange(cin_tedges.size), day_of_infiltration_edges) 

344 cout_arr = np.diff(cin_sum_edges, axis=0) / np.diff(n_measurements, axis=0) 

345 

346 with warnings.catch_warnings(): 

347 warnings.filterwarnings(action="ignore", message="Mean of empty slice") 

348 cout_data = np.nanmean(cout_arr, axis=0) 

349 

350 return pd.Series(data=cout_data, index=cout_time, name="cout") 

351 

352 

353def distribution_backward(cout, flow, aquifer_pore_volume_edges, retardation_factor=1.0): 

354 """ 

355 Compute the concentration of the infiltrating water from the extracted water concentration considering a distribution of aquifer pore volumes. 

356 

357 This function represents a backward operation (equivalent to deconvolution). 

358 

359 Parameters 

360 ---------- 

361 cout : pandas.Series 

362 Concentration of the compound in the extracted water [ng/m3]. 

363 flow : pandas.Series 

364 Flow rate of water in the aquifer [m3/day]. 

365 aquifer_pore_volume_edges : array-like 

366 Edges of the bins that define the distribution of the aquifer pore volume. 

367 Of size nbins + 1 [m3]. 

368 retardation_factor : float 

369 Retardation factor of the compound in the aquifer. 

370 

371 Returns 

372 ------- 

373 NotImplementedError 

374 This function is not yet implemented. 

375 """ 

376 msg = "Backward advection distribution (deconvolution) is not implemented yet" 

377 raise NotImplementedError(msg)