Coverage for src/gwtransport/diffusion_fast_fast.py: 100%

1""" 

2Fast *approximate* 1D advection-dispersion transport (Kreft-Zuber flux concentration). 

3 

4This module shares the conceptual model of :mod:`gwtransport.diffusion` and 

5:mod:`gwtransport.diffusion_fast` -- advection with microdispersion (``alpha_L``) and molecular 

6diffusion (``D_m``) along orthogonal (Cartesian) flow paths, one independent streamtube per aquifer 

7pore volume, the spread across the pore volume distribution providing macrodispersion, and linear 

8sorption via the retardation factor. It 

9targets the bin-averaged Kreft-Zuber (1978) flux concentration ``C_F`` on the streamtube bundle, but 

10trades exactness for a single fast (~1.5 ms) native-grid evaluation that does not depend on the flow 

11being constant. 

12It is **approximate**: where :mod:`gwtransport.diffusion_fast` reproduces the quadrature reference to 

13machine precision, this module is accurate to ~3e-4 in the common regime and degrades in a 

14documented corner (below). When you need machine precision, use :mod:`gwtransport.diffusion_fast`. 

15 

16How it works -- an operator split in two coordinates 

17---------------------------------------------------- 

18 

19The moving-frame dispersion product ``D_t = D_m*tau + alpha_L*xi`` mixes a *time* term (molecular 

20diffusion ``D_m*tau``) and a *volume* term (microdispersion ``alpha_L*xi``). The two are split into 

21the coordinate each is stationary in, so the dominant part is built once and is flow-independent: 

22 

231. **Advection + macrodispersion + microdispersion** are the *exact* skewed ``D_m=0`` Kreft-Zuber 

24 breakthrough, applied banded on the **native cumulative-volume grid**. The whole aquifer pore 

25 volume distribution (APVD) is pre-summed into a single 1D antiderivative ``Ibar(dV)`` -- exact 

26 for any APVD shape -- finely sampled once and read back by interpolation. This part is 

27 volume-stationary, hence **flow-independent** (constant and strongly variable flow alike). 

282. **Molecular diffusion** is a symmetric **time-domain Gaussian** applied to the outlet signal 

29 (variance ``2*D_m*tau_bt*(R*Vbar/L)^2/Q^2``). This is the only modelling approximation: the true 

30 Kreft-Zuber molecular breakthrough is skewed, and at realistic (sub-bin) spreading the Gaussian 

31 is nearly a no-op, so the molecular term is dropped rather than skewed. 

32 

33``tedges`` need **not** be regularly spaced and ``cout_tedges`` need not equal ``tedges`` (supply 

34``flow_out`` when they differ): step 1 runs on the native cumulative-volume grid for any spacing. 

35Only the molecular Gaussian assumes a roughly regular grid -- it convolves in bin-index space using 

36the mean bin width -- so a strongly irregular grid adds a small extra error to the (usually 

37sub-dominant) molecular term; use :mod:`gwtransport.diffusion_fast` for the molecular-dominated + 

38irregular-grid corner. 

39 

40Accuracy (vs :mod:`gwtransport.diffusion_fast`, flow-independent unless noted) 

41------------------------------------------------------------------------------ 

42 

43- **~3e-4 whenever microdispersion is present** (``alpha_L > 0`` -- the typical groundwater 

44 regime, Peclet number >> 1), constant *and* variable flow, for realistic solute diffusivities 

45 (``D_m`` ~ 1e-4) or ``R = 1``. Here molecular diffusion is sub-dominant, so approximating it 

46 barely matters. This survives retardation for a typical APVD (measured <~1e-3 up to ``R = 3``). 

47- In the **molecular-diffusion-dominated** corner (``alpha_L`` ~ 0): ~1e-4 for smooth inputs, but 

48 degrading to ~1e-2 for sharp inputs (and ~5e-2 for sharp inputs with a very wide / bimodal APVD or 

49 a large single pore volume), because the symmetric time-Gaussian cannot reproduce the skewed 

50 molecular breakthrough. **Retardation enlarges this corner:** the Gaussian's variance scales as 

51 ``sigma_t^2 ~ D_m * R^3``, so ``R > 1`` reaches the ~1e-2 looseness at a smaller ``D_m`` -- a 

52 sharp input at ``D_m = 0.01`` degrades from ~8e-4 at ``R = 1`` to ~1.7e-2 at ``R = 2`` and ~3.5e-2 

53 at ``R = 3``. **Use** :mod:`gwtransport.diffusion_fast` **for exact results in this regime** -- 

54 in particular for heat transport (``R > 1`` with a large ``D_m``). 

55 

56The inverse (:func:`extraction_to_infiltration`) deconvolves the *same* approximate operator the 

57forward applies. It assembles ``W = G . M`` directly in banded form (one ``Ibar`` gather plus a 

58sparse ``G . M`` product -- no per-pore-volume closed-form loop, no dense ``(n_cout, n_cin)`` matrix) 

59and solves it with banded Tikhonov regularisation (banded Cholesky, ``O(n * band**2)``), so it is 

60much faster than :mod:`gwtransport.diffusion_fast`'s reverse, especially for many streamtubes. 

61Inverting exactly the forward operator makes a round trip self-consistent: it recovers the input up 

62to the deconvolution conditioning, with the only error being the forward operator's approximation of 

63:mod:`gwtransport.diffusion_fast` (use that module when the approximation is unacceptable). 

64 

65Available functions: 

66 

67- :func:`infiltration_to_extraction` -- forward transport (approximate). 

68- :func:`extraction_to_infiltration` -- inverse via banded Tikhonov regularisation (approximate). 

69- :func:`gamma_infiltration_to_extraction` -- gamma-distributed APVD (forward). 

70- :func:`gamma_extraction_to_infiltration` -- same, inverse. 

71 

72References 

73---------- 

74Kreft, A., & Zuber, A. (1978). On the physical meaning of the dispersion equation and its 

75solutions for different initial and boundary conditions. Chemical Engineering Science, 

7633(11), 1471-1480. 

77 

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

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

80""" 

81 

82import numpy as np 

83import numpy.typing as npt 

84import pandas as pd 

