Coverage for src/gwtransport/diffusion.py: 0%

1r""" 

2Analytical solutions for 1D advection-dispersion transport. 

3 

4Water infiltrates and is transported in parallel along multiple aquifer pore volumes to 

5extraction. For each aquifer pore volume, transport is 1D advection with microdispersion, 

6molecular diffusion, and linear sorption; the spread across aquifer pore volumes provides 

7macrodispersion. Forward and backward modeling are supported. The flow is assumed orthogonal. 

8 

9The orthogonal-flow (Cartesian) geometry is what makes the Kreft-Zuber breakthrough the exact 

101D solution used below. 

11 

12Choosing among the three diffusion modules 

13------------------------------------------ 

14 

15This is the reference implementation: it evaluates the bin-averaged Kreft-Zuber flux 

16concentration by resolution-aware composite Gauss-Legendre quadrature (splitting at 

17flow-bin boundaries, with extra front-centred panels wherever a sharp breakthrough front 

18is otherwise under-resolved). 

19Prefer it only when the output grid is coarser than the flow detail -- it integrates the 

20full within-bin flow, which the closed-form :mod:`gwtransport.diffusion_fast` approximates as 

21constant per output bin. Otherwise that module computes the same physics to machine 

22precision for *every* parameter regime (including ``retardation_factor != 1`` with non-zero 

23molecular diffusivity, whose flux correction it also evaluates in closed form) and is 

24~80-90x faster (no quadrature, no residence-time inversion). Both modules accept 

25per-streamtube ``streamline_length`` / ``molecular_diffusivity`` / 

26``longitudinal_dispersivity`` arrays (heterogeneous flow paths -- partially-penetrating 

27wells, wedge-shaped capture zones). 

28 

29:mod:`gwtransport.diffusion_fast_fast` is the third rung and is approximate by design: it folds 

30molecular diffusion into an effective dispersivity ``alpha_eff = alpha_L + D_m * R * V_pore / 

31(L * q_mean)`` per streamtube, so the whole streamtube bundle collapses to one banded breakthrough 

32on the native cumulative-volume grid -- the fastest of the three. At constant flow it reproduces 

33:mod:`gwtransport.diffusion_fast` to ~1e-6 for smooth inputs and ~1e-4 for sharp ones, and it loses 

34accuracy in the molecular-diffusion-dominated corner (``alpha_L`` ~ 0) under strongly variable flow, 

35where the frozen record-mean flow leaves a commutator residual. This module is the ground truth both 

36fast modules are validated against. 

37 

38Reported outlet concentration: Kreft-Zuber (1978) flux concentration 

39--------------------------------------------------------------------- 

40 

41The outlet concentration reported by this module is the **flux concentration** 

42 

43 C_F(L, t) = C_R(L, t) - (D_s / v_s) * dC_R/dx \|_{x=L} 

44 

45with the solute-front (retarded-frame) velocity v_s = Q L / (R V_pore) and the 

46dispersion D_s = D_m + alpha_L * v_s, so the flux coefficient is 

47D_s / v_s = D_m / v_s + alpha_L = R D_m / v_fluid + alpha_L (with the fluid 

48velocity v_fluid = Q L / V_pore). The resident profile C_R solves the retarded 

49ADE with advection v_s and dispersion D_s, so its flux-vs-resident correction 

50must use v_s — not v_fluid; pairing v_s with the moving-frame variance below is 

51what conserves mass for R > 1 with D_m > 0. 

52 

53— the solute mass flux at the outlet divided by the volumetric fluid flux. This 

54is what is measured when sampling the extracted fluid. The resident 

55concentration ``C_R`` is Bear (1972) eq. 10.6.4, the variable-flow moving-frame 

56Ogata-Banks solution 

57 

58 C_R(L, V; t_j) = 0.5 * erfc((L - xi_j(V)) / (2 * sqrt(D_t(V)))) 

59 

60with the dispersion variance accumulated in the moving (Lagrangian) frame: 

61 

62 D_t(V) = sigma^2(V) / 2 = D_m * tau(V) + alpha_L * xi(V) 

63 

64where: 

65 

66- D_m is the effective molecular (or thermal) diffusivity [m²/day] 

67- alpha_L is the longitudinal dispersivity [m] 

68- tau(V) is the elapsed time since infiltration [day], with V the cumulative 

69 extracted volume 

70- xi(V) = L (V - V_j) / (R V_pore) is the distance the parcel has actually 

71 travelled [m] 

72 

73The K-Z flux-correction term is what makes the column-sum invariant 

74``integral Q c_out dt = integral Q c_in dt`` hold under arbitrary variable Q. 

75Without it, the leading-order C_R loses O(1/Pe) per column under variable Q + 

76pure D_m. 

77 

78Implementation: the bin-averaged C_F is computed by resolution-aware composite 

79Gauss-Legendre quadrature in volume space, split at flow-bin boundaries so each 

80sub-interval sees a linear t(V). Within a sub-interval the erf-like front has 

81width ``sqrt(4*D_t)`` (in volume units); for near-zero dispersivity this can be 

82orders of magnitude below the flow-bin width, so a single fixed-order rule 

83cannot resolve it. Sub-intervals whose front is under-resolved are therefore 

84tiled with front-centred panels (fine near the front, flat tails outside), 

85which restores the column-mass invariant to ~1e-11 for every dispersion regime; 

86smooth/already-resolved sub-intervals keep the plain single 16-point rule. The 

87variance is evaluated at each quadrature node from the parcel's own tau and xi 

88histories — never capped at the residence time. The K-Z identity requires 

89Bear's formula to satisfy the variable-coefficient ADE exactly, which holds only 

90when D_t is allowed to keep growing past breakthrough. 

91 

92Macrodispersion vs microdispersion 

93---------------------------------- 

94 

95This module adds microdispersion (alpha_L) and molecular diffusion (D_m) on top of 

96macrodispersion captured by the pore volume distribution (APVD). Both represent velocity 

97heterogeneity at different scales. Microdispersion is an aquifer property; macrodispersion 

98depends additionally on hydrological boundary conditions. See :ref:`concept-dispersion-scales` 

99for guidance on when to use each approach and how to avoid double-counting spreading effects. 

100 

101Streamtube assumption (no cross-sectional area parameter) 

102--------------------------------------------------------- 

103 

104Each entry in ``aquifer_pore_volumes`` is treated as an independent 1D streamtube. There is 

105no cross-sectional area parameter: the variance budget uses ``2 D_m tau`` (molecular 

106diffusion in time) and ``2 alpha_L xi`` (microdispersion in travelled distance), with 

107the streamline length ``L`` and the pore volume ``V_pore`` together fixing the implicit 

108streamtube cross-section ``A = V_pore / L``. Callers who need distributed-area effects must 

109provide multiple streamtubes (via ``aquifer_pore_volumes`` or the gamma-parameterised 

110wrappers). 

111 

112Available functions: 

113 

114- :func:`infiltration_to_extraction` - Forward transport from an explicit ``aquifer_pore_volumes`` 

115 distribution: returns the bin-averaged Kreft-Zuber flux concentration on ``cout_tedges``, integrated 

116 over the full ``tedges``-resolution flow within each output bin and averaged with equal weight over 

117 the streamtubes. ``streamline_length``, ``molecular_diffusivity`` and ``longitudinal_dispersivity`` 

118 are either scalars shared by all streamtubes or one value per pore volume. Output bins that no 

119 infiltration has reached, or whose look-back leaves the record, are NaN. 

120 

121- :func:`extraction_to_infiltration` - Reverse direction: builds the same forward coefficient matrix 

122 and solves ``W @ cin = cout`` by Tikhonov regularization, returning the bin-averaged infiltration 

123 concentration on ``tedges``. NaN entries in ``cout`` mark measurement gaps; their rows are dropped 

124 from the solve and cin bins that nothing else constrains are returned as NaN. 

125 

126- :func:`gamma_infiltration_to_extraction` - :func:`infiltration_to_extraction` with the pore volume 

127 distribution given as a (shifted) gamma -- either (mean, std) or (alpha, beta), plus ``loc`` -- 

128 discretized into ``n_bins`` equal-probability streamtubes that share one ``streamline_length`` and 

129 one pair of dispersion parameters. 

130 

131- :func:`gamma_extraction_to_infiltration` - :func:`extraction_to_infiltration` with the same gamma 

132 parameterization of the pore volume distribution: reconstructs ``cin`` on ``tedges`` from ``cout``. 

133 

134References 

135---------- 

136Bear, J. (1972). Dynamics of Fluids in Porous Media. American Elsevier 

137Publishing Company. Equation 10.6.4 (variable-flow Ogata-Banks form). Provides 

138the resident concentration ``C_R``. 

139 

140Kreft, A., & Zuber, A. (1978). On the physical meaning of the dispersion 

141equation and its solutions for different initial and boundary conditions. 

142Chemical Engineering Science, 33(11), 1471-1480. Eq. 2 gives the resident-to- 

143flux concentration transformation; Eq. 1 is the mass-balance identity that 

144makes the column-sum invariant exact. 

145""" 

