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

138 statements  

« prev     ^ index     » next       coverage.py v7.13.4, created at 2026-03-04 18:51 +0000

1""" 

2Advective Transport Modeling for 1D Aquifer Systems. 

3 

4This module provides functions to model compound transport by advection in one-dimensional 

5aquifer systems, enabling prediction of solute or temperature concentrations in extracted 

6water based on infiltration data and aquifer properties. The model assumes one-dimensional 

7groundwater flow where water infiltrates with concentration ``cin``, flows through the aquifer 

8with pore volume distribution, compounds are transported with retarded velocity (retardation 

9factor >= 1.0), and water is extracted with concentration ``cout``. 

10 

11Available functions: 

12 

13- :func:`infiltration_to_extraction_series` - Single pore volume, time-shift only. Shifts 

14 infiltration time edges forward by residence time. Concentration values remain unchanged 

15 (cout = cin). No support for custom output time edges. Use case: Deterministic transport 

16 with single flow path. 

17 

18- :func:`infiltration_to_extraction` - Arbitrary pore volume distribution, convolution. 

19 Supports explicit distribution of aquifer pore volumes with flow-weighted averaging. 

20 Flexible output time resolution via cout_tedges. Use case: Known pore volume distribution 

21 from streamline analysis. 

22 

23- :func:`gamma_infiltration_to_extraction` - Gamma-distributed pore volumes, convolution. 

24 Models aquifer heterogeneity with 2-parameter gamma distribution. Parameterizable via 

25 (alpha, beta) or (mean, std). Discretizes gamma distribution into equal-probability bins. 

26 Use case: Heterogeneous aquifer with calibrated gamma parameters. 

27 

28Note on dispersion: The spreading from the pore volume distribution (APVD) represents 

29macrodispersion—aquifer-scale velocity heterogeneity that depends on both aquifer 

30properties and hydrological boundary conditions. If APVD was calibrated from 

31measurements, this spreading already includes microdispersion and molecular diffusion. 

32To add microdispersion and molecular diffusion separately (when APVD comes from 

33streamline analysis), use :mod:`gwtransport.diffusion`. 

34See :ref:`concept-dispersion-scales` for details. 

35 

36Note on cross-compound calibration: When APVD is calibrated from measurements of one 

37compound (e.g., temperature with D_m ~ 0.1 m²/day) and used to predict another (e.g., a 

38solute with D_m ~ 1e-4 m²/day), the molecular diffusion contribution baked into the 

39calibrated std may need correction. See :doc:`/examples/05_Diffusion_Dispersion` for 

40the correction procedure. 

41 

42- :func:`extraction_to_infiltration_series` - Single pore volume, time-shift only 

43 (deconvolution). Shifts extraction time edges backward by residence time. Concentration 

44 values remain unchanged (cin = cout). Symmetric inverse of infiltration_to_extraction_series. 

45 Use case: Backward tracing with single flow path. 

46 

47- :func:`extraction_to_infiltration` - Arbitrary pore volume distribution, deconvolution. 

48 Inverts forward transport for arbitrary pore volume distributions. Symmetric inverse of 

49 infiltration_to_extraction. Flow-weighted averaging in reverse direction. Use case: 

50 Estimating infiltration history from extraction data. 

51 

52- :func:`gamma_extraction_to_infiltration` - Gamma-distributed pore volumes, deconvolution. 

53 Inverts forward transport for gamma-distributed pore volumes. Symmetric inverse of 

54 gamma_infiltration_to_extraction. Use case: Calibrating infiltration conditions from 

55 extraction measurements. 

56 

57- :func:`infiltration_to_extraction_front_tracking` - Exact front tracking with Freundlich sorption. 

58 Event-driven algorithm that solves 1D advective transport with Freundlich isotherm using 

59 analytical integration of shock and rarefaction waves. Machine-precision physics (no numerical 

60 dispersion). Returns bin-averaged concentrations. Use case: Sharp concentration fronts with 

61 exact mass balance required, single deterministic flow path. 

62 

63- :func:`infiltration_to_extraction_front_tracking_detailed` - Front tracking with piecewise structure. 

64 Same as infiltration_to_extraction_front_tracking but also returns complete piecewise analytical 

65 structure including all events, segments, and callable analytical forms C(t). Use case: Detailed 

66 analysis of shock and rarefaction wave dynamics. 

67 

68This file is part of gwtransport which is released under AGPL-3.0 license. 

69See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details. 

70""" 

71 

72import numpy as np 

73import numpy.typing as npt 

74import pandas as pd 

75 

76from gwtransport import gamma 

77from gwtransport.advection_utils import ( 

78 _extraction_to_infiltration_weights, 

79 _infiltration_to_extraction_weights, 

80) 

81from gwtransport.fronttracking.math import EPSILON_FREUNDLICH_N, ConstantRetardation, FreundlichSorption 

82from gwtransport.fronttracking.output import compute_bin_averaged_concentration_exact 

83from gwtransport.fronttracking.solver import FrontTracker 

84from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave 

85from gwtransport.residence_time import residence_time 

86 

87 

88def infiltration_to_extraction_series( 

89 *, 

90 flow: npt.ArrayLike, 

91 tedges: pd.DatetimeIndex, 

92 aquifer_pore_volume: float, 

93 retardation_factor: float = 1.0, 

94) -> pd.DatetimeIndex: 

95 """ 

96 Compute extraction time edges from infiltration time edges using residence time shifts. 

97 

98 This function shifts infiltration time edges forward in time based on residence 

99 times computed from flow rates and aquifer properties. The concentration values remain 

100 unchanged (cout equals cin), only the time edges are shifted. This assumes a single pore 

101 volume (no distribution) and deterministic advective transport. 

102 

103 NOTE: This function is specifically designed for single aquifer pore volumes and does not 

104 support custom output time edges (cout_tedges). For distributions of aquifer pore volumes 

105 or custom output time grids, use `infiltration_to_extraction` instead. 

106 

107 Parameters 

108 ---------- 

109 flow : array-like 

110 Flow rate values in the aquifer [m3/day]. 

111 Length must match the number of time bins defined by tedges (len(tedges) - 1). 

112 tedges : pandas.DatetimeIndex 

113 Time edges defining bins for both cin and flow data. Has length of len(flow) + 1. 

114 aquifer_pore_volume : float 

115 Single aquifer pore volume [m3] used to compute residence times. 

116 retardation_factor : float, optional 

117 Retardation factor of the compound in the aquifer (default 1.0). 

118 Values > 1.0 indicate slower transport due to sorption/interaction. 

119 

120 Returns 

121 ------- 

122 pandas.DatetimeIndex 

123 Time edges for the extracted water concentration. Same length as tedges. 

124 The concentration values in the extracted water (cout) equal cin, but are 

125 aligned with these shifted time edges. 

126 

127 Examples 

128 -------- 

129 Basic usage with constant flow: 

130 

131 >>> import pandas as pd 

132 >>> import numpy as np 

133 >>> from gwtransport.utils import compute_time_edges 

134 >>> from gwtransport.advection import infiltration_to_extraction_series 

135 >>> 

136 >>> # Create input data 

137 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-10", freq="D") 

138 >>> tedges = compute_time_edges( 

139 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates) 

140 ... ) 

141 >>> 

142 >>> # Constant concentration and flow 

143 >>> cin = np.ones(len(dates)) * 10.0 

144 >>> flow = np.ones(len(dates)) * 100.0 # 100 m3/day 

145 >>> 

146 >>> # Run infiltration_to_extraction_series with 500 m3 pore volume 

147 >>> tedges_out = infiltration_to_extraction_series( 

148 ... flow=flow, 

149 ... tedges=tedges, 

150 ... aquifer_pore_volume=500.0, 

151 ... ) 

152 >>> len(tedges_out) 

153 11 

154 >>> # Time shift should be residence time = pore_volume / flow = 500 / 100 = 5 days 

155 >>> tedges_out[0] - tedges[0] 

156 Timedelta('5 days 00:00:00') 

157 

158 Plotting the input and output concentrations: 

159 

160 >>> import matplotlib.pyplot as plt 

161 >>> # Prepare data for step plot (repeat values for visualization) 

162 >>> xplot_in = np.repeat(tedges, 2)[1:-1] 

163 >>> yplot_in = np.repeat(cin, 2) 

164 >>> plt.plot( 

165 ... xplot_in, yplot_in, label="Concentration of infiltrated water" 

166 ... ) # doctest: +SKIP 

167 >>> 

168 >>> # cout equals cin, just with shifted time edges 

169 >>> xplot_out = np.repeat(tedges_out, 2)[1:-1] 

170 >>> yplot_out = np.repeat(cin, 2) 

171 >>> plt.plot( 

172 ... xplot_out, yplot_out, label="Concentration of extracted water" 

173 ... ) # doctest: +SKIP 

174 >>> plt.xlabel("Time") # doctest: +SKIP 

175 >>> plt.ylabel("Concentration") # doctest: +SKIP 

176 >>> plt.legend() # doctest: +SKIP 

177 >>> plt.show() # doctest: +SKIP 

178 

179 With retardation factor: 

180 

181 >>> tedges_out = infiltration_to_extraction_series( 

182 ... flow=flow, 

183 ... tedges=tedges, 

184 ... aquifer_pore_volume=500.0, 

185 ... retardation_factor=2.0, # Doubles residence time 

186 ... ) 

187 >>> # Time shift is doubled: 500 * 2.0 / 100 = 10 days 

188 >>> tedges_out[0] - tedges[0] 

189 Timedelta('10 days 00:00:00') 

190 """ 