85from scipy.ndimage import gaussian_filter1d 

86from scipy.sparse import coo_array 

87 

88from gwtransport import gamma 

89from gwtransport._diffusion_shared import ( 

90 _DT_FLOOR, 

91 _EPSILON_COEFF_SUM, 

92 _breakthrough_antideriv, 

93 _broadcast_to_pore_volumes, 

94 _cout_cumulative_volume, 

95 _extend_tedges_flag, 

96 _solve_reverse_banded, 

97 _validate_inputs, 

98) 

99from gwtransport._time import dt_to_days, tedges_to_days 

100from gwtransport.diffusion_fast import _DEFAULT_SATURATION_THRESHOLD 

101from gwtransport.residence_time import fraction_explained_full 

102from gwtransport.utils import cumulative_flow_volume 

103 

104# Samples per native bin used to discretise the 1D breakthrough antiderivative ``Ibar``. Higher = 

105# more accurate ``Ibar`` interpolation (advection+micro error ~ O(1/_KERNEL_FINE^2)) at the cost of 

106# a larger one-time precompute; 16 gives ~1e-4 in the advection+micro part, which is below the 

107# molecular-Gaussian floor, so it is not exposed as a user knob. 

108_KERNEL_FINE = 16 

109 

110# Upper bound on the number of fine ``Ibar`` samples. Caps the one-time precompute when the 

111# breakthrough band spans far more than the record (e.g. tiny flow -> enormous front offset); the 

112# affected output bins are masked by residence time anyway, so coarsening the band there is benign. 

113_MAX_KERNEL_SAMPLES = 20000 

114 

115 

116def _summed_antideriv( 

117 *, 

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

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

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

121 retardation_factor: float, 

122 mean_bin_volume: float, 

123 saturation_threshold: float, 

124) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating], float, float, float]: 

125 r"""Precompute the APVD-summed breakthrough antiderivative ``Ibar(dV)`` on a fine 1D grid. 

126 

127 For one streamtube the bin-averaged ``D_m=0`` flux fraction over a cout bin equals the second 

128 difference of the antiderivative ``I(x)`` (:func:`gwtransport._diffusion_shared._breakthrough_antideriv`) 

129 of the resident concentration in the breakthrough coordinate ``x = (dV - r_vpv)*L/r_vpv`` 

130 (``r_vpv = R*V_pore``). The antiderivative *with respect to cumulative volume* ``dV`` is 

131 ``(r_vpv/L)*I(x(dV))``; averaging it over the APVD gives a single 1D function 

132 

133 .. math:: 

134 

135 \bar I(\Delta V) = \operatorname{mean}_{pv}\Bigl[\tfrac{r_{vpv}}{L}\, 

136 I\bigl((\Delta V - r_{vpv})L/r_{vpv}\bigr)\Bigr],\quad D_t = \alpha_L\,\xi, 

137 

138 whose edge-differences reproduce the per-streamtube-averaged ``C_F`` **exactly** for any APVD 

139 (the ``r_vpv/L`` Jacobian and the cout-bin volume normalisation cancel the per-streamtube ``dx``). 

140 Only the interpolation of ``Ibar`` is approximate. 

141 

142 The grid is **uniform** over the breakthrough band ``[off_lo, off_hi]`` (front ``r_vpv`` plus the 

143 conservative ``alpha_L`` dispersion smear, unioned over streamtubes) plus a margin, sampled at 

144 ``mean_bin_volume / _KERNEL_FINE`` (uniformity lets :func:`_eval_antideriv` interpolate by 

145 fractional indexing instead of a per-point search). For ``alpha_L > 0`` ``Ibar`` is smooth and 

146 the interpolation error is ``O(1/_KERNEL_FINE^2)``; at ``alpha_L = 0`` it has a kink at each 

147 ``r_vpv`` and the error is ``O(1/_KERNEL_FINE)`` (sub-1e-3, well below the molecular floor). 

148 

149 Returns 

150 ------- 

151 grid : ndarray 

152 Cumulative-volume offsets ``dV`` (uniformly spaced, strictly increasing). 

153 ibar : ndarray 

154 ``Ibar`` sampled at ``grid``. 

155 mean_r_vpv : float 

156 ``R * mean(V_pore)`` -- the saturated offset (``Ibar -> dV - mean_r_vpv`` above the band). 

157 off_lo, off_hi : float 

158 Lower / upper cumulative-volume offset of the breakthrough band. 

159 """ 

160 r_vpv = retardation_factor * aquifer_pore_volumes 

161 mean_r_vpv = float(r_vpv.mean()) 

162 u = saturation_threshold 

163 # D_m=0 dispersion half-widths in the breakthrough coordinate x (front D_t = alpha_L*L): the 

164 # pre side shrinks (|x| = U*2*sqrt(alpha_L*L)); the post side grows with slope alpha_L, giving 

165 # the quadratic root x = 2U^2*alpha_L + 2U*sqrt(U^2*alpha_L^2 + alpha_L*L). Mapped to volume 

166 # offsets via r_vpv/L and unioned over streamtubes (conservative, never under-covers the band). 

167 pre_x = u * 2.0 * np.sqrt(longitudinal_dispersivity * streamline_length) 

168 post_x = 2.0 * u * u * longitudinal_dispersivity + 2.0 * u * np.sqrt( 

169 u * u * longitudinal_dispersivity**2 + longitudinal_dispersivity * streamline_length 

170 ) 

171 off_lo = float(np.min(r_vpv - (r_vpv / streamline_length) * pre_x)) 

172 off_hi = float(np.max(r_vpv + (r_vpv / streamline_length) * post_x)) 

173 

174 margin = 4.0 * mean_bin_volume 

175 span = (off_hi - off_lo) + 2.0 * margin 

176 dv_fine = mean_bin_volume / _KERNEL_FINE 

177 if span / dv_fine > _MAX_KERNEL_SAMPLES: 

178 dv_fine = span / _MAX_KERNEL_SAMPLES 

179 grid = np.arange(off_lo - margin, off_hi + margin + dv_fine, dv_fine) 

180 

181 x = (grid[None, :] - r_vpv[:, None]) * streamline_length[:, None] / r_vpv[:, None] 