146 

147import numpy as np 

148import numpy.typing as npt 

149import pandas as pd 

150from scipy import special 

151 

152from gwtransport import gamma 

153from gwtransport._diffusion_shared import _broadcast_to_pore_volumes, _extend_tedges_flag 

154from gwtransport._time import dt_to_days, tedges_to_days 

155from gwtransport._validation import ( 

156 _validate_no_nan, 

157 _validate_non_negative_array, 

158 _validate_positive_array, 

159 _validate_retardation_factor, 

160 _validate_scalar_or_matching_length, 

161 _validate_tedges_parity, 

162) 

163from gwtransport.residence_time import fraction_explained_full 

164from gwtransport.utils import cumulative_flow_volume, solve_inverse_transport 

165 

166# Numerical tolerance for coefficient sum to determine valid output bins 

167_EPSILON_COEFF_SUM = 1e-10 

168 

169# Gauss-Legendre quadrature nodes and weights for volume-space integration 

170_GL_NODES, _GL_WEIGHTS = np.polynomial.legendre.leggauss(16) 

171 

172# Resolution-aware composite-quadrature template for the erf-like breakthrough front. 

173# When the front width sqrt(4*D_t) (in volume units) is much smaller than a 

174# (cell, flow-bin) sub-interval, plain 16-point GL cannot resolve it. Cells flagged 

175# as under-resolved get panels placed at ``V_front + front_width * _FRONT_OFFSETS`` 

176# (clipped to the sub-interval); the two outer panels then cover the flat erf tails, 

177# where 16-point GL is already exact. Panels near the front are ~1 front-width wide, 

178# spanning +-6 front-widths (erf and the flux-correction Gaussian are flat to ~1e-16 

179# beyond that). A cell is refined only when its front lies inside the sub-interval and 

180# the sub-interval is wider than _REFINE_RATIO front-widths, so adaptivity triggers 

181# only near sharp fronts (smooth regimes keep the plain single-panel cost and answer). 

182_FRONT_REACH = 6.0 

183_FRONT_OFFSETS = np.arange(-_FRONT_REACH, _FRONT_REACH + 0.5, 1.0) 

184_REFINE_RATIO = 4.0 

185 

186 