191 rt_array = residence_time( 

192 flow=flow, 

193 flow_tedges=tedges, 

194 index=tedges, 

195 aquifer_pore_volumes=aquifer_pore_volume, 

196 retardation_factor=retardation_factor, 

197 direction="infiltration_to_extraction", 

198 ) 

199 return tedges + pd.to_timedelta(rt_array[0], unit="D", errors="coerce") 

200 

201 

202def extraction_to_infiltration_series( 

203 *, 

204 flow: npt.ArrayLike, 

205 tedges: pd.DatetimeIndex, 

206 aquifer_pore_volume: float, 

207 retardation_factor: float = 1.0, 

208) -> pd.DatetimeIndex: 

209 """ 

210 Compute infiltration time edges from extraction time edges (deconvolution). 

211 

212 This function shifts extraction time edges backward in time based on residence 

213 times computed from flow rates and aquifer properties. The concentration values remain 

214 unchanged (cin equals cout), only the time edges are shifted. This assumes a single pore 

215 volume (no distribution) and deterministic advective transport. This is the inverse 

216 operation of infiltration_to_extraction_series. 

217 

218 NOTE: This function is specifically designed for single aquifer pore volumes and does not 

219 support custom output time edges (cin_tedges). For distributions of aquifer pore volumes 

220 or custom output time grids, use `extraction_to_infiltration` instead. 

221 

222 Parameters 

223 ---------- 

224 flow : array-like 

225 Flow rate values in the aquifer [m3/day]. 

226 Length must match the number of time bins defined by tedges (len(tedges) - 1). 

227 tedges : pandas.DatetimeIndex 

228 Time edges defining bins for both cout and flow data. Has length of len(flow) + 1. 

229 aquifer_pore_volume : float 

230 Single aquifer pore volume [m3] used to compute residence times. 

231 retardation_factor : float, optional 

232 Retardation factor of the compound in the aquifer (default 1.0). 

233 Values > 1.0 indicate slower transport due to sorption/interaction. 

234 

235 Returns 

236 ------- 

237 pandas.DatetimeIndex 

238 Time edges for the infiltrating water concentration. Same length as tedges. 

239 The concentration values in the infiltrating water (cin) equal cout, but are 

240 aligned with these shifted time edges. 

241 

242 Examples 

243 -------- 

244 Basic usage with constant flow: 

245 

246 >>> import pandas as pd 

247 >>> import numpy as np 

248 >>> from gwtransport.utils import compute_time_edges 

249 >>> from gwtransport.advection import extraction_to_infiltration_series 

250 >>> 

251 >>> # Create input data 

252 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-10", freq="D") 

253 >>> tedges = compute_time_edges( 

254 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates) 

255 ... ) 

256 >>> 

257 >>> # Constant concentration and flow 

258 >>> cout = np.ones(len(dates)) * 10.0 

259 >>> flow = np.ones(len(dates)) * 100.0 # 100 m3/day 

260 >>> 

261 >>> # Run extraction_to_infiltration_series with 500 m3 pore volume 

262 >>> tedges_out = extraction_to_infiltration_series( 

263 ... flow=flow, 

264 ... tedges=tedges, 

265 ... aquifer_pore_volume=500.0, 

266 ... ) 

267 >>> len(tedges_out) 

268 11 

269 >>> # Time shift should be residence time = pore_volume / flow = 500 / 100 = 5 days (backward) 

270 >>> # First few elements are NaT due to insufficient history, check a valid index 

271 >>> tedges[5] - tedges_out[5] 

272 Timedelta('5 days 00:00:00') 

273 

274 Plotting the input and output concentrations: 

275 

276 >>> import matplotlib.pyplot as plt 

277 >>> # Prepare data for step plot (repeat values for visualization) 

278 >>> xplot_in = np.repeat(tedges, 2)[1:-1] 

279 >>> yplot_in = np.repeat(cout, 2) 

280 >>> plt.plot( 

281 ... xplot_in, yplot_in, label="Concentration of extracted water" 

282 ... ) # doctest: +SKIP 

283 >>> 

284 >>> # cin equals cout, just with shifted time edges 

285 >>> xplot_out = np.repeat(tedges_out, 2)[1:-1] 

286 >>> yplot_out = np.repeat(cout, 2) 

287 >>> plt.plot( 

288 ... xplot_out, yplot_out, label="Concentration of infiltrated water" 

289 ... ) # doctest: +SKIP 

290 >>> plt.xlabel("Time") # doctest: +SKIP 

291 >>> plt.ylabel("Concentration") # doctest: +SKIP 

292 >>> plt.legend() # doctest: +SKIP 

293 >>> plt.show() # doctest: +SKIP 

294 

295 With retardation factor: 

296 

297 >>> tedges_out = extraction_to_infiltration_series( 

298 ... flow=flow, 

299 ... tedges=tedges, 

300 ... aquifer_pore_volume=500.0, 

301 ... retardation_factor=2.0, # Doubles residence time 

302 ... ) 

303 >>> # Time shift is doubled: 500 * 2.0 / 100 = 10 days (backward) 

304 >>> # With longer residence time, more elements are NaT, check the last valid index 

305 >>> tedges[10] - tedges_out[10] 

306 Timedelta('10 days 00:00:00') 

307 """ 

308 rt_array = residence_time( 

309 flow=flow, 

310 flow_tedges=tedges, 

311 index=tedges, 

312 aquifer_pore_volumes=aquifer_pore_volume, 

313 retardation_factor=retardation_factor, 

314 direction="extraction_to_infiltration", 

315 ) 