182 dt_var = np.maximum(longitudinal_dispersivity[:, None] * np.maximum(x + streamline_length[:, None], 0.0), _DT_FLOOR) 

183 antideriv = _breakthrough_antideriv(x, dt_var) 

184 ibar = ((r_vpv[:, None] / streamline_length[:, None]) * antideriv).mean(axis=0) 

185 return grid, ibar, mean_r_vpv, off_lo, off_hi 

186 

187 

188def _eval_antideriv( 

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

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

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

192 mean_r_vpv: float, 

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

194 """Evaluate ``Ibar`` at arbitrary cumulative-volume offsets. 

195 

196 Linear interpolation on the *uniform* precomputed grid by fractional indexing (no per-point 

197 search -- the gather is evaluated at ``O(N*band)`` points, so this dominates the runtime). 

198 Below the grid ``Ibar = 0`` (not broken through); above it ``Ibar = dV - mean_r_vpv`` (saturated, 

199 the exact linear asymptote). 

200 

201 Returns 

202 ------- 

203 ndarray 

204 ``Ibar(dv)`` with the same shape as ``dv``. 

205 """ 

206 g0 = grid[0] 

207 dstep = grid[1] - grid[0] 

208 f = (dv - g0) / dstep 

209 # astype truncates toward zero; equals floor on f >= 0 (the only regime that survives the 

210 # dv < g0 override below and the lower clip), so the floor call is redundant. Precompute the 

211 # segment slopes once so the O(N*band) gather reads a single array instead of differencing two. 

212 slopes = np.diff(ibar) 

213 i0 = np.clip(f.astype(np.intp), 0, grid.size - 2) 

214 out = ibar[i0] + (f - i0) * slopes[i0] 

215 out = np.where(dv < g0, 0.0, out) 

216 return np.where(dv > grid[-1], dv - mean_r_vpv, out) 

217 

218 

219def _advection_micro_band( 

220 *, 

221 cumulative_volume_at_cin: npt.NDArray[np.floating], 

222 cumulative_volume_at_cout: npt.NDArray[np.floating], 

223 grid: npt.NDArray[np.floating], 

224 ibar: npt.NDArray[np.floating], 

225 mean_r_vpv: float, 

226 off_lo: float, 

227 off_hi: float, 

228 extend: bool, 

229) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.intp]]: 

230 """Banded ``D_m=0`` advection+macro+micro coefficients on the native cumulative-volume grid. 

231 

232 For each cout bin, gathers the band of cin edges whose breakthrough offset falls in 

233 ``[off_lo, off_hi]`` (a conservative fixed window from ``searchsorted``, mirroring 

234 :func:`gwtransport.diffusion_fast._closed_form_coeff_matrix`), evaluates ``Ibar`` at the native 

235 edge offsets, and telescopes the edge-differences into per-cin-bin coefficients. With 

236 ``extend`` the cin axis is extended by one wide virtual bin each side carrying the constant 

237 boundary value (``cin[0]`` / ``cin[-1]``), reproducing the 100-year warm-start. 

238 

239 The coefficients depend only on the volume grid (not on ``cin``), so this band is built once and 

240 shared: :func:`infiltration_to_extraction` applies it to the (extended) ``cin``, and 

241 :func:`extraction_to_infiltration` folds it onto the real cin axis to assemble the banded 

242 operator it deconvolves -- guaranteeing forward and reverse use the same operator. 

243 

244 Returns 

245 ------- 

246 coeff : ndarray, shape (n_cout, band) 

247 Per-(cout bin, band slot) coefficient. 

248 cin_bin : ndarray of int, shape (n_cout, band) 

249 Column index of each coefficient on the (warm-start-extended when ``extend``) cin axis: 

250 ``cin_ext[cin_bin]`` for the forward, ``clip(cin_bin - int(extend), 0, n_cin - 1)`` for the 

251 real cin axis in the reverse. 

252 """ 

253 vi = cumulative_volume_at_cin 

254 vc = cumulative_volume_at_cout 

255 if extend: 

256 big = abs(off_hi) + abs(off_lo) + (vc[-1] - vc[0]) + (vi[-1] - vi[0]) + 1.0 

257 vi_ext = np.concatenate([[vi[0] - big], vi, [vi[-1] + big]]) 

258 n_cin_ext = vi.size + 1 # (vi.size - 1) real bins + 2 virtual boundary bins 

259 else: 

260 vi_ext = vi 

261 n_cin_ext = vi.size - 1 

262 

263 vc_lo, vc_hi = vc[:-1], vc[1:] 

264 w = vc_hi - vc_lo 

265 vc_mid = 0.5 * (vc_lo + vc_hi) 

266 # Conservative two-sided fixed window: center on the front, widen to the farthest cin edge whose 

267 # offset is still inside the band on either side. The +1 absorbs the edge consumed by the 

268 # telescoping difference and searchsorted rounding (an under-sized window silently drops mass). 

269 center = np.searchsorted(vi_ext, vc_mid - mean_r_vpv) 

270 j_lo = np.searchsorted(vi_ext, vc_lo - off_hi, side="left") 

271 j_hi = np.searchsorted(vi_ext, vc_hi - off_lo, side="right") 

272 hw_lo = int(np.max(center - j_lo)) + 1 

273 hw_hi = int(np.max(j_hi - center)) + 1 

274 

275 cols = center[:, None] + np.arange(-hw_lo, hw_hi + 1)[None, :] 

276 vi_band = vi_ext[np.clip(cols, 0, len(vi_ext) - 1)] 

277 ibar_hi = _eval_antideriv(vc_hi[:, None] - vi_band, grid, ibar, mean_r_vpv) 

278 ibar_lo = _eval_antideriv(vc_lo[:, None] - vi_band, grid, ibar, mean_r_vpv) 

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

280 frac_edge = np.where(w[:, None] > 0.0, (ibar_hi - ibar_lo) / w[:, None], 0.0) 

281 coeff = frac_edge[:, :-1] - frac_edge[:, 1:] 

282 cin_bin = np.clip(cols[:, :-1], 0, n_cin_ext - 1) 

283 return coeff, cin_bin 