187def _cfrac_mean_volume( 

188 *, 

189 step_widths: npt.NDArray[np.floating], 

190 cumulative_volume_at_cout_tedges: npt.NDArray[np.floating], 

191 cumulative_volume_at_cin_tedges: npt.NDArray[np.floating], 

192 tedges_days: npt.NDArray[np.floating], 

193 molecular_diffusivity: float, 

194 longitudinal_dispersivity: float, 

195 r_vpv: float, 

196 streamline_len: float, 

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

198 r"""Compute bin-averaged flux concentration at the outlet for each cell. 

199 

200 For each cell (cout-bin *i*, cin-edge *j*), computes the flow-weighted 

201 average of the Kreft-Zuber (1978) **flux concentration** at the outlet: 

202 

203 .. math:: 

204 

205 \text{frac}_{i,j} = \frac{1}{\Delta V_i} 

206 \int_{V_i}^{V_{i+1}} C_F\!\left(L,\,V;\,t_j\right) dV 

207 

208 with :math:`C_F = C_R - (D_s / v_s) \, \partial_x C_R\big|_{x=L}`, Bear's 

209 resident concentration :math:`C_R`, the moving-frame variance 

210 :math:`D_t = D_m \tau + \alpha_L \xi` and the solute-front flux coefficient 

211 :math:`D_s/v_s`, all as derived in the module docstring. The correction term 

212 added to :math:`C_R` at each quadrature node is 

213 

214 .. math:: 

215 

216 \frac{D_s}{v_s(t(V))} \cdot 

217 \frac{1}{\sqrt{4\pi\,D_t(V)}}\, 

218 \exp\!\left( -\frac{(L - \xi_j(V))^2}{4\,D_t(V)} \right). 

219 

220 Implementation: resolution-aware composite Gauss-Legendre quadrature in 

221 volume space, split at flow-bin boundaries so that within each sub-interval 

222 :math:`t(V)` is linear. The erf-like front has width :math:`\sqrt{4 D_t}` 

223 (in volume units); a sub-interval whose front is under-resolved by a single 

224 16-point rule (front width far below the sub-interval width) is tiled with 

225 front-centred panels (see ``_FRONT_OFFSETS``), while smooth/already-resolved 

226 sub-intervals keep the plain single 16-point rule (bit-identical to it *per 

227 sub-interval*; a cell split at flow-bin boundaries still integrates more 

228 accurately than one un-split rule over the whole cell). The moving-frame 

229 variance keeps growing past breakthrough and is never capped: the K-Z 

230 identity holds only while Bear's formula satisfies the variable-coefficient 

231 ADE exactly. 

232 

233 Parameters 

234 ---------- 

235 step_widths : ndarray, shape (n_cout_edges, n_cin_edges) 

236 x-position ``x(V_cout, V_cin) = (V_cout - V_cin - r_vpv) * L / r_vpv`` 

237 at each (cout-edge, cin-edge). NaN for inactive cells. Equals 

238 :math:`\xi - L`. 

239 cumulative_volume_at_cout_tedges : ndarray, shape (n_cout_edges,) 

240 Cumulative extracted volume at each cout time edge [m³]. 

241 cumulative_volume_at_cin_tedges : ndarray, shape (n_cin_edges,) 

242 Cumulative volume at each cin (flow) time edge [m³]. 

243 tedges_days : ndarray, shape (n_cin_edges,) 

244 Flow time edges in days. 

245 molecular_diffusivity : float 

246 Effective (retarded-frame) molecular diffusivity D_m [m²/day]. 

247 Contributes ``D_m * tau`` to the dispersion product ``D_t``. 

248 longitudinal_dispersivity : float 

249 Longitudinal dispersivity alpha_L [m]. Contributes ``alpha_L * xi`` 

250 to the dispersion product ``D_t``. 

251 r_vpv : float 

252 Retardation factor times pore volume = R * V_pore [m³]. 

253 streamline_len : float 

254 Streamline length L [m]. 

255 

256 Returns 

257 ------- 

258 ndarray, shape (n_cout_bins, n_cin_edges) 

259 Bin-averaged flux concentration for each cell. NaN for inactive cells. 

260 """ 

261 n_cout_edges, n_cin_edges = step_widths.shape 

262 n_cout_bins = n_cout_edges - 1 

263 

264 x_lo = step_widths[:-1] 

265 x_hi = step_widths[1:] 

266 dx = x_hi - x_lo 

267 

268 v_lo_arr = cumulative_volume_at_cout_tedges[:-1] 

269 v_hi_arr = cumulative_volume_at_cout_tedges[1:] 

270 

271 is_valid = ~np.isnan(x_lo) & ~np.isnan(x_hi) 

272 

273 frac = np.full((n_cout_bins, n_cin_edges), np.nan) 

274 

275 # --- No dispersion: C_F = C_R = step function (no dispersive flux) --- 

276 if molecular_diffusivity == 0.0 and longitudinal_dispersivity == 0.0: 

277 with np.errstate(divide="ignore", invalid="ignore"): 

278 cr_no_disp = 0.5 + 0.5 * (np.abs(x_hi) - np.abs(x_lo)) / dx 

279 cr_no_disp = np.where(dx == 0.0, 0.5 + 0.5 * np.sign(x_lo), cr_no_disp) 

280 return np.where(is_valid, cr_no_disp, frac) 

281 

282 # --- Pre-compute solute-front velocity and K-Z coefficient (D_s/v_s) per flow bin --- 

283 dv_per_bin = np.diff(cumulative_volume_at_cin_tedges) 

284 dt_per_bin = np.diff(tedges_days) 

285 with np.errstate(divide="ignore", invalid="ignore"): 

286 q_per_bin = np.where(dt_per_bin > 0, dv_per_bin / dt_per_bin, 0.0) 

287 # Solute-front velocity v_s = Q L / (R V_pore), the advection speed of the retarded ADE 

288 # that C_R solves -- not the fluid velocity Q L / V_pore. 

289 v_per_bin = q_per_bin * streamline_len / r_vpv 

290 # (D_s/v_s) = D_m/v_s + alpha_L. At v_s=0 the bin has dV=0 and is skipped below; the 

291 # surrounding errstate suppresses the divide warning for those lanes. 

292 dl_over_v_per_bin = np.where( 

293 v_per_bin > 0, 

294 molecular_diffusivity / v_per_bin + longitudinal_dispersivity, 

295 0.0, 

296 ) 

297 

298 # --- Resolution-aware composite Gauss-Legendre quadrature, split by flow bins --- 

299 # The integration window per (cell, flow-bin) is the intersection of 

300 # [V_lo, V_hi] (cell), [ve_lo, ve_hi] (flow bin), and [V_j, infty) (parcel 

301 # entered). Within each sub-interval, t(V) is linear. Where the sub-interval is 

302 # much wider than the front width sqrt(4*D_t), the erf-like front is under-resolved 

303 # by a single 16-point GL rule; such sub-intervals are tiled with front-centred 

304 # panels (see _FRONT_OFFSETS). Smooth sub-intervals keep the single panel and are 

305 # bit-identical to the plain 16-point rule. 

306 idx_i, idx_j = np.nonzero(is_valid) 

307 if len(idx_i) == 0: 

308 return frac 

309 

310 v_lo_cells = v_lo_arr[idx_i] 

311 v_hi_cells = v_hi_arr[idx_i] 

312 v_cin_cells = cumulative_volume_at_cin_tedges[idx_j] 

313 t_j_cells = tedges_days[idx_j] 

314 total_dv = v_hi_cells - v_lo_cells 

315 valid_cells = total_dv > 0 

316 

317 # Post-injection lower bound, loop-invariant: max(V_lo, V_cin). 

318 v_lo_or_cin = np.maximum(v_lo_cells, v_cin_cells) 

319 

320 integral_cf = np.zeros(len(idx_i)) 

321 

322 vol_edges = cumulative_volume_at_cin_tedges 

323 for k in range(len(vol_edges) - 1): 

324 ve_lo, ve_hi = vol_edges[k], vol_edges[k + 1] 

325 if v_per_bin[k] <= 0.0: 

326 continue 

327 dl_over_v_k = dl_over_v_per_bin[k] 

328 dt_sub_bin = tedges_days[k + 1] - tedges_days[k] 

329 dv_sub_edge = ve_hi - ve_lo 

330 

331 # Intersection of cell, flow-bin, and post-injection range 

332 sub_lo = np.maximum(v_lo_or_cin, ve_lo) 

333 sub_hi = np.minimum(v_hi_cells, ve_hi) 

334 overlap = (sub_hi > sub_lo) & valid_cells 

335 if not np.any(overlap): 

336 continue 

337 

338 lo = sub_lo[overlap] 

339 hi = sub_hi[overlap] 

340 vcin = v_cin_cells[overlap] 

341 tj = t_j_cells[overlap] 

342 

343 # Front centre (x = 0 => xi = L => V = V_cin + r_vpv) and its width in volume. 

344 # front_width = sqrt(4*D_t_front) * r_vpv / L, with D_t_front = D_m*tau_front + 

345 # alpha_L*L (xi = L at the front). t(V) is linear within flow bin k. 

346 v_front = vcin + r_vpv 

347 t_front = tedges_days[k] + (v_front - ve_lo) * (dt_sub_bin / dv_sub_edge) 

348 tau_front = np.maximum(t_front - tj, 0.0) 

349 dt_front = molecular_diffusivity * tau_front + longitudinal_dispersivity * streamline_len 

350 front_width = np.sqrt(4.0 * dt_front) * r_vpv / streamline_len 

351 

352 # Refine only where a sharp front region intersects this sub-interval and is 

353 # under-resolved by a single 16-point rule; elsewhere the integrand is flat 

354 # or already resolved, so one panel is exact. The front-region test (not just 

355 # "centre in this bin") also catches a front whose tail spills across a 

356 # flow-bin boundary. A spuriously large extrapolated front_width far from the 

357 # front fails ``underresolved`` and so cannot trigger refinement. 

358 front_hits = (v_front + _FRONT_REACH * front_width > lo) & (v_front - _FRONT_REACH * front_width < hi) 

359 underresolved = (hi - lo) > _REFINE_RATIO * front_width 

360 if np.any(front_hits & underresolved): 

361 inner = np.clip( 

362 v_front[:, np.newaxis] + front_width[:, np.newaxis] * _FRONT_OFFSETS[np.newaxis, :], 

363 lo[:, np.newaxis], 

364 hi[:, np.newaxis], 

365 ) 

366 edges = np.concatenate([lo[:, np.newaxis], inner, hi[:, np.newaxis]], axis=1) 

367 else: 

368 edges = np.stack([lo, hi], axis=1) 

369 

370 p_lo = edges[:, :-1] 

371 p_hi = edges[:, 1:] 

372 p_mid = 0.5 * (p_lo + p_hi) 

373 p_half = 0.5 * (p_hi - p_lo) 

374 

375 # GL nodes over every panel: shape (n_cell, n_panel, n_gl) 

376 v_nodes = p_mid[:, :, np.newaxis] + p_half[:, :, np.newaxis] * _GL_NODES[np.newaxis, np.newaxis, :] 

377 

378 # Geometry: x = xi - L = (V - V_j - r_vpv) * L / r_vpv (parcel position 

379 # relative to outlet); xi = parcel travel distance. 

380 x_nodes = (v_nodes - vcin[:, np.newaxis, np.newaxis] - r_vpv) * streamline_len / r_vpv 

381 xi_nodes = x_nodes + streamline_len 

382 

383 t_nodes = tedges_days[k] + (v_nodes - ve_lo) * (dt_sub_bin / dv_sub_edge) 

384 # tau >= 0 by construction (lo >= v_cin); clip for safety. 

385 tau_nodes = np.maximum(t_nodes - tj[:, np.newaxis, np.newaxis], 0.0) 

386 

387 # Bear's variance accumulator (sigma^2/2) — NO capping at RT/L 

388 dt_nodes = molecular_diffusivity * tau_nodes + longitudinal_dispersivity * xi_nodes 

389 

390 with np.errstate(divide="ignore", invalid="ignore"): 

391 arg = x_nodes / (2.0 * np.sqrt(dt_nodes)) 

392 # C_R = 0.5 * (1 + erf(arg)) = 0.5 * erfc((L-xi)/(2*sqrt(Dt))) 

393 erf_vals = np.where(np.isfinite(arg), special.erf(arg), np.sign(x_nodes)) 

394 cr_vals = 0.5 * (1.0 + erf_vals) 

395 

396 # K-Z flux correction: FC = (D_s/v_s) * (1/sqrt(4 pi D_t)) * exp(-arg^2) 

397 with np.errstate(divide="ignore", invalid="ignore"): 

398 gauss_vals = np.where( 

399 dt_nodes > 0.0, 

400 np.exp(-(arg**2)) / np.sqrt(4.0 * np.pi * dt_nodes), 

401 0.0, 

402 ) 

403 cf_vals = cr_vals + dl_over_v_k * gauss_vals 

404 

405 # Integrate: GL-weight over nodes, then sum panel contributions per cell. 

406 # The weight contraction is done in 2D so a single-panel (non-refined) 

407 # sub-interval is bit-identical to a plain 16-point rule. 

408 cf_weighted = (cf_vals.reshape(-1, cf_vals.shape[-1]) @ _GL_WEIGHTS).reshape(cf_vals.shape[:-1]) 

409 integral_cf[overlap] += (p_half * cf_weighted).sum(axis=1) 

410 

411 with np.errstate(divide="ignore", invalid="ignore"): 

412 frac_cells = np.where(valid_cells, integral_cf / total_dv, np.nan) 

413 frac[idx_i, idx_j] = frac_cells 

414 

415 return frac 

416 

417 

418def _validate_diffusion_inputs( 

419 *, 

420 tedges: pd.DatetimeIndex, 

421 flow: npt.NDArray[np.floating], 

422 aquifer_pore_volumes: npt.NDArray[np.floating], 

423 streamline_length: npt.NDArray[np.floating], 

424 molecular_diffusivity: npt.NDArray[np.floating], 

425 longitudinal_dispersivity: npt.NDArray[np.floating], 

426 retardation_factor: float, 

427 cin_values: npt.NDArray[np.floating] | None = None, 

428 cout_values: npt.NDArray[np.floating] | None = None, 

429 cout_tedges: pd.DatetimeIndex | None = None, 

430) -> None: 

431 """Validate inputs common to diffusion forward / reverse entry points. 

432 

433 The caller supplies the kwargs of its own direction; the parity checks that apply to 

434 that direction then run: 

435 

436 - ``cin_values`` (forward): ``tedges`` parities cin and flow. 

437 - ``cout_values`` + ``cout_tedges`` (reverse): ``tedges`` parities flow, ``cout_tedges`` 

438 parities cout. 

439 

440 The physical-parameter checks below run either way. 

441 

442 Raises 

443 ------ 

444 ValueError 

445 If any check fails. The message identifies which invariant was violated. 

446 """ 

447 n_pore_volumes = len(aquifer_pore_volumes) 

448 

449 if cin_values is not None: 

450 _validate_tedges_parity(tedges, cin_values, tedges_name="tedges", values_name="cin") 

451 _validate_tedges_parity(tedges, flow, tedges_name="tedges", values_name="flow") 

452 elif cout_values is not None and cout_tedges is not None: 

453 _validate_tedges_parity(tedges, flow, tedges_name="tedges", values_name="flow") 

454 _validate_tedges_parity(cout_tedges, cout_values, tedges_name="cout_tedges", values_name="cout") 

455 if len(aquifer_pore_volumes) != len(streamline_length): 

456 msg = "aquifer_pore_volumes and streamline_length must have the same length" 

457 raise ValueError(msg) 

458 _validate_scalar_or_matching_length( 

459 molecular_diffusivity, 

460 name="molecular_diffusivity", 

461 expected_len=n_pore_volumes, 

462 ref_name="aquifer_pore_volumes", 

463 ) 

464 _validate_scalar_or_matching_length( 

465 longitudinal_dispersivity, 

466 name="longitudinal_dispersivity", 

467 expected_len=n_pore_volumes, 

468 ref_name="aquifer_pore_volumes", 

469 ) 

470 _validate_non_negative_array(molecular_diffusivity, name="molecular_diffusivity") 

471 _validate_non_negative_array(longitudinal_dispersivity, name="longitudinal_dispersivity") 

472 # Reverse cout may contain NaN (measurement gaps); gapped rows are excluded from the solve. 

473 if cin_values is not None: 

474 _validate_no_nan(cin_values, name="cin") 

475 _validate_no_nan(flow, name="flow") 

476 _validate_non_negative_array(flow, name="flow", message="flow must be non-negative (negative flow not supported)") 

477 _validate_positive_array(aquifer_pore_volumes, name="aquifer_pore_volumes") 

478 _validate_positive_array(streamline_length, name="streamline_length") 

479 _validate_retardation_factor(retardation_factor) 

480 

481 

482def _prepare_diffusion_arrays( 

483 *, 

484 flow: npt.ArrayLike, 

485 aquifer_pore_volumes: npt.ArrayLike, 

486 streamline_length: npt.ArrayLike, 

487 molecular_diffusivity: npt.ArrayLike, 

488 longitudinal_dispersivity: npt.ArrayLike, 

489) -> tuple[ 

490 npt.NDArray[np.floating], 

491 npt.NDArray[np.floating], 

492 npt.NDArray[np.floating], 

493 npt.NDArray[np.floating], 

494 npt.NDArray[np.floating], 

495]: 

496 """Coerce flow / geometry / dispersion inputs to broadcasted float arrays. 

497 

498 Each per-streamtube parameter (``streamline_length``, ``molecular_diffusivity``, 

499 ``longitudinal_dispersivity``) may be passed as a scalar; it is broadcast to one 

500 value per pore volume. The returned arrays are read-only views when broadcast (none 

501 is mutated downstream). 

502 

503 Returns 

504 ------- 

505 tuple of ndarray 

506 ``(flow, aquifer_pore_volumes, streamline_length, molecular_diffusivity, 

507 longitudinal_dispersivity)`` as float arrays. 

508 """ 

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

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

511 n_pore_volumes = len(aquifer_pore_volumes) 

512 return ( 

513 flow, 

514 aquifer_pore_volumes, 

515 _broadcast_to_pore_volumes(streamline_length, n_pore_volumes), 

516 _broadcast_to_pore_volumes(molecular_diffusivity, n_pore_volumes), 

517 _broadcast_to_pore_volumes(longitudinal_dispersivity, n_pore_volumes), 

518 ) 

519 

520 

521def _infiltration_to_extraction_coeff_matrix( 

522 *, 

523 flow: npt.NDArray[np.floating], 

524 tedges: pd.DatetimeIndex, 

525 cout_tedges: pd.DatetimeIndex, 

526 aquifer_pore_volumes: npt.NDArray[np.floating], 

527 streamline_length: npt.NDArray[np.floating], 

528 molecular_diffusivity: npt.NDArray[np.floating], 

529 longitudinal_dispersivity: npt.NDArray[np.floating], 

530 retardation_factor: float, 

531 extend_tedges: bool = True, 

532) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.bool_]]: 