316 return tedges - pd.to_timedelta(rt_array[0], unit="D", errors="coerce") 

317 

318 

319def gamma_infiltration_to_extraction( 

320 *, 

321 cin: npt.ArrayLike, 

322 flow: npt.ArrayLike, 

323 tedges: pd.DatetimeIndex, 

324 cout_tedges: pd.DatetimeIndex, 

325 alpha: float | None = None, 

326 beta: float | None = None, 

327 mean: float | None = None, 

328 std: float | None = None, 

329 n_bins: int = 100, 

330 retardation_factor: float = 1.0, 

331) -> npt.NDArray[np.floating]: 

332 """ 

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

334 

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

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

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

338 parameters alpha and beta. 

339 

340 This function represents infiltration to extraction modeling (equivalent to convolution). 

341 

342 Provide either alpha and beta or mean and std. 

343 

344 Parameters 

345 ---------- 

346 cin : array-like 

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

348 water. 

349 flow : array-like 

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

351 tedges : pandas.DatetimeIndex 

352 Time edges for both cin and flow data. Used to compute the cumulative concentration. 

353 Has a length of one more than `cin` and `flow`. 

354 cout_tedges : pandas.DatetimeIndex 

355 Time edges for the output data. Used to compute the cumulative concentration. 

356 Has a length of one more than the desired output length. 

357 alpha : float, optional 

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

359 beta : float, optional 

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

361 mean : float, optional 

362 Mean of the gamma distribution. 

363 std : float, optional 

364 Standard deviation of the gamma distribution. 

365 n_bins : int 

366 Number of bins to discretize the gamma distribution. 

367 retardation_factor : float 

368 Retardation factor of the compound in the aquifer. 

369 

370 Returns 

371 ------- 

372 numpy.ndarray 

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

374 

375 See Also 

376 -------- 

377 infiltration_to_extraction : Transport with explicit pore volume distribution 

378 gamma_extraction_to_infiltration : Reverse operation (deconvolution) 

379 gwtransport.gamma.bins : Create gamma distribution bins 

380 gwtransport.residence_time.residence_time : Compute residence times 

381 gwtransport.diffusion.infiltration_to_extraction : Add microdispersion and molecular diffusion 

382 :ref:`concept-gamma-distribution` : Two-parameter pore volume model 

383 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate 

384 

385 Notes 

386 ----- 

387 The APVD is only time-invariant under the steady-streamlines assumption 

388 (see :ref:`assumption-steady-streamlines`). 

389 

390 The spreading from the gamma-distributed pore volumes represents macrodispersion 

391 (aquifer-scale heterogeneity). When ``std`` comes from calibration on measurements, 

392 it absorbs all mixing: macrodispersion, microdispersion, and an average molecular 

393 diffusion contribution. When calibrating with the diffusion module, these three 

394 components are taken into account separately. When ``std`` comes from streamline 

395 analysis, it represents macrodispersion only; microdispersion and molecular diffusion 

396 can be added via the diffusion module or by combining variances 

397 (see :doc:`/examples/05_Diffusion_Dispersion`). 

398 

399 If calibrating with one compound (e.g., temperature) and predicting for another 

400 (e.g., a solute), the baked-in molecular diffusion contribution may need 

401 correction — see :doc:`/examples/05_Diffusion_Dispersion`. 

402 See :ref:`concept-dispersion-scales` for guidance on when to add microdispersion 

403 using the diffusion module. 

404 

405 Examples 

406 -------- 

407 Basic usage with alpha and beta parameters: 

408 

409 >>> import pandas as pd 

410 >>> import numpy as np 

411 >>> from gwtransport.utils import compute_time_edges 

412 >>> from gwtransport.advection import gamma_infiltration_to_extraction 

413 >>> 

414 >>> # Create input data with aligned time edges 

415 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D") 

416 >>> tedges = compute_time_edges( 

417 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates) 

418 ... ) 

419 >>> 

420 >>> # Create output time edges (can be different alignment) 

421 >>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D") 

422 >>> cout_tedges = compute_time_edges( 

423 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates) 

424 ... ) 

425 >>> 

426 >>> # Input concentration and flow (same length, aligned with tedges) 

427 >>> cin = pd.Series(np.ones(len(dates)), index=dates) 

428 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day 

429 >>> 

430 >>> # Run gamma_infiltration_to_extraction with alpha/beta parameters 

431 >>> cout = gamma_infiltration_to_extraction( 

432 ... cin=cin, 

433 ... flow=flow, 

434 ... tedges=tedges, 

435 ... cout_tedges=cout_tedges, 

436 ... alpha=10.0, 

437 ... beta=10.0, 

438 ... n_bins=5, 

439 ... ) 

440 >>> cout.shape 

441 (11,) 

442 

443 Using mean and std parameters instead: 

444 

445 >>> cout = gamma_infiltration_to_extraction( 

446 ... cin=cin, 

447 ... flow=flow, 

448 ... tedges=tedges, 

449 ... cout_tedges=cout_tedges, 

450 ... mean=100.0, 

451 ... std=20.0, 

452 ... n_bins=5, 

453 ... ) 

454 

455 With retardation factor: 

456 

457 >>> cout = gamma_infiltration_to_extraction( 

458 ... cin=cin, 

459 ... flow=flow, 

460 ... tedges=tedges, 

461 ... cout_tedges=cout_tedges, 

462 ... alpha=10.0, 

463 ... beta=10.0, 

464 ... retardation_factor=2.0, # Doubles residence time 

465 ... ) 

466 """ 

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

468 return infiltration_to_extraction( 

469 cin=cin, 

470 flow=flow, 

471 tedges=tedges, 

472 cout_tedges=cout_tedges, 

473 aquifer_pore_volumes=bins["expected_values"], 

474 retardation_factor=retardation_factor, 

475 ) 

476 

477 

478def gamma_extraction_to_infiltration( 

479 *, 

480 cout: npt.ArrayLike, 

481 flow: npt.ArrayLike, 

482 tedges: pd.DatetimeIndex, 

483 cin_tedges: pd.DatetimeIndex, 

484 alpha: float | None = None, 

485 beta: float | None = None, 

486 mean: float | None = None, 

487 std: float | None = None, 

488 n_bins: int = 100, 

489 retardation_factor: float = 1.0, 

490) -> npt.NDArray[np.floating]: 