284 

285 

286def _build_forward_operator( 

287 *, 

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

289 tedges: pd.DatetimeIndex, 

290 cout_tedges: pd.DatetimeIndex, 

291 flow_out: npt.NDArray[np.floating] | None, 

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

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

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

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

296 retardation_factor: float, 

297 extend: bool, 

298 saturation_threshold: float, 

299) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.intp], float, npt.NDArray[np.bool_]] | None: 

300 r"""Build the cin-independent pieces of the approximate banded forward operator ``W = G . M``. 

301 

302 Both directions share this build: the advection+macro+micro band ``M`` (``coeff`` / ``cin_bin``, 

303 :func:`_advection_micro_band`), the molecular time-Gaussian width ``sigma_bins`` (the operator 

304 ``G``), and the residence-time ``valid`` mask. Because the band depends only on the volume grid 

305 (not on ``cin`` / ``cout``), forward transport and reverse deconvolution operate on exactly the 

306 same operator. 

307 

308 Returns 

309 ------- 

310 coeff : ndarray, shape (n_cout, band) 

311 Banded ``M`` coefficients. 

312 cin_bin : ndarray of int, shape (n_cout, band) 

313 Column indices of ``coeff`` on the warm-start-extended cin axis. 

314 sigma_bins : float 

315 Molecular Gaussian width in output-bin units (0 -> ``G`` is the identity). 

316 valid : ndarray of bool, shape (n_cout,) 

317 Output bins with residence time finite at both edges (complete breakthrough). 

318 

319 None 

320 Returned instead of the tuple when there is no through-flow (nothing breaks through). 

321 """ 

322 tedges_days = tedges_to_days(tedges) 

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

324 cumulative_volume_at_cin = cumulative_flow_volume(flow, dt_to_days(tedges)) 

325 total_volume = float(cumulative_volume_at_cin[-1]) 

326 if total_volume <= 0.0: 

327 return None 

328 

329 cumulative_volume_at_cout = _cout_cumulative_volume( 

330 flow_out=flow_out, 

331 cout_tedges=cout_tedges, 

332 cout_tedges_days=cout_tedges_days, 

333 tedges_days=tedges_days, 

334 cumulative_volume_at_cin=cumulative_volume_at_cin, 

335 ) 

336 

337 n_cout = len(cout_tedges) - 1 

338 mean_bin_volume = total_volume / len(flow) 

339 grid, ibar, mean_r_vpv, off_lo, off_hi = _summed_antideriv( 

340 aquifer_pore_volumes=aquifer_pore_volumes, 

341 streamline_length=streamline_length, 

342 longitudinal_dispersivity=longitudinal_dispersivity, 

343 retardation_factor=retardation_factor, 

344 mean_bin_volume=mean_bin_volume, 

345 saturation_threshold=saturation_threshold, 

346 ) 

347 coeff, cin_bin = _advection_micro_band( 

348 cumulative_volume_at_cin=cumulative_volume_at_cin, 

349 cumulative_volume_at_cout=cumulative_volume_at_cout, 

350 grid=grid, 

351 ibar=ibar, 

352 mean_r_vpv=mean_r_vpv, 

353 off_lo=off_lo, 

354 off_hi=off_hi, 

355 extend=extend, 

356 ) 

357 

358 # Molecular diffusion: a single mean-streamtube time-domain Gaussian on the outlet signal. 

359 # sigma_t^2 = 2*D_m*tau_bt*(r_vpv/L)^2/Q^2, with tau_bt = R*mean(V_pore)/Q the mean breakthrough 

360 # time and Q = total_volume/total_time the flow-weighted mean throughflow. sigma_t (days) is 

361 # converted with the OUTPUT-grid mean bin width (not the flow grid -- they differ for a coarse 

362 # cout grid); a smear wider than the record cannot be resolved. Cap sigma_bins so the forward's 

363 # truncation radius int(6*sigma_bins + 0.5) cannot exceed n_cout - 1 -- this matches the explicit 

364 # min(..., n_cout - 1) clip the reverse applies in _banded_forward_matrix, so forward and reverse 

365 # always use the identical molecular kernel (the cap is unreachable for realistic D_m). 

366 total_days = float(tedges_days[-1] - tedges_days[0]) 

367 q_mean = total_volume / total_days 

368 tau_bt = retardation_factor * float(aquifer_pore_volumes.mean()) / q_mean 

369 sigma_t2 = ( 

370 2.0 * float(molecular_diffusivity.mean()) * tau_bt * (mean_r_vpv / float(streamline_length.mean())) ** 2 

371 ) / q_mean**2 

372 mean_cout_dt = (cout_tedges_days[-1] - cout_tedges_days[0]) / n_cout 

373 sigma_bins = min(float(np.sqrt(sigma_t2)) / mean_cout_dt, (n_cout - 1) / 6.0) 

374 

375 # Mask bins beyond the data range (and, without warm-start, incompletely-broken-through spin-up 

376 # bins). residence_time uses the extended grid when warm-starting so spin-up bins stay valid. 

377 work_tedges = tedges 

378 if extend: 

379 edge_values = tedges.to_numpy().copy() 

380 delta = np.timedelta64(36500, "D") 

381 edge_values[0] -= delta 

382 edge_values[-1] += delta 

383 work_tedges = pd.DatetimeIndex(edge_values) 

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

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

386 valid = np.all( 

387 fraction_explained_full( 

388 flow=flow, 

389 tedges=work_tedges, 

390 cout_tedges=cout_tedges, 

391 aquifer_pore_volumes=aquifer_pore_volumes, 

392 retardation_factor=retardation_factor, 

393 direction="extraction_to_infiltration", 

394 ) 

395 >= 1.0, 

396 axis=0, 

397 ) 

398 return coeff, cin_bin, sigma_bins, valid 

399 

400 

401def _banded_forward_matrix( 

402 *, 

403 coeff: npt.NDArray[np.floating], 

404 cin_bin: npt.NDArray[np.intp], 

405 extend: bool, 

406 n_cin: int, 

407 sigma_bins: float, 

408) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.intp]]: 