533 """Build the forward coefficient matrix for diffusion transport. 

534 

535 Constructs the matrix W such that ``cout = W @ cin``, accounting for 

536 advection, microdispersion, and molecular diffusion. NaN entries in the raw coefficient 

537 matrix are replaced with zero. 

538 

539 Parameters 

540 ---------- 

541 flow : ndarray 

542 Flow rate of water [m³/day]. Already validated. 

543 tedges : DatetimeIndex 

544 Cin/flow time edges (not yet extended for spin-up). 

545 cout_tedges : DatetimeIndex 

546 Cout time edges. 

547 aquifer_pore_volumes : ndarray 

548 Pore volumes [m³]. Already validated. 

549 streamline_length : ndarray 

550 Travel distances [m]. Already validated. 

551 molecular_diffusivity : ndarray 

552 Effective molecular diffusivities [m²/day]. Already broadcasted. 

553 See :func:`infiltration_to_extraction` for physical interpretation. 

554 longitudinal_dispersivity : ndarray 

555 Longitudinal dispersivities [m]. Already broadcasted. 

556 retardation_factor : float 

557 Retardation factor. 

558 

559 Returns 

560 ------- 

561 coeff_matrix : ndarray 

562 Filled coefficient matrix of shape (n_cout, n_cin). NaN replaced 

563 with zero. 

564 valid_cout_bins : ndarray 

565 Boolean mask of shape (n_cout,) indicating valid output bins. 

566 """ 