491 """ 

492 Compute the concentration of the infiltrating water from extracted water (deconvolution). 

493 

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

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

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

497 parameters alpha and beta. 

498 

499 This function represents extraction to infiltration modeling (equivalent to deconvolution). 

500 It is symmetric to gamma_infiltration_to_extraction. 

501 

502 Provide either alpha and beta or mean and std. 

503 

504 Parameters 

505 ---------- 

506 cout : array-like 

507 Concentration of the compound in extracted water or temperature of extracted 

508 water. 

509 tedges : pandas.DatetimeIndex 

510 Time edges for the cout and flow data. Used to compute the cumulative concentration. 

511 Has a length of one more than `cout` and `flow`. 

512 cin_tedges : pandas.DatetimeIndex 

513 Time edges for the output (infiltration) data. Used to compute the cumulative concentration. 

514 Has a length of one more than the desired output length. 

515 flow : array-like 

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

517 alpha : float, optional 

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

519 beta : float, optional 

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

521 mean : float, optional 

522 Mean of the gamma distribution. 

523 std : float, optional 

524 Standard deviation of the gamma distribution. 

525 n_bins : int 

526 Number of bins to discretize the gamma distribution. 

527 retardation_factor : float, optional 

528 Retardation factor of the compound in the aquifer (default 1.0). 

529 Values > 1.0 indicate slower transport due to sorption/interaction. 

530 

531 Returns 

532 ------- 

533 numpy.ndarray 

534 Concentration of the compound in the infiltrating water [ng/m3] or temperature. 

535 

536 See Also 

537 -------- 

538 extraction_to_infiltration : Deconvolution with explicit pore volume distribution 

539 gamma_infiltration_to_extraction : Forward operation (convolution) 

540 gwtransport.gamma.bins : Create gamma distribution bins 

541 gwtransport.diffusion.extraction_to_infiltration : Deconvolution with microdispersion and molecular diffusion 

542 :ref:`concept-gamma-distribution` : Two-parameter pore volume model 

543 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate 

544 

545 Notes 

546 ----- 

547 The APVD is only time-invariant under the steady-streamlines assumption 

548 (see :ref:`assumption-steady-streamlines`). 

549 

550 The spreading from the gamma-distributed pore volumes represents macrodispersion 

551 (aquifer-scale heterogeneity). When ``std`` comes from calibration on measurements, 

552 it absorbs all mixing: macrodispersion, microdispersion, and an average molecular 

553 diffusion contribution. When calibrating with the diffusion module, these three 

554 components are taken into account separately. When ``std`` comes from streamline 

555 analysis, it represents macrodispersion only; microdispersion and molecular diffusion 

556 can be added via the diffusion module or by combining variances 

557 (see :doc:`/examples/05_Diffusion_Dispersion`). 

558 

559 If calibrating with one compound (e.g., temperature) and predicting for another 

560 (e.g., a solute), the baked-in molecular diffusion contribution may need 

561 correction — see :doc:`/examples/05_Diffusion_Dispersion`. 

562 See :ref:`concept-dispersion-scales` for guidance on when to add microdispersion 

563 using the diffusion module. 

564 

565 Examples 

566 -------- 

567 Basic usage with alpha and beta parameters: 

568 

569 >>> import pandas as pd 

570 >>> import numpy as np 

571 >>> from gwtransport.utils import compute_time_edges 

572 >>> from gwtransport.advection import gamma_extraction_to_infiltration 

573 >>> 

574 >>> # Create input data with aligned time edges 

575 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D") 

576 >>> tedges = compute_time_edges( 

577 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates) 

578 ... ) 

579 >>> 

580 >>> # Create output time edges (can be different alignment) 

581 >>> cin_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D") 

582 >>> cin_tedges = compute_time_edges( 

583 ... tedges=None, tstart=None, tend=cin_dates, number_of_bins=len(cin_dates) 

584 ... ) 

585 >>> 

586 >>> # Input concentration and flow (same length, aligned with tedges) 

587 >>> cout = pd.Series(np.ones(len(dates)), index=dates) 

588 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day 

589 >>> 

590 >>> # Run gamma_extraction_to_infiltration with alpha/beta parameters 

591 >>> cin = gamma_extraction_to_infiltration( 

592 ... cout=cout, 

593 ... tedges=tedges, 

594 ... cin_tedges=cin_tedges, 

595 ... flow=flow, 

596 ... alpha=10.0, 

597 ... beta=10.0, 

598 ... n_bins=5, 

599 ... ) 

600 >>> cin.shape 

601 (22,) 

602 

603 Using mean and std parameters instead: 

604 

605 >>> cin = gamma_extraction_to_infiltration( 

606 ... cout=cout, 

607 ... tedges=tedges, 

608 ... cin_tedges=cin_tedges, 

609 ... flow=flow, 

610 ... mean=100.0, 

611 ... std=20.0, 

612 ... n_bins=5, 

613 ... ) 

614 

615 With retardation factor: 

616 

617 >>> cin = gamma_extraction_to_infiltration( 

618 ... cout=cout, 

619 ... tedges=tedges, 

620 ... cin_tedges=cin_tedges, 

621 ... flow=flow, 

622 ... alpha=10.0, 

623 ... beta=10.0, 

624 ... retardation_factor=2.0, # Doubles residence time 

625 ... ) 

626 """ 

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

628 return extraction_to_infiltration( 

629 cout=cout, 

630 flow=flow, 

631 tedges=tedges, 

632 cin_tedges=cin_tedges, 

633 aquifer_pore_volumes=bins["expected_values"], 

634 retardation_factor=retardation_factor, 

635 ) 

636 

637 

638def infiltration_to_extraction( 

639 *, 

640 cin: npt.ArrayLike, 

641 flow: npt.ArrayLike, 

642 tedges: pd.DatetimeIndex, 

643 cout_tedges: pd.DatetimeIndex, 

644 aquifer_pore_volumes: npt.ArrayLike, 

645 retardation_factor: float = 1.0, 

646) -> npt.NDArray[np.floating]: 