409 """Assemble ``W = G . M`` as a per-row contiguous band for ``_solve_reverse_banded``. 

410 

411 ``M`` (advection + macro + microdispersion) is scattered from the native band onto the real cin 

412 axis, folding the warm-start virtual columns into the boundary bins 

413 (``clip(cin_bin - int(extend), 0, n_cin - 1)``, so ``M @ cin`` equals the forward's 

414 ``coeff @ cin_ext`` exactly). ``G`` is the molecular time-Gaussian along the output-bin axis (the 

415 same ``mode="nearest"`` kernel the forward applies with :func:`scipy.ndimage.gaussian_filter1d`). 

416 The returned band carries the forward operator verbatim; ``_solve_reverse_banded`` masks the 

417 spin-up rows/columns and normalizes, so a forward-then-inverse round trip is self-consistent. 

418 

419 Returns 

420 ------- 

421 band_vals : ndarray, shape (n_cout, full_band) 

422 Forward weights in banded layout (explicit zeros outside each row's support). 

423 col_start : ndarray of int, shape (n_cout,) 

424 First real-cin column of each row's band. 

425 """ 

426 n_cout = coeff.shape[0] 

427 rows = np.broadcast_to(np.arange(n_cout)[:, None], coeff.shape) 

428 real_col = np.clip(cin_bin - int(extend), 0, n_cin - 1) 

429 m_mat = coo_array((coeff.ravel(), (rows.ravel(), real_col.ravel())), shape=(n_cout, n_cin)).tocsr() 

430 

431 lw = min(int(6.0 * sigma_bins + 0.5), n_cout - 1) if sigma_bins > 0.0 else 0 

432 if lw > 0: 

433 offsets = np.arange(-lw, lw + 1) 

434 kernel = np.exp(-0.5 * (offsets / sigma_bins) ** 2) 

435 kernel /= kernel.sum() 

436 g_rows = np.repeat(np.arange(n_cout), offsets.size) 

437 g_cols = np.clip(np.arange(n_cout)[:, None] + offsets[None, :], 0, n_cout - 1).ravel() 

438 g_mat = coo_array((np.tile(kernel, n_cout), (g_rows, g_cols)), shape=(n_cout, n_cout)).tocsr() 

439 w_mat = (g_mat @ m_mat).tocsr() 

440 else: 

441 w_mat = m_mat 

442 

443 # CSR -> contiguous per-row band. Each W row is the union of overlapping contiguous M bands, so 

444 # it stays contiguous (a flow spike that briefly opens a gap only adds explicit interior zeros). 

445 w_mat.sort_indices() 

446 indptr, indices, data = w_mat.indptr, w_mat.indices, w_mat.data 

447 row_counts = np.diff(indptr) 

448 nonempty = row_counts > 0 

449 col_start = np.zeros(n_cout, dtype=np.intp) 

450 last_col = np.zeros(n_cout, dtype=np.intp) 

451 col_start[nonempty] = indices[indptr[:-1][nonempty]] 

452 last_col[nonempty] = indices[indptr[1:][nonempty] - 1] 

453 full_band = int((last_col[nonempty] - col_start[nonempty] + 1).max()) if nonempty.any() else 1 

454 

455 band_vals = np.zeros((n_cout, full_band)) 

456 rows_of_nz = np.repeat(np.arange(n_cout), row_counts) 

457 band_vals[rows_of_nz, indices - col_start[rows_of_nz]] = data 

458 return band_vals, col_start 

459 

460 