567 if extend_tedges: 

568 # Extend tedges by 100 years on each side to provide warm-start cin 

569 # and post-data flow for cout bins near the boundaries. Equivalent to 

570 # the public ``spinup="constant"`` policy in other modules. 

571 tedges = pd.DatetimeIndex([ 

572 tedges[0] - pd.Timedelta("36500D"), 

573 *list(tedges[1:-1]), 

574 tedges[-1] + pd.Timedelta("36500D"), 

575 ]) 

576 

577 # Compute the cumulative flow at tedges 

578 cumulative_volume_at_cin_tedges = cumulative_flow_volume(flow, dt_to_days(tedges)) # m³ 

579 

580 # Compute the cumulative flow at cout_tedges. Both edge arrays are first reduced to a shared 

581 # day axis: np.interp coerces each datetime64 operand to int64 in its own resolution, so a 

582 # cout_tedges / tedges unit mismatch (e.g. ns vs us) would send every query out of range and 

583 # silently return all-NaN. 

584 tedges_days_arr = tedges_to_days(tedges) 

585 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0]) 

586 cumulative_volume_at_cout_tedges = np.interp( 

587 cout_tedges_days, tedges_days_arr, cumulative_volume_at_cin_tedges 

588 ).astype(float) 

589 

590 # Output bin valid where every streamtube's advective look-back is in-record across the whole 