647 """ 

648 Compute the concentration of the extracted water using flow-weighted advection. 

649 

650 This function implements an infiltration to extraction advection model where cin and flow values 

651 correspond to the same aligned time bins defined by tedges. 

652 

653 The algorithm: 

654 1. Computes residence times for each pore volume at cout time edges 

655 2. Calculates infiltration time edges by subtracting residence times 

656 3. Determines temporal overlaps between infiltration and cin time windows 

657 4. Creates flow-weighted overlap matrices normalized by total weights 

658 5. Computes weighted contributions and averages across pore volumes 

659 

660 

661 Parameters 

662 ---------- 

663 cin : array-like 

664 Concentration values of infiltrating water or temperature [concentration units]. 

665 Length must match the number of time bins defined by tedges. 

666 flow : array-like 

667 Flow rate values in the aquifer [m3/day]. 

668 Length must match cin and the number of time bins defined by tedges. 

669 tedges : pandas.DatetimeIndex 

670 Time edges defining bins for both cin and flow data. Has length of 

671 len(cin) + 1 and len(flow) + 1. 

672 cout_tedges : pandas.DatetimeIndex 

673 Time edges for output data bins. Has length of desired output + 1. 

674 Can have different time alignment and resolution than tedges. 

675 aquifer_pore_volumes : array-like 

676 Array of aquifer pore volumes [m3] representing the distribution 

677 of residence times in the aquifer system. 

678 retardation_factor : float, optional 

679 Retardation factor of the compound in the aquifer (default 1.0). 

680 Values > 1.0 indicate slower transport due to sorption/interaction. 

681 

682 Returns 

683 ------- 

684 numpy.ndarray 

685 Flow-weighted concentration in the extracted water. Same units as cin. 

686 Length equals len(cout_tedges) - 1. NaN values indicate time periods 

687 with no valid contributions from the infiltration data. 

688 

689 Raises 

690 ------ 

691 ValueError 

692 If tedges length doesn't match cin/flow arrays plus one, or if 

693 infiltration time edges become non-monotonic (invalid input conditions). 

694 

695 See Also 

696 -------- 

697 gamma_infiltration_to_extraction : Transport with gamma-distributed pore volumes 

698 extraction_to_infiltration : Reverse operation (deconvolution) 

699 infiltration_to_extraction_series : Simple time-shift for single pore volume 

700 gwtransport.residence_time.residence_time : Compute residence times from flow and pore volume 

701 gwtransport.residence_time.freundlich_retardation : Compute concentration-dependent retardation 

702 :ref:`concept-pore-volume-distribution` : Background on aquifer heterogeneity modeling 

703 :ref:`concept-transport-equation` : Flow-weighted averaging approach 

704 

705 Examples 

706 -------- 

707 Basic usage with pandas Series: 

708 

709 >>> import pandas as pd 

710 >>> import numpy as np 

711 >>> from gwtransport.utils import compute_time_edges 

712 >>> from gwtransport.advection import infiltration_to_extraction 

713 >>> 

714 >>> # Create input data 

715 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D") 

716 >>> tedges = compute_time_edges( 

717 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates) 

718 ... ) 

719 >>> 

720 >>> # Create output time edges (different alignment) 

721 >>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D") 

722 >>> cout_tedges = compute_time_edges( 

723 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates) 

724 ... ) 

725 >>> 

726 >>> # Input concentration and flow 

727 >>> cin = pd.Series(np.ones(len(dates)), index=dates) 

728 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day 

729 >>> 

730 >>> # Define distribution of aquifer pore volumes 

731 >>> aquifer_pore_volumes = np.array([50, 100, 200]) # m3 

732 >>> 

733 >>> # Run infiltration_to_extraction 

734 >>> cout = infiltration_to_extraction( 

735 ... cin=cin, 

736 ... flow=flow, 

737 ... tedges=tedges, 

738 ... cout_tedges=cout_tedges, 

739 ... aquifer_pore_volumes=aquifer_pore_volumes, 

740 ... ) 

741 >>> cout.shape 

742 (11,) 

743 

744 Using array inputs instead of pandas Series: 

745 

746 >>> # Convert to arrays 

747 >>> cin_values = cin.values 

748 >>> flow_values = flow.values 

749 >>> 

750 >>> cout = infiltration_to_extraction( 

751 ... cin=cin_values, 

752 ... flow=flow, 

753 ... tedges=tedges, 

754 ... cout_tedges=cout_tedges, 

755 ... aquifer_pore_volumes=aquifer_pore_volumes, 

756 ... ) 

757 

758 With constant retardation factor (linear sorption): 

759 

760 >>> cout = infiltration_to_extraction( 

761 ... cin=cin, 

762 ... flow=flow, 

763 ... tedges=tedges, 

764 ... cout_tedges=cout_tedges, 

765 ... aquifer_pore_volumes=aquifer_pore_volumes, 

766 ... retardation_factor=2.0, # Compound moves twice as slowly 

767 ... ) 

768 

769 Note: For concentration-dependent retardation (nonlinear sorption), 

770 use `infiltration_to_extraction_front_tracking_detailed` instead, as this 

771 function only supports constant (float) retardation factors. 

772 

773 Using single pore volume: 

774 

775 >>> single_volume = np.array([100]) # Single 100 m3 pore volume 

776 >>> cout = infiltration_to_extraction( 

777 ... cin=cin, 

778 ... flow=flow, 

779 ... tedges=tedges, 

780 ... cout_tedges=cout_tedges, 

781 ... aquifer_pore_volumes=single_volume, 

782 ... ) 

783 """ 

784 tedges = pd.DatetimeIndex(tedges) 

785 cout_tedges = pd.DatetimeIndex(cout_tedges) 

786 

787 # Convert to arrays for vectorized operations 

788 cin = np.asarray(cin) 

789 flow = np.asarray(flow) 

790 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes) 

791 

792 if len(tedges) != len(cin) + 1: 

793 msg = "tedges must have one more element than cin" 

794 raise ValueError(msg) 

795 if len(tedges) != len(flow) + 1: 

796 msg = "tedges must have one more element than flow" 

797 raise ValueError(msg) 

798 

799 # Validate inputs do not contain NaN values 

800 if np.any(np.isnan(cin)): 

801 msg = "cin contains NaN values, which are not allowed" 

802 raise ValueError(msg) 

803 if np.any(np.isnan(flow)): 

804 msg = "flow contains NaN values, which are not allowed" 

805 raise ValueError(msg) 

806 if retardation_factor < 1.0: 

807 msg = "retardation_factor must be >= 1.0" 

808 raise ValueError(msg) 

809 

810 # Compute normalized weights (includes all pre-computation) 

811 normalized_weights = _infiltration_to_extraction_weights( 

812 tedges=tedges, 

813 cout_tedges=cout_tedges, 

814 aquifer_pore_volumes=aquifer_pore_volumes, 

815 flow=flow, 

816 retardation_factor=retardation_factor, 

817 ) 

818 

819 # Apply to concentrations and handle NaN for periods with no contributions 

820 out = normalized_weights.dot(cin) 

821 # Set NaN where no valid pore volumes contributed 

822 total_weights = np.sum(normalized_weights, axis=1) 

823 out[total_weights == 0] = np.nan 

824 

825 return out 

826 

827 

828def extraction_to_infiltration( 

829 *, 

830 cout: npt.ArrayLike, 

831 flow: npt.ArrayLike, 

832 tedges: pd.DatetimeIndex, 

833 cin_tedges: pd.DatetimeIndex, 

834 aquifer_pore_volumes: npt.ArrayLike, 

835 retardation_factor: float = 1.0, 

836) -> npt.NDArray[np.floating]: 