461def infiltration_to_extraction( 

462 *, 

463 cin: npt.ArrayLike, 

464 flow: npt.ArrayLike, 

465 tedges: pd.DatetimeIndex, 

466 cout_tedges: pd.DatetimeIndex, 

467 aquifer_pore_volumes: npt.ArrayLike, 

468 streamline_length: npt.NDArray[np.floating] | float, 

469 molecular_diffusivity: npt.NDArray[np.floating] | float, 

470 longitudinal_dispersivity: npt.NDArray[np.floating] | float, 

471 retardation_factor: float = 1.0, 

472 flow_out: npt.ArrayLike | None = None, 

473 spinup: str | None = "constant", 

474 saturation_threshold: float = _DEFAULT_SATURATION_THRESHOLD, 

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

476 """Compute extracted concentration with advection, microdispersion, and molecular diffusion (approximate). 

477 

478 Fast *approximate* counterpart of :func:`gwtransport.diffusion_fast.infiltration_to_extraction`. 

479 The advection + macrodispersion + microdispersion (``alpha_L``) part is the exact skewed 

480 ``D_m=0`` Kreft-Zuber breakthrough applied on the native cumulative-volume grid; molecular 

481 diffusion (``D_m``) is a symmetric time-domain Gaussian. The result is flow-independent and 

482 accurate to ~3e-4 whenever ``alpha_L > 0`` (the typical regime) for realistic solute ``D_m`` 

483 (~1e-4) or ``R = 1``. It loosens to ~1e-2 (sharp inputs) in the molecular-diffusion-dominated 

484 corner (``alpha_L`` ~ 0), and retardation enlarges that corner because the Gaussian variance 

485 grows as ``D_m * R^3`` (a sharp input at ``D_m = 0.01`` reaches ~1.7e-2 at ``R = 2`` and ~3.5e-2 

486 at ``R = 3``). For machine precision -- or heat transport with ``R > 1`` and large ``D_m`` -- use 

487 :mod:`gwtransport.diffusion_fast`. 

488 

489 Parameters 

490 ---------- 

491 cin : array-like 

492 Concentration of the compound in the infiltrating water. Length ``len(tedges) - 1``. 

493 flow : array-like 

494 Flow rate of water in the aquifer [m³/day]. Length ``len(tedges) - 1``. 

495 tedges : pandas.DatetimeIndex 

496 Time edges for cin and flow data. Length ``len(cin) + 1``. 

497 cout_tedges : pandas.DatetimeIndex 

498 Time edges for output data bins. Length ``len(output) + 1``. 

499 aquifer_pore_volumes : array-like 

500 Aquifer pore volumes [m³] -- one independent streamtube per entry. Any distribution shape 

501 (the APVD is pre-summed exactly). 

502 streamline_length : float or ndarray 

503 Travel distance L [m]: a scalar (shared by all streamtubes) or an array with one 

504 value per aquifer pore volume. Must be positive. 

505 molecular_diffusivity : float or ndarray 

506 Effective molecular diffusivity D_m [m²/day]: scalar or one value per pore volume. 

507 Must be non-negative. 

508 longitudinal_dispersivity : float or ndarray 

509 Longitudinal dispersivity alpha_L [m] (microdispersion): scalar or one value per pore volume. 

510 Must be non-negative. 

511 retardation_factor : float, optional 

512 Retardation factor (default 1.0). Values > 1.0 indicate slower transport. 

513 flow_out : array-like or None, optional 

514 Extraction flow rate [m³/day] on the output grid (aligned to ``cout_tedges``, 

515 length ``len(cout_tedges) - 1``). Required when ``cout_tedges`` differs from ``tedges``; 

516 may be omitted only when ``cout_tedges`` equals ``tedges``. Default None. 

517 spinup : {"constant"} | None, optional 

518 ``"constant"`` (default) extends ``tedges`` by 100 years on each side so a constant 

519 warm-start fills the left-edge spin-up region; ``None`` leaves spin-up cout as NaN. 

520 saturation_threshold : float, optional 

521 Breakthrough-band cutoff ``U`` (default 7.0). Sets how far into the breakthrough tail the 

522 banded build reaches; see :func:`gwtransport.diffusion_fast.infiltration_to_extraction`. 

523 

524 Returns 

525 ------- 

526 numpy.ndarray 

527 Bin-averaged Kreft-Zuber flux concentration ``C_F`` in the extracted water. Length 

528 ``len(cout_tedges) - 1``. NaN where no infiltration data has broken through. 

529 

530 See Also 

531 -------- 

532 gwtransport.diffusion_fast.infiltration_to_extraction : Exact (machine-precision) counterpart; 

533 use it when approximation is unacceptable, especially in the molecular-dominant regime. 

534 gwtransport.diffusion.infiltration_to_extraction : Quadrature reference. 

535 extraction_to_infiltration : Inverse operation (deconvolves this same operator). 

536 :ref:`concept-dispersion-scales` : Macrodispersion vs microdispersion. 

537 """ 

538 tedges = pd.DatetimeIndex(tedges) 

539 cout_tedges = pd.DatetimeIndex(cout_tedges) 

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

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

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

543 if flow_out is not None: 

544 flow_out = np.asarray(flow_out, dtype=float) 

545 

546 _validate_inputs( 

547 cin_or_cout=cin, 

548 flow=flow, 

549 tedges=tedges, 

550 cout_tedges=cout_tedges, 

551 aquifer_pore_volumes=aquifer_pore_volumes, 

552 streamline_length=streamline_length, 

553 molecular_diffusivity=molecular_diffusivity, 

554 longitudinal_dispersivity=longitudinal_dispersivity, 

555 retardation_factor=retardation_factor, 

556 is_forward=True, 

557 flow_out=flow_out, 

558 ) 

559 

560 n_pv = len(aquifer_pore_volumes) 

561 streamline_length = _broadcast_to_pore_volumes(streamline_length, n_pv) 

562 molecular_diffusivity = _broadcast_to_pore_volumes(molecular_diffusivity, n_pv) 

563 longitudinal_dispersivity = _broadcast_to_pore_volumes(longitudinal_dispersivity, n_pv) 

564 

565 n_cout = len(cout_tedges) - 1 

566 extend = _extend_tedges_flag(spinup) 

567 operator = _build_forward_operator( 

568 flow=flow, 

569 tedges=tedges, 

570 cout_tedges=cout_tedges, 

571 flow_out=flow_out, 

572 aquifer_pore_volumes=aquifer_pore_volumes, 

573 streamline_length=streamline_length, 

574 molecular_diffusivity=molecular_diffusivity, 

575 longitudinal_dispersivity=longitudinal_dispersivity, 

576 retardation_factor=retardation_factor, 

577 extend=extend, 

578 saturation_threshold=saturation_threshold, 

579 ) 

580 if operator is None: 

581 # No through-flow: nothing breaks through (matches diffusion_fast's all-NaN result). 

582 return np.full(n_cout, np.nan) 

583 coeff, cin_bin, sigma_bins, valid = operator 

584 

585 # Apply the advection+macro+micro band M to the (warm-start-extended) cin, then the molecular 

586 # time-Gaussian G. Output bins with no through-flow carry a hard 0 in cout_micro; a plain 

587 # Gaussian would smear those zeros into the valid neighbours. Convolve mask-aware instead -- 

588 # G(cout_micro*support)/G(support) -- so gap bins act as missing (not zero) data. This equals a 

589 # plain G where support is all-True (constants stay constant, incl. exactly across a gap). 

590 cin_ext = np.concatenate([[cin[0]], cin, [cin[-1]]]) if extend else cin 

591 cout_micro = np.einsum("kb,kb->k", coeff, cin_ext[cin_bin]) 

592 support = coeff.sum(axis=1) >= _EPSILON_COEFF_SUM 

593 if sigma_bins == 0.0: 

594 cout = cout_micro 

595 else: 

596 num = gaussian_filter1d(cout_micro * support, sigma_bins, mode="nearest", truncate=6.0) 

597 den = gaussian_filter1d(support.astype(float), sigma_bins, mode="nearest", truncate=6.0) 

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

599 cout = num / den # den > 0 wherever support is True; 0/0 gap bins are masked out below 

600 return np.where(support & valid, cout, np.nan) 

601 

602 

603def extraction_to_infiltration( 

604 *, 

605 cout: npt.ArrayLike, 

606 flow: npt.ArrayLike, 

607 tedges: pd.DatetimeIndex, 

608 cout_tedges: pd.DatetimeIndex, 

609 aquifer_pore_volumes: npt.ArrayLike, 

610 streamline_length: npt.NDArray[np.floating] | float, 

611 molecular_diffusivity: npt.NDArray[np.floating] | float, 

612 longitudinal_dispersivity: npt.NDArray[np.floating] | float, 

613 retardation_factor: float = 1.0, 

614 regularization_strength: float = 1e-10, 

615 flow_out: npt.ArrayLike | None = None, 

616 spinup: str | None = "constant", 

617 saturation_threshold: float = _DEFAULT_SATURATION_THRESHOLD, 

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

619 """Reconstruct infiltration concentration from extracted water (fast approximate deconvolution). 

620 

621 Inverts the **same** approximate operator the forward applies: it assembles ``W = G . M`` (the 

622 advection+macro+micro band ``M`` times the molecular time-Gaussian ``G``) directly in banded form 

623 and deconvolves it with banded Tikhonov regularization (``_solve_reverse_banded`` -- banded 

624 Cholesky on the normal equations, ``O(n * band**2)``). It builds ``W`` from one ``Ibar`` gather 

625 plus a sparse ``G . M`` product -- no per-pore-volume closed-form loop and no dense 

626 ``(n_cout, n_cin)`` matrix -- so it is much faster than 

627 :func:`gwtransport.diffusion_fast.extraction_to_infiltration` (which evaluates the exact 

628 breakthrough per streamtube), especially for many streamtubes. Because the deconvolved operator 

629 is exactly the forward operator, a forward-then-inverse round trip recovers ``cin`` up to the 

630 deconvolution conditioning and regularization; the approximation lives entirely in the forward 

631 operator vs :mod:`gwtransport.diffusion_fast`. 

632 

633 Parameters 

634 ---------- 

635 cout : array-like 

636 Concentration of the compound in extracted water. Length ``len(cout_tedges) - 1``. 

637 flow : array-like 

638 Flow rate of water in the aquifer [m³/day]. Length ``len(tedges) - 1``. 

639 tedges : pandas.DatetimeIndex 

640 Time edges for cin (output) and flow data. Length ``len(flow) + 1``. 

641 cout_tedges : pandas.DatetimeIndex 

642 Time edges for cout data bins. Length ``len(cout) + 1``. 

643 aquifer_pore_volumes : array-like 

644 Aquifer pore volumes [m³] -- one independent streamtube per entry. 

645 streamline_length : float or ndarray 

646 Travel distance L [m]: scalar or one value per pore volume. Must be positive. 

647 molecular_diffusivity : float or ndarray 

648 Effective molecular diffusivity D_m [m²/day]: scalar or one value per pore volume. 

649 Must be non-negative. 

650 longitudinal_dispersivity : float or ndarray 

651 Longitudinal dispersivity alpha_L [m] (microdispersion): scalar or one value per pore volume. 

652 Must be non-negative. 

653 retardation_factor : float, optional 

654 Retardation factor (default 1.0). 

655 regularization_strength : float, optional 

656 Tikhonov regularization parameter (default 1e-10). Must be strictly positive: the banded 

657 solver relies on it to make the normal equations positive-definite (it cannot return the 

658 dense ``lambda = 0`` minimum-norm solution). 

659 flow_out : array-like or None, optional 

660 Extraction flow rate [m³/day] on the output grid (aligned to ``cout_tedges``). See 

661 :func:`infiltration_to_extraction`. Default None. 

662 spinup : {"constant"} | None, optional 

663 See :func:`infiltration_to_extraction`. Default ``"constant"``. 

664 saturation_threshold : float, optional 

665 See :func:`infiltration_to_extraction`. Default 7.0. 

666 

667 Returns 

668 ------- 

669 numpy.ndarray 

670 Bin-averaged concentration in the infiltrating water. Length ``len(tedges) - 1``. 

671 NaN where no extraction data constrains the bin. 

672 

673 See Also 

674 -------- 

675 infiltration_to_extraction : Forward operation (the operator inverted here). 

676 gwtransport.diffusion_fast.extraction_to_infiltration : Exact (machine-precision) counterpart; 

677 use it when the approximation is unacceptable. 

678 :ref:`concept-dispersion-scales` : Macrodispersion vs microdispersion. 

679 """ 

680 tedges = pd.DatetimeIndex(tedges) 

681 cout_tedges = pd.DatetimeIndex(cout_tedges) 

682 cout = np.asarray(cout, dtype=float) 

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

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

685 if flow_out is not None: 

686 flow_out = np.asarray(flow_out, dtype=float) 

687 

688 _validate_inputs( 

689 cin_or_cout=cout, 

690 flow=flow, 

691 tedges=tedges, 

692 cout_tedges=cout_tedges, 

693 aquifer_pore_volumes=aquifer_pore_volumes, 

694 streamline_length=streamline_length, 

695 molecular_diffusivity=molecular_diffusivity, 

696 longitudinal_dispersivity=longitudinal_dispersivity, 

697 retardation_factor=retardation_factor, 

698 is_forward=False, 

699 flow_out=flow_out, 

700 ) 

701 

702 n_pv = len(aquifer_pore_volumes) 

703 streamline_length = _broadcast_to_pore_volumes(streamline_length, n_pv) 

704 molecular_diffusivity = _broadcast_to_pore_volumes(molecular_diffusivity, n_pv) 

705 longitudinal_dispersivity = _broadcast_to_pore_volumes(longitudinal_dispersivity, n_pv) 

706 

707 n_cin = len(tedges) - 1 

708 extend = _extend_tedges_flag(spinup) 

709 operator = _build_forward_operator( 

710 flow=flow, 

711 tedges=tedges, 

712 cout_tedges=cout_tedges, 

713 flow_out=flow_out, 

714 aquifer_pore_volumes=aquifer_pore_volumes, 

715 streamline_length=streamline_length, 

716 molecular_diffusivity=molecular_diffusivity, 

717 longitudinal_dispersivity=longitudinal_dispersivity, 

718 retardation_factor=retardation_factor, 

719 extend=extend, 

720 saturation_threshold=saturation_threshold, 

721 ) 

722 if operator is None: 

723 # No through-flow: nothing constrains the infiltration signal. 

724 return np.full(n_cin, np.nan) 

725 coeff, cin_bin, sigma_bins, valid = operator 

726 

727 band_vals, col_start = _banded_forward_matrix( 

728 coeff=coeff, cin_bin=cin_bin, extend=extend, n_cin=n_cin, sigma_bins=sigma_bins 

729 ) 

730 return _solve_reverse_banded( 

731 band_vals=band_vals, 

732 col_start=col_start, 

733 valid_cout_bins=valid, 

734 cout=cout, 

735 n_cin=n_cin, 

736 regularization_strength=regularization_strength, 

737 ) 

738 

739 

740def gamma_infiltration_to_extraction( 

741 *, 

742 cin: npt.ArrayLike, 

743 flow: npt.ArrayLike, 

744 tedges: pd.DatetimeIndex, 

745 cout_tedges: pd.DatetimeIndex, 

746 mean: float | None = None, 

747 std: float | None = None, 

748 loc: float = 0.0, 

749 alpha: float | None = None, 

750 beta: float | None = None, 

751 n_bins: int = 100, 

752 streamline_length: float, 

753 molecular_diffusivity: float, 

754 longitudinal_dispersivity: float, 

755 retardation_factor: float = 1.0, 

756 flow_out: npt.ArrayLike | None = None, 

757 spinup: str | None = "constant", 

758 saturation_threshold: float = _DEFAULT_SATURATION_THRESHOLD, 

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

760 """Compute extracted concentration for a gamma-distributed pore volume distribution (approximate). 

761 

762 Convenience wrapper around :func:`infiltration_to_extraction` that discretizes a (shifted) 

763 gamma aquifer pore-volume distribution into ``n_bins`` equal-probability streamtubes. Provide 

764 either (mean, std) or (alpha, beta); ``loc`` defaults to 0. Approximate -- see 

765 :func:`infiltration_to_extraction`. 

766 

767 Parameters 

768 ---------- 

769 cin : array-like 

770 Concentration of the compound in infiltrating water. 

771 flow : array-like 

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

773 tedges : pandas.DatetimeIndex 

774 Time edges for cin and flow data. Length ``len(cin) + 1``. 

775 cout_tedges : pandas.DatetimeIndex 

776 Time edges for output data bins. Length ``len(result) + 1``. 

777 mean, std : float, optional 

778 Mean and standard deviation of the gamma pore-volume distribution. 

779 loc : float, optional 

780 Location (minimum pore volume), ``0 <= loc < mean``. Default 0.0. 

781 alpha, beta : float, optional 

782 Shape and scale parameters of the gamma distribution (alternative to mean/std). 

783 n_bins : int, optional 

784 Number of equal-probability streamtubes. Default 100. 

785 streamline_length : float 

786 Travel distance L [m], applied to all gamma streamtubes. Must be positive. 

787 molecular_diffusivity : float 

788 Effective molecular diffusivity D_m [m²/day], applied to all streamtubes. Must be 

789 non-negative. 

790 longitudinal_dispersivity : float 

791 Longitudinal dispersivity alpha_L [m] (microdispersion), applied to all streamtubes. Must be non-negative. 

792 retardation_factor : float, optional 

793 Retardation factor (default 1.0). 

794 flow_out : array-like or None, optional 

795 Extraction flow rate [m³/day] on the output grid. See 

796 :func:`infiltration_to_extraction`. Default None. 

797 spinup : {"constant"} | None, optional 

798 See :func:`infiltration_to_extraction`. Default ``"constant"``. 

799 saturation_threshold : float, optional 

800 See :func:`infiltration_to_extraction`. Default 7.0. 

801 

802 Returns 

803 ------- 

804 numpy.ndarray 

805 Bin-averaged Kreft-Zuber flux concentration ``C_F`` in the extracted water. 

806 Length ``len(cout_tedges) - 1``. NaN where no infiltration data has broken through. 

807 

808 See Also 

809 -------- 

810 infiltration_to_extraction : Transport with an explicit pore volume distribution. 

811 gamma_extraction_to_infiltration : Reverse operation. 

812 gwtransport.gamma.bins : Create gamma distribution bins. 

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

814 """ 

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

816 return infiltration_to_extraction( 

817 cin=cin, 

818 flow=flow, 

819 tedges=tedges, 

820 cout_tedges=cout_tedges, 

821 aquifer_pore_volumes=bins["expected_values"], 

822 streamline_length=streamline_length, 

823 molecular_diffusivity=molecular_diffusivity, 

824 longitudinal_dispersivity=longitudinal_dispersivity, 

825 retardation_factor=retardation_factor, 

826 flow_out=flow_out, 

827 spinup=spinup, 

828 saturation_threshold=saturation_threshold, 

829 ) 

830 

831 

832def gamma_extraction_to_infiltration( 

833 *, 

834 cout: npt.ArrayLike, 

835 flow: npt.ArrayLike, 

836 tedges: pd.DatetimeIndex, 

837 cout_tedges: pd.DatetimeIndex, 

838 mean: float | None = None, 

839 std: float | None = None, 

840 loc: float = 0.0, 

841 alpha: float | None = None, 

842 beta: float | None = None, 

843 n_bins: int = 100, 

844 streamline_length: float, 

845 molecular_diffusivity: float, 

846 longitudinal_dispersivity: float, 

847 retardation_factor: float = 1.0, 

848 regularization_strength: float = 1e-10, 

849 flow_out: npt.ArrayLike | None = None, 

850 spinup: str | None = "constant", 

851 saturation_threshold: float = _DEFAULT_SATURATION_THRESHOLD, 

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

853 """Reconstruct infiltration concentration for a gamma-distributed pore volume distribution. 

854 

855 Convenience wrapper around :func:`extraction_to_infiltration` that discretizes a (shifted) 

856 gamma aquifer pore-volume distribution into ``n_bins`` equal-probability streamtubes. Provide 

857 either (mean, std) or (alpha, beta); ``loc`` defaults to 0. Fast approximate banded deconvolution 

858 (see :func:`extraction_to_infiltration`). 

859 

860 Parameters 

861 ---------- 

862 cout : array-like 

863 Concentration of the compound in extracted water. 

864 flow : array-like 

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

866 tedges : pandas.DatetimeIndex 

867 Time edges for cin (output) and flow data. Length ``len(flow) + 1``. 

868 cout_tedges : pandas.DatetimeIndex 

869 Time edges for cout data bins. Length ``len(cout) + 1``. 

870 mean, std : float, optional