591 # bin -- i.e. advective coverage == 1 for all pore volumes (NaN outside the record -> invalid). 

592 # This is the advective validity gate only; the dispersive informedness is the captured kernel 

593 # mass (total_coeff) applied downstream. 

594 valid_cout_bins = np.all( 

595 fraction_explained_full( 

596 flow=flow, 

597 tedges=tedges, 

598 cout_tedges=cout_tedges, 

599 aquifer_pore_volumes=aquifer_pore_volumes, 

600 retardation_factor=retardation_factor, 

601 direction="extraction_to_infiltration", 

602 ) 

603 >= 1.0, 

604 axis=0, 

605 ) 

606 

607 # Initialize coefficient matrix accumulator 

608 n_cout_bins = len(cout_tedges) - 1 

609 n_cin_bins = len(flow) 

610 accumulated_coeff = np.zeros((n_cout_bins, n_cin_bins)) 

611 

612 # Loop over each pore volume 

613 for i_pv in range(len(aquifer_pore_volumes)): 

614 r_vpv = retardation_factor * aquifer_pore_volumes[i_pv] 

615 

616 delta_volume = cumulative_volume_at_cout_tedges[:, None] - cumulative_volume_at_cin_tedges[None, :] - r_vpv 

617 

618 step_widths = delta_volume / r_vpv * streamline_length[i_pv] 

619 

620 frac = _cfrac_mean_volume( 

621 step_widths=step_widths, 

622 cumulative_volume_at_cout_tedges=cumulative_volume_at_cout_tedges, 

623 cumulative_volume_at_cin_tedges=cumulative_volume_at_cin_tedges, 

624 tedges_days=tedges_days_arr, 

625 molecular_diffusivity=float(molecular_diffusivity[i_pv]), 

626 longitudinal_dispersivity=float(longitudinal_dispersivity[i_pv]), 

627 r_vpv=r_vpv, 

628 streamline_len=streamline_length[i_pv], 

629 ) 

630 

631 accumulated_coeff += frac[:, :-1] - frac[:, 1:] 

632 

633 coeff_matrix_filled = np.nan_to_num(accumulated_coeff / len(aquifer_pore_volumes), nan=0.0) 

634 

635 return coeff_matrix_filled, valid_cout_bins 

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 streamline_length: npt.ArrayLike, 

646 molecular_diffusivity: npt.ArrayLike, 

647 longitudinal_dispersivity: npt.ArrayLike, 

648 retardation_factor: float = 1.0, 

649 spinup: str | None = "constant", 

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