837 """ 

838 Compute the concentration of the infiltrating water from extracted water (deconvolution). 

839 

840 This function implements an extraction to infiltration advection model (inverse of infiltration_to_extraction) 

841 where cout and flow values correspond to the same aligned time bins defined by tedges. 

842 

843 SYMMETRIC RELATIONSHIP: 

844 - infiltration_to_extraction: cin + tedges → cout + cout_tedges 

845 - extraction_to_infiltration: cout + tedges → cin + cin_tedges 

846 

847 The algorithm (symmetric to infiltration_to_extraction): 

848 1. Computes residence times for each pore volume at cint time edges 

849 2. Calculates extraction time edges by adding residence times (reverse of infiltration_to_extraction) 

850 3. Determines temporal overlaps between extraction and cout time windows 

851 4. Creates flow-weighted overlap matrices normalized by total weights 

852 5. Computes weighted contributions and averages across pore volumes 

853 

854 

855 Parameters 

856 ---------- 

857 cout : array-like 

858 Concentration values of extracted water [concentration units]. 

859 Length must match the number of time bins defined by tedges. 

860 flow : array-like 

861 Flow rate values in the aquifer [m3/day]. 

862 Length must match cout and the number of time bins defined by tedges. 

863 tedges : pandas.DatetimeIndex 

864 Time edges defining bins for both cout and flow data. Has length of 

865 len(cout) + 1 and len(flow) + 1. 

866 cin_tedges : pandas.DatetimeIndex 

867 Time edges for output (infiltration) data bins. Has length of desired output + 1. 

868 Can have different time alignment and resolution than tedges. 

869 aquifer_pore_volumes : array-like 

870 Array of aquifer pore volumes [m3] representing the distribution 

871 of residence times in the aquifer system. 

872 retardation_factor : float, optional 

873 Retardation factor of the compound in the aquifer (default 1.0). 

874 Values > 1.0 indicate slower transport due to sorption/interaction. 

875 

876 Returns 

877 ------- 

878 numpy.ndarray 

879 Flow-weighted concentration in the infiltrating water. Same units as cout. 

880 Length equals len(cin_tedges) - 1. NaN values indicate time periods 

881 with no valid contributions from the extraction data. 

882 

883 Raises 

884 ------ 

885 ValueError 

886 If tedges length doesn't match cout/flow arrays plus one, or if 

887 extraction time edges become non-monotonic (invalid input conditions). 

888 

889 See Also 

890 -------- 

891 gamma_extraction_to_infiltration : Deconvolution with gamma-distributed pore volumes 

892 infiltration_to_extraction : Forward operation (convolution) 

893 extraction_to_infiltration_series : Simple time-shift for single pore volume 

894 gwtransport.residence_time.residence_time : Compute residence times from flow and pore volume 

895 gwtransport.residence_time.freundlich_retardation : Compute concentration-dependent retardation 

896 :ref:`concept-pore-volume-distribution` : Background on aquifer heterogeneity modeling 

897 :ref:`concept-transport-equation` : Flow-weighted averaging approach 

898 

899 Examples 

900 -------- 

901 Basic usage with pandas Series: 

902 

903 >>> import pandas as pd 

904 >>> import numpy as np 

905 >>> from gwtransport.utils import compute_time_edges 

906 >>> from gwtransport.advection import extraction_to_infiltration 

907 >>> 

908 >>> # Create input data 

909 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D") 

910 >>> tedges = compute_time_edges( 

911 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates) 

912 ... ) 

913 >>> 

914 >>> # Create output time edges (different alignment) 

915 >>> cint_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D") 

916 >>> cin_tedges = compute_time_edges( 

917 ... tedges=None, tstart=None, tend=cint_dates, number_of_bins=len(cint_dates) 

918 ... ) 

919 >>> 

920 >>> # Input concentration and flow 

921 >>> cout = pd.Series(np.ones(len(dates)), index=dates) 

922 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day 

923 >>> 

924 >>> # Define distribution of aquifer pore volumes 

925 >>> aquifer_pore_volumes = np.array([50, 100, 200]) # m3 

926 >>> 

927 >>> # Run extraction_to_infiltration 

928 >>> cin = extraction_to_infiltration( 

929 ... cout=cout, 

930 ... flow=flow, 

931 ... tedges=tedges, 

932 ... cin_tedges=cin_tedges, 

933 ... aquifer_pore_volumes=aquifer_pore_volumes, 

934 ... ) 

935 >>> cin.shape 

936 (22,) 

937 

938 Using array inputs instead of pandas Series: 

939 

940 >>> # Convert to arrays 

941 >>> cout = cout.values 

942 >>> flow = flow.values 

943 >>> 

944 >>> cin = extraction_to_infiltration( 

945 ... cout=cout, 

946 ... flow=flow, 

947 ... tedges=tedges, 

948 ... cin_tedges=cin_tedges, 

949 ... aquifer_pore_volumes=aquifer_pore_volumes, 

950 ... ) 

951 

952 With retardation factor: 

953 

954 >>> cin = extraction_to_infiltration( 

955 ... cout=cout, 

956 ... flow=flow, 

957 ... tedges=tedges, 

958 ... cin_tedges=cin_tedges, 

959 ... aquifer_pore_volumes=aquifer_pore_volumes, 

960 ... retardation_factor=2.0, # Compound moves twice as slowly 

961 ... ) 

962 

963 Using single pore volume: 

964 

965 >>> single_volume = np.array([100]) # Single 100 m3 pore volume 

966 >>> cin = extraction_to_infiltration( 

967 ... cout=cout, 

968 ... flow=flow, 

969 ... tedges=tedges, 

970 ... cin_tedges=cin_tedges, 

971 ... aquifer_pore_volumes=single_volume, 

972 ... ) 

973 

974 """ 

975 tedges = pd.DatetimeIndex(tedges) 

976 cin_tedges = pd.DatetimeIndex(cin_tedges) 

977 

978 # Convert to arrays for vectorized operations 

979 cout = np.asarray(cout) 

980 flow = np.asarray(flow) 

981 

982 if len(tedges) != len(cout) + 1: 

983 msg = "tedges must have one more element than cout" 

984 raise ValueError(msg) 

985 if len(tedges) != len(flow) + 1: 

986 msg = "tedges must have one more element than flow" 

987 raise ValueError(msg) 

988 

989 # Validate inputs do not contain NaN values 

990 if np.any(np.isnan(cout)): 

991 msg = "cout contains NaN values, which are not allowed" 

992 raise ValueError(msg) 

993 if np.any(np.isnan(flow)): 

994 msg = "flow contains NaN values, which are not allowed" 

995 raise ValueError(msg) 

996 if retardation_factor < 1.0: 

997 msg = "retardation_factor must be >= 1.0" 

998 raise ValueError(msg) 

999 

1000 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes) 

1001 

1002 # Compute normalized weights (includes all pre-computation) 

1003 normalized_weights = _extraction_to_infiltration_weights( 

1004 tedges=tedges, 

1005 cin_tedges=cin_tedges, 

1006 aquifer_pore_volumes=aquifer_pore_volumes, 

1007 flow=flow, 

1008 retardation_factor=retardation_factor, 

1009 ) 

1010 

1011 # Apply to concentrations and handle NaN for periods with no contributions 

1012 out = normalized_weights.dot(cout) 

1013 # Set NaN where no valid pore volumes contributed 

1014 total_weights = np.sum(normalized_weights, axis=1) 

1015 out[total_weights == 0] = np.nan 

1016 

1017 return out 

1018 

1019 

1020def _validate_front_tracking_inputs( 

1021 *, 

1022 cin: npt.ArrayLike, 

1023 flow: npt.ArrayLike, 

1024 tedges: pd.DatetimeIndex, 

1025 cout_tedges: pd.DatetimeIndex, 

1026 aquifer_pore_volumes: npt.ArrayLike, 

1027 freundlich_k: float | None, 

1028 freundlich_n: float | None, 

1029 bulk_density: float | None, 

1030 porosity: float | None, 

1031 retardation_factor: float | None, 

1032) -> tuple: 

1033 """Validate inputs and create sorption object for front tracking functions.""" 

1034 cin = np.asarray(cin, dtype=float) 

1035 flow = np.asarray(flow, dtype=float) 

1036 tedges = pd.DatetimeIndex(tedges) 

1037 cout_tedges = pd.DatetimeIndex(cout_tedges) 

1038 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes, dtype=float) 

1039 

1040 if len(tedges) != len(cin) + 1: 

1041 msg = "tedges must have length len(cin) + 1" 

1042 raise ValueError(msg) 

1043 if len(flow) != len(cin): 

1044 msg = "flow must have same length as cin" 

1045 raise ValueError(msg) 

1046 if np.any(cin < 0): 

1047 msg = "cin must be non-negative" 

1048 raise ValueError(msg) 

1049 if np.any(flow <= 0): 

1050 msg = "flow must be positive" 

1051 raise ValueError(msg) 

1052 if np.any(np.isnan(cin)) or np.any(np.isnan(flow)): 

1053 msg = "cin and flow must not contain NaN" 

1054 raise ValueError(msg) 

1055 if np.any(aquifer_pore_volumes <= 0): 

1056 msg = "aquifer_pore_volumes must be positive" 

1057 raise ValueError(msg) 

1058 

1059 # Convert cout_tedges to days (relative to tedges[0]) for output computation 

1060 t_ref = tedges[0] 

1061 cout_tedges_days = ((cout_tedges - t_ref) / pd.Timedelta(days=1)).values 

1062 

1063 # Create sorption object 

1064 if retardation_factor is not None: 

1065 if retardation_factor < 1.0: 

1066 msg = "retardation_factor must be >= 1.0" 

1067 raise ValueError(msg) 

1068 

1069 sorption = ConstantRetardation(retardation_factor=retardation_factor) 

1070 else: 

1071 if freundlich_k is None or freundlich_n is None or bulk_density is None or porosity is None: 

1072 msg = ( 

1073 "Must provide either retardation_factor or all Freundlich parameters " 

1074 "(freundlich_k, freundlich_n, bulk_density, porosity)" 

1075 ) 

1076 raise ValueError(msg) 

1077 if freundlich_k <= 0 or freundlich_n <= 0: 

1078 msg = "Freundlich parameters must be positive" 

1079 raise ValueError(msg) 

1080 if abs(freundlich_n - 1.0) < EPSILON_FREUNDLICH_N: 

1081 msg = "freundlich_n = 1 not supported (use retardation_factor for linear case)" 

1082 raise ValueError(msg) 

1083 if bulk_density <= 0 or not 0 < porosity < 1: 

1084 msg = "Invalid physical parameters" 

1085 raise ValueError(msg) 

1086 

1087 sorption = FreundlichSorption( 

1088 k_f=freundlich_k, 

1089 n=freundlich_n, 

1090 bulk_density=bulk_density, 

1091 porosity=porosity, 

1092 ) 

1093 

1094 return cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days 

1095 

1096 

1097def infiltration_to_extraction_front_tracking( 

1098 *, 

1099 cin: npt.ArrayLike, 

1100 flow: npt.ArrayLike, 

1101 tedges: pd.DatetimeIndex, 

1102 cout_tedges: pd.DatetimeIndex, 

1103 aquifer_pore_volumes: npt.ArrayLike, 

1104 freundlich_k: float | None = None, 

1105 freundlich_n: float | None = None, 

1106 bulk_density: float | None = None, 

1107 porosity: float | None = None, 

1108 retardation_factor: float | None = None, 

1109 max_iterations: int = 10000, 

1110) -> npt.NDArray[np.floating]: 

1111 """ 

1112 Compute extracted concentration using exact front tracking with nonlinear sorption. 

1113 

1114 Uses event-driven analytical algorithm that tracks shock waves, rarefaction waves, 

1115 and characteristics with machine precision. No numerical dispersion, exact mass 

1116 balance to floating-point precision. 

1117 

1118 Parameters 

1119 ---------- 

1120 cin : array-like 

1121 Infiltration concentration [mg/L or any units]. 

1122 Length = len(tedges) - 1. 

1123 flow : array-like 

1124 Flow rate [m³/day]. Must be positive. 

1125 Length = len(tedges) - 1. 

1126 tedges : pandas.DatetimeIndex 

1127 Time bin edges. Length = len(cin) + 1. 

1128 cout_tedges : pandas.DatetimeIndex 

1129 Output time bin edges. Can be different from tedges. 

1130 Length determines output array size. 

1131 aquifer_pore_volumes : array-like 

1132 Array of aquifer pore volumes [m³] representing the distribution 

1133 of residence times in the aquifer system. Each pore volume must be positive. 

1134 freundlich_k : float, optional 

1135 Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive. 

1136 Used if retardation_factor is None. 

1137 freundlich_n : float, optional 

1138 Freundlich exponent [-]. Must be positive and != 1. 

1139 Used if retardation_factor is None. 

1140 bulk_density : float, optional 

1141 Bulk density [kg/m³]. Must be positive. 

1142 Used if retardation_factor is None. 

1143 porosity : float, optional 

1144 Porosity [-]. Must be in (0, 1). 

1145 Used if retardation_factor is None. 

1146 retardation_factor : float, optional 

1147 Constant retardation factor [-]. If provided, uses linear retardation 

1148 instead of Freundlich sorption. Must be >= 1.0. 

1149 max_iterations : int, optional 

1150 Maximum number of events. Default 10000. 

1151 

1152 Returns 

1153 ------- 

1154 cout : numpy.ndarray 

1155 Bin-averaged extraction concentration averaged across all pore volumes. 

1156 Length = len(cout_tedges) - 1. 

1157 

1158 Notes 

1159 ----- 

1160 **Spin-up Period**: 

1161 The function computes the first arrival time t_first. Concentrations 

1162 before t_first are affected by unknown initial conditions and should 

1163 not be used for analysis. Use `infiltration_to_extraction_front_tracking_detailed` 

1164 to access t_first. 

1165 

1166 **Machine Precision**: 

1167 All calculations use exact analytical formulas. Mass balance is conserved 

1168 to floating-point precision (~1e-14 relative error). No numerical tolerances 

1169 are used for time/position calculations. 

1170 

1171 **Physical Correctness**: 

1172 - All shocks satisfy Lax entropy condition 

1173 - Rarefaction waves use self-similar solutions 

1174 - Causality is strictly enforced 

1175 

1176 Examples 

1177 -------- 

1178 >>> import numpy as np 

1179 >>> import pandas as pd 

1180 >>> 

1181 >>> # Pulse injection with single pore volume 

1182 >>> tedges = pd.date_range("2020-01-01", periods=4, freq="10D") 

1183 >>> cin = np.array([0.0, 10.0, 0.0]) 

1184 >>> flow = np.array([100.0, 100.0, 100.0]) 

1185 >>> cout_tedges = pd.date_range("2020-01-01", periods=10, freq="5D") 

1186 >>> 

1187 >>> cout = infiltration_to_extraction_front_tracking( 

1188 ... cin=cin, 

1189 ... flow=flow, 

1190 ... tedges=tedges, 

1191 ... cout_tedges=cout_tedges, 

1192 ... aquifer_pore_volumes=np.array([500.0]), 

1193 ... freundlich_k=0.01, 

1194 ... freundlich_n=2.0, 

1195 ... bulk_density=1500.0, 

1196 ... porosity=0.3, 

1197 ... ) 

1198 

1199 With multiple pore volumes (distribution): 

1200 

1201 >>> aquifer_pore_volumes = np.array([400.0, 500.0, 600.0]) 

1202 >>> cout = infiltration_to_extraction_front_tracking( 

1203 ... cin=cin, 

1204 ... flow=flow, 

1205 ... tedges=tedges, 

1206 ... cout_tedges=cout_tedges, 

1207 ... aquifer_pore_volumes=aquifer_pore_volumes, 

1208 ... freundlich_k=0.01, 

1209 ... freundlich_n=2.0, 

1210 ... bulk_density=1500.0, 

1211 ... porosity=0.3, 

1212 ... ) 

1213 

1214 See Also 

1215 -------- 

1216 infiltration_to_extraction_front_tracking_detailed : Returns detailed structure 

1217 infiltration_to_extraction : Convolution-based approach for linear case 

1218 gamma_infiltration_to_extraction : For distributions of pore volumes 

1219 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and front-tracking theory 

1220 :ref:`assumption-advection-dominated` : When diffusion/dispersion is negligible 

1221 """ 

1222 cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days = _validate_front_tracking_inputs( 

1223 cin=cin, 

1224 flow=flow, 

1225 tedges=tedges, 

1226 cout_tedges=cout_tedges, 

1227 aquifer_pore_volumes=aquifer_pore_volumes, 

1228 freundlich_k=freundlich_k, 

1229 freundlich_n=freundlich_n, 

1230 bulk_density=bulk_density, 

1231 porosity=porosity, 

1232 retardation_factor=retardation_factor, 

1233 ) 

1234 

1235 # Loop over each pore volume and compute concentration 

1236 cout_all = np.zeros((len(aquifer_pore_volumes), len(cout_tedges) - 1)) 

1237 

1238 for i, aquifer_pore_volume in enumerate(aquifer_pore_volumes): 

1239 # Create tracker and run simulation for this pore volume 

1240 tracker = FrontTracker( 

1241 cin=cin, 

1242 flow=flow, 

1243 tedges=tedges, 

1244 aquifer_pore_volume=aquifer_pore_volume, 

1245 sorption=sorption, 

1246 ) 

1247 