651 """ 

652 Compute extracted concentration with advection, microdispersion, and molecular diffusion. 

653 

654 This function models 1D solute transport through an aquifer system along orthogonal 

655 (Cartesian) flow paths. Each aquifer pore volume is an independent streamline carrying 

656 advection with microdispersion (alpha_L) and molecular diffusion (D_m); the spread across 

657 the pore volume distribution provides macrodispersion. Linear sorption enters via the 

658 retardation factor. 

659 

660 The physical model assumes: 

661 

662 1. Water infiltrates with concentration cin at time t_in 

663 2. Water travels distance L through aquifer with residence time tau = V_pore / Q 

664 3. During transport, microdispersion and molecular diffusion spread each streamline, 

665 while the spread across pore volumes provides macrodispersion 

666 4. At extraction, the concentration is a dispersed breakthrough curve 

667 

668 The reported extracted concentration is the Kreft-Zuber (1978) **flux 

669 concentration** at the outlet -- the solute mass flux divided by the 

670 volumetric fluid flux, which is what is measured when sampling the 

671 outflowing fluid. Microdispersion and molecular diffusion enter as the 

672 moving-frame variance ``sigma^2(V) = 2*D_m*tau(V) + 2*alpha_L*xi(V)``, with 

673 ``tau(V)`` the elapsed time since infiltration and ``xi(V)`` the distance the 

674 parcel has actually travelled. See the module docstring for the derivation 

675 and for why the flux form is what makes the column-sum invariant 

676 ``integral Q c_out dt = integral Q c_in dt`` hold exactly under variable flow. 

677 

678 Parameters 

679 ---------- 

680 cin : array-like 

681 Concentration of the compound in infiltrating water [concentration units]. 

682 Length must match the number of time bins defined by tedges. The model assumes 

683 this value is constant over each interval ``[tedges[i], tedges[i+1])``. 

684 flow : array-like 

685 Flow rate of water in the aquifer [m³/day]. 

686 Length must match cin and the number of time bins defined by tedges. The model 

687 assumes this value is constant over each interval ``[tedges[i], tedges[i+1])``. 

688 tedges : pandas.DatetimeIndex 

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

690 len(cin) + 1. 

691 cout_tedges : pandas.DatetimeIndex 

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

693 The output concentration is averaged over each bin. 

694 aquifer_pore_volumes : array-like 

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

696 of flow paths. Each pore volume determines the residence time for 

697 that flow path: tau = V_pore / Q. 

698 streamline_length : array-like 

699 Array of travel distances [m] corresponding to each pore volume. 

700 Must have the same length as aquifer_pore_volumes. 

701 molecular_diffusivity : float or array-like 

702 Effective (retarded-frame) molecular diffusivity [m²/day]. Can be a 

703 scalar (same for all pore volumes) or an array with the same length as 

704 aquifer_pore_volumes. Must be non-negative. For solute transport, this is 

705 the molecular diffusion coefficient D_m [m²/day] — typically ~1e-5 m²/day, 

706 negligible compared to microdispersion. For heat transport, pass the 

707 thermal diffusivity D_th = lambda / (rho*c)_eff [m²/day], typically 

708 0.01-0.1 m²/day. 

709 

710 Internally, this contributes ``2 * molecular_diffusivity * tau`` to the 

711 variance, where ``tau`` is the elapsed time in days (no extra factor of 

712 R). The retardation factor instead enters the flux coefficient 

713 ``D_s/v_s = R D_m / v_fluid + alpha_L`` through the solute-front velocity 

714 ``v_s = Q L / (R V_pore)``. For heat transport, the thermal diffusivity 

715 already represents the effective diffusivity D_eff in the porous matrix; 

716 for solutes the contribution is typically negligible. 

717 longitudinal_dispersivity : float or array-like 

718 Longitudinal dispersivity [m]. Can be a scalar (same for all pore 

719 volumes) or an array with the same length as aquifer_pore_volumes. 

720 Must be non-negative. Represents microdispersion from pore-scale velocity variations. 

721 Set to 0 for pure molecular diffusion. 

722 retardation_factor : float, optional 

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

724 Values > 1.0 indicate slower transport due to sorption. 

725 spinup : {'constant'} or None, optional 

726 Spin-up policy (default ``'constant'``). ``'constant'`` extends tedges by 

727 100 years on each side so that output bins near the boundary are always 

728 informed. ``None`` disables the extension; output bins without sufficient 

729 upstream data become NaN. Float fraction-threshold mode is not implemented 

730 and raises ``NotImplementedError``. 

731 

732 Returns 

733 ------- 

734 numpy.ndarray 

735 Bin-averaged concentration in the extracted water. Same units as cin. 

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

737 with no valid contributions from the infiltration data. 

738 

739 Raises 

740 ------ 

741 ValueError 

742 If input dimensions are inconsistent, if diffusivity is negative, 

743 or if aquifer_pore_volumes and streamline_length have different lengths. 

744 

745 See Also 

746 -------- 

747 extraction_to_infiltration : Inverse operation (deconvolution) 

748 gwtransport.advection.infiltration_to_extraction : Pure advection (no dispersion) 

749 gwtransport.diffusion_fast.infiltration_to_extraction : Fast closed-form equivalent 

750 :ref:`concept-dispersion-scales` : Macrodispersion vs microdispersion 

751 

752 Notes 

753 ----- 

754 The algorithm constructs a coefficient matrix W where cout = W @ cin: 

755 

756 1. For each pore volume, build a cell grid in cumulative volume space: 

757 

758 - cells span ``(V_cout[i], V_cout[i+1]) x V_cin[j]`` for each 

759 (cout-bin i, cin-edge j) 

760 - delta_volume = V_cout - V_cin - r_vpv encodes the parcel's offset 

761 from the outlet at each (cout-edge, cin-edge) 

762 

763 2. For each cell, compute the bin-averaged Kreft-Zuber flux concentration 

764 ``frac[i, j] = (1/dV_i) * integral C_F(L, V; t_j) dV`` by resolution-aware 

765 composite Gauss-Legendre quadrature in volume space, split at flow-bin 

766 boundaries so that ``t(V)`` is linear within each sub-interval. Where the 

767 erf-like front (width ``sqrt(4*D_t)`` in volume units) is under-resolved by 

768 a single 16-point rule -- as for near-zero dispersivity -- the sub-interval 

769 is tiled with front-centred panels; smooth sub-intervals keep the single 

770 rule. The moving-frame variance ``D_t = D_m*tau + alpha_L*xi`` is evaluated 

771 at each quadrature node (never capped at the residence time). 

772 

773 3. Coefficient for bin: ``coeff[i,j] = frac[i, j] - frac[i, j+1]``. This 

774 is the contribution of cin[j] to cout[i] in the W matrix. 

775 

776 4. Average coefficients across all pore volumes. 

777 

778 Examples 

779 -------- 

780 Basic usage with constant flow: 

781 

782 >>> import pandas as pd 

783 >>> import numpy as np 

784 >>> from gwtransport.diffusion import infiltration_to_extraction 

785 >>> 

786 >>> # Create time edges 

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

788 >>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D") 

789 >>> 

790 >>> # Input concentration (step function) and constant flow 

791 >>> cin = np.zeros(len(tedges) - 1) 

792 >>> cin[5:10] = 1.0 # Pulse of concentration 

793 >>> flow = np.ones(len(tedges) - 1) * 100.0 # 100 m³/day 

794 >>> 

795 >>> # Single pore volume of 500 m³, travel distance 100 m 

796 >>> aquifer_pore_volumes = np.array([500.0]) 

797 >>> streamline_length = np.array([100.0]) 

798 >>> 

799 >>> # Compute with dispersion (molecular diffusion + dispersivity) 

800 >>> # Scalar values broadcast to all pore volumes 

801 >>> cout = infiltration_to_extraction( 

802 ... cin=cin, 

803 ... flow=flow, 

804 ... tedges=tedges, 

805 ... cout_tedges=cout_tedges, 

806 ... aquifer_pore_volumes=aquifer_pore_volumes, 

807 ... streamline_length=streamline_length, 

808 ... molecular_diffusivity=1e-4, # m²/day, same for all pore volumes 

809 ... longitudinal_dispersivity=1.0, # m, same for all pore volumes 

810 ... ) 

811 

812 With multiple pore volumes (heterogeneous aquifer): 

813 

814 >>> # Distribution of pore volumes and corresponding travel distances 

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

816 >>> streamline_length = np.array([80.0, 100.0, 120.0]) 

817 >>> 

818 >>> # Scalar diffusion parameters broadcast to all pore volumes 

819 >>> cout = infiltration_to_extraction( 

820 ... cin=cin, 

821 ... flow=flow, 

822 ... tedges=tedges, 

823 ... cout_tedges=cout_tedges, 

824 ... aquifer_pore_volumes=aquifer_pore_volumes, 

825 ... streamline_length=streamline_length, 

826 ... molecular_diffusivity=1e-4, # m²/day 

827 ... longitudinal_dispersivity=1.0, # m 

828 ... ) 

829 """ 

830 cout_tedges = pd.DatetimeIndex(cout_tedges) 

831 tedges = pd.DatetimeIndex(tedges) 

832 

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

834 flow, aquifer_pore_volumes, streamline_length, molecular_diffusivity, longitudinal_dispersivity = ( 

835 _prepare_diffusion_arrays( 

836 flow=flow, 

837 aquifer_pore_volumes=aquifer_pore_volumes, 

838 streamline_length=streamline_length, 

839 molecular_diffusivity=molecular_diffusivity, 

840 longitudinal_dispersivity=longitudinal_dispersivity, 

841 ) 

842 ) 

843 

844 _validate_diffusion_inputs( 

845 tedges=tedges, 

846 flow=flow, 

847 aquifer_pore_volumes=aquifer_pore_volumes, 

848 streamline_length=streamline_length, 

849 molecular_diffusivity=molecular_diffusivity, 

850 longitudinal_dispersivity=longitudinal_dispersivity, 

851 retardation_factor=retardation_factor, 

852 cin_values=cin, 

853 ) 

854 

855 extend_tedges = _extend_tedges_flag(spinup) 

856 coeff_matrix, valid_cout_bins = _infiltration_to_extraction_coeff_matrix( 

857 flow=flow, 

858 tedges=tedges, 

859 cout_tedges=cout_tedges, 

860 aquifer_pore_volumes=aquifer_pore_volumes, 

861 streamline_length=streamline_length, 

862 molecular_diffusivity=molecular_diffusivity, 

863 longitudinal_dispersivity=longitudinal_dispersivity, 

864 retardation_factor=retardation_factor, 

865 extend_tedges=extend_tedges, 

866 ) 