1248 tracker.run(max_iterations=max_iterations) 

1249 

1250 # Extract bin-averaged concentrations at outlet for this pore volume 

1251 cout_all[i, :] = compute_bin_averaged_concentration_exact( 

1252 t_edges=cout_tedges_days, 

1253 v_outlet=aquifer_pore_volume, 

1254 waves=tracker.state.waves, 

1255 sorption=sorption, 

1256 ) 

1257 

1258 # Return average across all pore volumes 

1259 return np.mean(cout_all, axis=0) 

1260 

1261 

1262def infiltration_to_extraction_front_tracking_detailed( 

1263 *, 

1264 cin: npt.ArrayLike, 

1265 flow: npt.ArrayLike, 

1266 tedges: pd.DatetimeIndex, 

1267 cout_tedges: pd.DatetimeIndex, 

1268 aquifer_pore_volumes: npt.ArrayLike, 

1269 freundlich_k: float | None = None, 

1270 freundlich_n: float | None = None, 

1271 bulk_density: float | None = None, 

1272 porosity: float | None = None, 

1273 retardation_factor: float | None = None, 

1274 max_iterations: int = 10000, 

1275) -> tuple[npt.NDArray[np.floating], list[dict]]: 

1276 """ 

1277 Compute extracted concentration with complete diagnostic information. 

1278 

1279 Returns both bin-averaged concentrations and detailed simulation structure for each pore volume. 

1280 

1281 Parameters 

1282 ---------- 

1283 cin : array-like 

1284 Infiltration concentration [mg/L or any units]. 

1285 Length = len(tedges) - 1. 

1286 flow : array-like 

1287 Flow rate [m³/day]. Must be positive. 

1288 Length = len(tedges) - 1. 

1289 tedges : pandas.DatetimeIndex 

1290 Time bin edges. Length = len(cin) + 1. 

1291 cout_tedges : pandas.DatetimeIndex 

1292 Output time bin edges. Can be different from tedges. 

1293 Length determines output array size. 

1294 aquifer_pore_volumes : array-like 

1295 Array of aquifer pore volumes [m³] representing the distribution 

1296 of residence times in the aquifer system. Each pore volume must be positive. 

1297 freundlich_k : float, optional 

1298 Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive. 

1299 Used if retardation_factor is None. 

1300 freundlich_n : float, optional 

1301 Freundlich exponent [-]. Must be positive and != 1. 

1302 Used if retardation_factor is None. 

1303 bulk_density : float, optional 

1304 Bulk density [kg/m³]. Must be positive. 

1305 Used if retardation_factor is None. 

1306 porosity : float, optional 

1307 Porosity [-]. Must be in (0, 1). 

1308 Used if retardation_factor is None. 

1309 retardation_factor : float, optional 

1310 Constant retardation factor [-]. If provided, uses linear retardation 

1311 instead of Freundlich sorption. Must be >= 1.0. 

1312 max_iterations : int, optional 

1313 Maximum number of events. Default 10000. 

1314 

1315 Returns 

1316 ------- 

1317 cout : numpy.ndarray 

1318 Bin-averaged concentrations averaged across all pore volumes. 

1319 

1320 structures : list of dict 

1321 List of detailed simulation structures, one for each pore volume, with keys: 

1322 

1323 - 'waves': List[Wave] - All wave objects created during simulation 

1324 - 'events': List[dict] - All events with times, types, and details 

1325 - 't_first_arrival': float - First arrival time (end of spin-up period) 

1326 - 'n_events': int - Total number of events 

1327 - 'n_shocks': int - Number of shocks created 

1328 - 'n_rarefactions': int - Number of rarefactions created 

1329 - 'n_characteristics': int - Number of characteristics created 

1330 - 'final_time': float - Final simulation time 

1331 - 'sorption': FreundlichSorption | ConstantRetardation - Sorption object 

1332 - 'tracker_state': FrontTrackerState - Complete simulation state 

1333 - 'aquifer_pore_volume': float - Pore volume for this simulation 

1334 

1335 Examples 

1336 -------- 

1337 :: 

1338 

1339 cout, structures = infiltration_to_extraction_front_tracking_detailed( 

1340 cin=cin, 

1341 flow=flow, 

1342 tedges=tedges, 

1343 cout_tedges=cout_tedges, 

1344 aquifer_pore_volumes=np.array([500.0]), 

1345 freundlich_k=0.01, 

1346 freundlich_n=2.0, 

1347 bulk_density=1500.0, 

1348 porosity=0.3, 

1349 ) 

1350 

1351 # Access spin-up period for first pore volume 

1352 print(f"First arrival: {structures[0]['t_first_arrival']:.2f} days") 

1353 

1354 # Analyze events for first pore volume 

1355 for event in structures[0]["events"]: 

1356 print(f"t={event['time']:.2f}: {event['type']}") 

1357 

1358 See Also 

1359 -------- 

1360 infiltration_to_extraction_front_tracking : Returns concentrations only (simpler interface) 

1361 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and front-tracking theory 

1362 :ref:`assumption-advection-dominated` : When diffusion/dispersion is negligible 

1363 """ 

1364 cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days = _validate_front_tracking_inputs( 

1365 cin=cin, 

1366 flow=flow, 

1367 tedges=tedges, 

1368 cout_tedges=cout_tedges, 

1369 aquifer_pore_volumes=aquifer_pore_volumes, 

1370 freundlich_k=freundlich_k, 

1371 freundlich_n=freundlich_n, 

1372 bulk_density=bulk_density, 

1373 porosity=porosity, 

1374 retardation_factor=retardation_factor, 

1375 ) 

1376 

1377 # Loop over each pore volume and compute concentration 

1378 cout_all = np.zeros((len(aquifer_pore_volumes), len(cout_tedges) - 1)) 

1379 structures = [] 

1380 

1381 for i, aquifer_pore_volume in enumerate(aquifer_pore_volumes): 

1382 # Create tracker and run simulation for this pore volume 

1383 tracker = FrontTracker( 

1384 cin=cin, 

1385 flow=flow, 

1386 tedges=tedges, 

1387 aquifer_pore_volume=aquifer_pore_volume, 

1388 sorption=sorption, 

1389 ) 

1390 

1391 tracker.run(max_iterations=max_iterations) 

1392 

1393 # Extract bin-averaged concentrations for this pore volume 

1394 cout_all[i, :] = compute_bin_averaged_concentration_exact( 

1395 t_edges=cout_tedges_days, 

1396 v_outlet=aquifer_pore_volume, 

1397 waves=tracker.state.waves, 

1398 sorption=sorption, 

1399 ) 

1400 

1401 # Build detailed structure dict for this pore volume 

1402 structure = { 

1403 "waves": tracker.state.waves, 

1404 "events": tracker.state.events, 

1405 "t_first_arrival": tracker.t_first_arrival, 

1406 "n_events": len(tracker.state.events), 

1407 "n_shocks": sum(1 for w in tracker.state.waves if isinstance(w, ShockWave)), 

1408 "n_rarefactions": sum(1 for w in tracker.state.waves if isinstance(w, RarefactionWave)), 

1409 "n_characteristics": sum(1 for w in tracker.state.waves if isinstance(w, CharacteristicWave)), 

1410 "final_time": tracker.state.t_current, 

1411 "sorption": sorption, 

1412 "tracker_state": tracker.state, 

1413 "aquifer_pore_volume": aquifer_pore_volume, 

1414 } 

1415 structures.append(structure) 

1416 

1417 # Return average concentrations and list of structures 

1418 return np.mean(cout_all, axis=0), structures