867 

868 cout = coeff_matrix @ cin 

869 

870 # Output bins are invalid where the coefficient sum is near zero (no cin has broken 

871 # through yet) or the bin extends beyond the input data range (valid_cout_bins). 

872 total_coeff = np.sum(coeff_matrix, axis=1) 

873 no_valid_contribution = (total_coeff < _EPSILON_COEFF_SUM) | ~valid_cout_bins 

874 cout[no_valid_contribution] = np.nan 

875 

876 return cout 

877 

878 

879def extraction_to_infiltration( 

880 *, 

881 cout: npt.ArrayLike, 

882 flow: npt.ArrayLike, 

883 tedges: pd.DatetimeIndex, 

884 cout_tedges: pd.DatetimeIndex, 

885 aquifer_pore_volumes: npt.ArrayLike, 

886 streamline_length: npt.ArrayLike, 

887 molecular_diffusivity: npt.ArrayLike, 

888 longitudinal_dispersivity: npt.ArrayLike, 

889 retardation_factor: float = 1.0, 

890 regularization_strength: float = 1e-10, 

891 spinup: str | None = "constant", 

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

893 """ 

894 Compute infiltration concentration from extracted water (deconvolution with dispersion). 

895 

896 Inverts the forward transport model by building the forward coefficient 

897 matrix ``W_forward`` from :func:`infiltration_to_extraction` and solving 

898 ``W_forward @ cin = cout`` via Tikhonov regularization. Well-determined 

899 modes are dominated by the data; poorly-determined modes are pulled 

900 toward the physically motivated target (transpose-and-normalize of the 

901 forward matrix). 

902 

903 Parameters 

904 ---------- 

905 cout : array-like 

906 Concentration of the compound in extracted water [concentration units]. 

907 Length must match the number of time bins defined by cout_tedges. 

908 flow : array-like 

909 Flow rate of water in the aquifer [m³/day]. 

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

911 tedges : pandas.DatetimeIndex 

912 Time edges defining bins for cin (output) and flow data. 

913 Has length of len(flow) + 1. Output cin has length len(tedges) - 1. 

914 cout_tedges : pandas.DatetimeIndex 

915 Time edges for cout data bins. Has length of len(cout) + 1. 

916 Can have different time alignment and resolution than tedges. 

917 aquifer_pore_volumes : array-like 

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

919 of flow paths. Each pore volume determines the residence time for 

920 that flow path: tau = V_pore / Q. 

921 streamline_length : array-like 

922 Array of travel distances [m] corresponding to each pore volume. 

923 Must have the same length as aquifer_pore_volumes. 

924 molecular_diffusivity : float or array-like 

925 Effective molecular diffusivity [m²/day]. Can be a scalar (same for all 

926 pore volumes) or an array with the same length as aquifer_pore_volumes. 

927 Must be non-negative. See :func:`infiltration_to_extraction` for 

928 details on the physical interpretation and the interaction with 

929 retardation_factor. 

930 longitudinal_dispersivity : float or array-like 

931 Longitudinal dispersivity [m]. Can be a scalar (same for all pore 

932 volumes) or an array with the same length as aquifer_pore_volumes. 

933 Must be non-negative. 

934 retardation_factor : float, optional 

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

936 Values > 1.0 indicate slower transport due to sorption. 

937 regularization_strength : float, optional 

938 Tikhonov regularization parameter λ. See 

939 :func:`gwtransport.advection.extraction_to_infiltration` for details. 

940 Default is 1e-10. 

941 spinup : {'constant'} or None, optional 

942 Spin-up policy (default ``'constant'``). ``'constant'`` extends tedges by 

943 100 years on each side so that output bins near the boundary are always 

944 informed. ``None`` disables the extension; output bins without sufficient 

945 upstream data become NaN. Float fraction-threshold mode is not implemented 

946 and raises ``NotImplementedError``. 

947 

948 Returns 

949 ------- 

950 numpy.ndarray 

951 Bin-averaged concentration in the infiltrating water. Same units as cout. 

952 Length equals len(tedges) - 1. NaN values indicate time periods 

953 with no valid contributions from the extraction data. 

954 

955 Raises 

956 ------ 

957 ValueError 

958 If input dimensions are inconsistent, if diffusivity is negative, 

959 or if aquifer_pore_volumes and streamline_length have different lengths. 

960 

961 See Also 

962 -------- 

963 infiltration_to_extraction : Forward operation (convolution) 

964 gwtransport.advection.extraction_to_infiltration : Pure advection (no dispersion) 

965 :ref:`concept-dispersion-scales` : Macrodispersion vs microdispersion 

966 

967 Notes 

968 ----- 

969 The algorithm builds the forward coefficient matrix ``W_forward`` (same as 

970 used by :func:`infiltration_to_extraction`) and solves ``W_forward @ cin = cout`` 

971 using :func:`gwtransport.utils.solve_inverse_transport` (which builds the 

972 reverse regularization target and then applies Tikhonov regularization). 

973 This ensures mathematical consistency between forward and inverse operations. 

974 

975 NaN values in ``cout`` mark measurement gaps (e.g. sparse lab samples). 

976 Their rows are excluded from the Tikhonov solve, matching 

977 :func:`gwtransport.deposition.extraction_to_deposition`; cin bins 

978 constrained only by gapped ``cout`` bins are returned as NaN. 

979 

980 Examples 

981 -------- 

982 Basic usage with constant flow: 

983 

984 >>> import pandas as pd 

985 >>> import numpy as np 

986 >>> from gwtransport.diffusion import extraction_to_infiltration 

987 >>> 

988 >>> # Create time edges: tedges for cin/flow, cout_tedges for cout 

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

990 >>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D") 

991 >>> 

992 >>> # Extracted concentration and constant flow 

993 >>> cout = np.zeros(len(cout_tedges) - 1) 

994 >>> cout[5:10] = 1.0 # Observed pulse at extraction 

995 >>> flow = np.ones(len(tedges) - 1) * 100.0 # 100 m³/day 

996 >>> 

997 >>> # Single pore volume of 500 m³, travel distance 100 m 

998 >>> aquifer_pore_volumes = np.array([500.0]) 

999 >>> streamline_length = np.array([100.0]) 

1000 >>> 

1001 >>> # Reconstruct infiltration concentration 

1002 >>> cin = extraction_to_infiltration( 

1003 ... cout=cout, 

1004 ... flow=flow, 

1005 ... tedges=tedges, 

1006 ... cout_tedges=cout_tedges, 

1007 ... aquifer_pore_volumes=aquifer_pore_volumes, 

1008 ... streamline_length=streamline_length, 

1009 ... molecular_diffusivity=1e-4,