Coverage for src/gwtransport/radial_asr.py: 88%
160 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1r"""Exact radial advection-dispersion transport for a single well (push-pull / ASR).
3Water is injected in an infinite aquifer at a single fully-penetrating well and later recovered at
4the same well under a signed flow schedule (push-pull / ASR). Transport is radial advection with
5microdispersion, molecular diffusion, and linear sorption; the spread of velocities across the well
6screen provides macrodispersion. Forward and backward modeling are supported.
8Computes the extracted flux concentration ``cout`` at a single fully-penetrating well driven by an
9arbitrary signed flow schedule (positive = injection, negative = extraction, zero = rest) and an
10arbitrary injected concentration ``cin``. The physics is the exact radial advection-dispersion of the
11radial ASR knowledge base: volume coordinate ``V(r) = pi b n (r^2 - r_w^2)``, Scheidegger
12velocity-dependent dispersion ``D = alpha_L |u| + D_m`` (microdispersion ``alpha_L |u|`` plus molecular diffusion ``D_m``), Kreft-Zuber flux boundary conditions, and
13the exact per-phase kernels (Airy for ``D_m = 0``; the log-derivative Riccati ODE for ``D_m > 0``).
14Nothing is reduced to a Gaussian; the exact
15non-Gaussian breakthrough (with the correct skewness) is carried.
17The forward map is **grid-free** end to end -- no PDE is discretized, so none of the finite-volume
18artefacts appear. A single inject-then-extract cycle with no intervening rest uses the closed-form echo operator
19(``gwtransport._radial_asr_compose``, KB Sec. 10a) -- exact for arbitrary within-phase variable flow,
20with the exact temporal moments. Any other signed-flow schedule (more reversals / multi-cycle ASR, or a
21single cycle with a rest under nonzero ``D_m``) uses the reused-propagator-matrix engine
22(``gwtransport._radial_asr_reuse``, KB addendum Sec. A1-A7), which composes the exact per-phase kernels
23(Airy / Riccati / Bessel) through the interior two-point Green's functions. Each per-reversal field
24hand-off ``f_out = P @ f`` is a bounded linear operator; its matrix ``P`` is built once per distinct
25``(direction, phase volume)`` from a single batched de Hoog inversion and reused at every recurrence, so
26the special-function + inversion cost is ``O(distinct phase volumes)`` rather than ``O(reversals)``. It is
27bit-equivalent, to the de Hoog floor, to the per-reversal grid-free composition. Molecular diffusion during
28pumping (the ``D_m > 0`` Whittaker kernel) is evaluated through the log-derivative Riccati ODE
29(``gwtransport._radial_asr_kernels.resolvent_riccati``) -- exact to the de Hoog inversion floor at any
30``A_0/D_m``, with no special-function precision cap, and reducing continuously to the Airy branch as
31``D_m -> 0``. During a **rest** (``Q = 0``) advection and microdispersion vanish and molecular
32diffusion acts alone on the wall-clock clock; it is carried exactly by the order-0 modified Bessel
33pure-diffusion kernel, the dominant mixing for seasonal storage / ATES. The only
34numerical steps are Gauss-Legendre quadrature and de Hoog Laplace inversion of exact special-function
35kernels. An independent finite-volume solve of the same PDE (``tests/src/_radial_asr_fv_oracle.py``,
36KB Sec. 9) is used only as a test oracle. The propagator matrices are assembled on the Bromwich
37contour (``Re s > 0``), where the field hand-off is well-conditioned at any Peclet. The engine is chosen
38automatically; cycles are expressed through the flow sign pattern, not
39an argument.
41The reported ``cout`` is the flow-weighted average over each output bin -- defined on extraction bins
42(``flow < 0``) and ``NaN`` on injection / rest bins (nothing is recovered there).
44Macrodispersion within the well screen
45--------------------------------------
46The well screen has a **known** height; macrodispersion is the spread of arrival times caused by
47*velocity heterogeneity across the screen*. It is modelled as parallel streamtubes (``pore_heights``):
48each streamtube is an independent radial cell carrying the full flow, with an effective pore height
49that sets its velocity, and the output is the weight-averaged breakthrough. A streamtube of effective
50height ``b`` has velocity ``proportional to 1/b`` (its pore volume to radius ``r`` is
51``pi b n (r^2 - r_w^2)``), so smaller ``b`` means faster breakthrough.
52:func:`gamma_infiltration_to_extraction` builds this ensemble from a gamma distribution of the layer
53**velocity** within the fixed screen height (see that function); the mean velocity is set by the screen
54height and the spread by a velocity coefficient of variation. The spread is a within-screen velocity
55distribution -- velocity heterogeneity across the well screen -- not an aquifer pore-volume distribution.
57Regional background flow (drift)
58--------------------------------
59With a steady uniform regional Darcy flux ``regional_flux`` (``U``, drift seepage ``v_d = U/n``) the well
60field is superimposed on a regional gradient, so the stored bubble drifts and recovery degrades. The
61radial symmetry is broken and the transport is solved by an **azimuthal Fourier-mode** expansion
62``c(r, theta) = sum_m c_m(r) e^{i m theta}`` (``m = 0`` is the radial engine; drift couples ``m`` to
63``m +- 1``), composed through the same per-phase interior Green's functions
64(``gwtransport._radial_asr_drift_kernels``). ``regional_flux = 0`` (default) dispatches to the radial path
65bit-for-bit. The engine is for the **slow-drift** envelope -- the plume (including its rest-phase drift
66displacement) must stay well inside the stagnation radius ``r_s = |A_0|/|v_d|`` (else a ``ValueError``).
67Rest phases (``flow == 0``) are propagated by the exact free-space drift kernel (translate + anisotropic
68spread). The drift-induced recovery loss is validated against an independent 2-D finite-volume oracle.
70Available functions:
72- :func:`infiltration_to_extraction` -- forward transport (cin -> cout).
73- :func:`extraction_to_infiltration` -- inverse via Tikhonov regularization (cout -> cin).
74- :func:`gamma_infiltration_to_extraction` -- gamma-distributed screen velocity (forward).
75- :func:`gamma_extraction_to_infiltration` -- same, inverse.
77References
78----------
79The references below give the published closed-form solutions for the **single-phase** radial *injection*
80problem (steady divergent flow from one well) -- the per-phase forward kernel this module composes. The
81convergent-extraction dual (KB Sec. 7) and the multi-cycle push-pull / ASR composition across flow
82reversals are built on top of those kernels here and are not in the single-injection references. All
83share the assumptions used here: a single fully-penetrating well in a homogeneous medium with steady
84divergent flow ``v = Q / (2 pi b n r)``, plus retardation.
86The ``D_m = 0`` kernel (velocity-proportional microdispersion ``D = alpha_L |u|``, Airy functions)
87is the classical radial-dispersion problem: Tang & Babu (1979) under a Dirichlet (resident-concentration)
88well boundary, and Chen (1987) under the Cauchy / third-type (flux) boundary used here -- explicitly the
89Kreft-Zuber flux concentration, with transfer function ``Ai(Y) / [Ai(Y0)/2 - p^(1/3) Ai'(Y0)]`` equal to
90the flux operator this module evaluates. The ``D_m > 0`` kernel (``D = alpha_L |u| + D_m``, Kummer /
91confluent-hypergeometric functions) under the same flux boundary, with retardation, is Aichi & Akitaya
92(2018) -- whose well operator ``U(a,b) + 2a U(a+1,b+1)`` is this module's Whittaker flux boundary; they
93record the ``D_m -> 0`` reduction to Chen (1987) as an open problem, which this module performs
94continuously -- the log-derivative Riccati kernel reduces smoothly to the Airy branch as ``D_m -> 0``.
95The ``alpha_L = 0`` limit (constant diffusion, drift-dominated radial transport, Whittaker equation) is
96Akanji & Falade (2019). Each is an injection-only solution; none treats extraction or multi-cycle push-pull.
98Kreft, A., & Zuber, A. (1978). On the physical meaning of the dispersion equation and its solutions
99for different initial and boundary conditions. Chemical Engineering Science, 33(11), 1471-1480.
101Tang, D. H., & Babu, D. K. (1979). Analytical solution of a velocity dependent dispersion problem.
102Water Resources Research, 15(6), 1471-1478.
104Chen, C.-S. (1987). Analytical solutions for radial dispersion with Cauchy boundary at injection well.
105Water Resources Research, 23(7), 1217-1224.
107Aichi, M., & Akitaya, K. (2018). Analytical solution for a radial advection-dispersion equation
108including both mechanical dispersion and molecular diffusion for a steady-state flow field in a
109horizontal aquifer caused by a constant rate injection from a well. Hydrological Research Letters,
11012(3), 23-27.
112Akanji, L. T., & Falade, G. K. (2019). Closed-form solution of radial transport of tracers in porous
113media influenced by linear drift. Energies, 12(1), 29.
115This file is part of gwtransport which is released under AGPL-3.0 license.
116See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
117"""
119import numpy as np
120import numpy.typing as npt
121import pandas as pd
123from gwtransport import gamma
124from gwtransport._radial_asr_compose import single_cycle_echo_matrix
125from gwtransport._radial_asr_drift_kernels import _RS_FRAC, block_cout_deviation
126from gwtransport._radial_asr_reuse import cout_deviation
127from gwtransport._time import dt_to_days
128from gwtransport._validation import _validate_retardation_factor
131def _is_single_cycle(flow: npt.NDArray[np.floating]) -> bool:
132 """Return True if the schedule is a single injection block followed by a single extraction block.
134 Such schedules (one flow reversal, injection first) use the exact closed-form echo operator; any
135 other signed-flow pattern (more reversals, extraction first) uses the reused-propagator-matrix engine.
137 Returns
138 -------
139 bool
140 Whether ``flow`` is a single inject-then-extract cycle.
141 """
142 signs = np.sign(flow[flow != 0.0])
143 n_changes = int(np.sum(np.diff(signs) != 0)) if signs.size else 0
144 return n_changes <= 1 and (signs.size == 0 or signs[0] > 0)
147def _validate(
148 *,
149 cin_or_cout: npt.NDArray[np.floating],
150 flow: npt.NDArray[np.floating],
151 tedges: pd.DatetimeIndex,
152 cout_tedges: pd.DatetimeIndex,
153 pore_heights: npt.NDArray[np.floating],
154 porosity: float,
155 well_radius: float,
156 longitudinal_dispersivity: float,
157 molecular_diffusivity: float,
158 retardation_factor: float,
159 weights: npt.NDArray[np.floating] | None,
160 regional_flux: float,
161 n_modes: int | None,
162) -> None:
163 """Validate inputs for the radial single-well transport functions (signed flow is allowed).
165 Raises
166 ------
167 ValueError
168 On inconsistent lengths, non-positive geometry, out-of-range porosity, negative dispersion,
169 ``retardation_factor < 1``, mismatched ``weights``, NaN in ``flow``, a non-finite ``regional_flux``,
170 an ``n_modes < 1``, or a ``cout_tedges`` that differs from ``tedges`` (a distinct output grid is not
171 yet supported).
172 """
173 if len(tedges) != len(flow) + 1:
174 msg = "tedges must have one more element than flow"
175 raise ValueError(msg)
176 if len(cin_or_cout) != len(flow):
177 msg = "cin/cout must have the same length as flow"
178 raise ValueError(msg)
179 if not tedges.equals(cout_tedges):
180 msg = "cout_tedges must equal tedges (a distinct output grid is not yet supported)"
181 raise ValueError(msg)
182 if np.any(np.isnan(flow)):
183 msg = "flow contains NaN values, which are not allowed"
184 raise ValueError(msg)
185 if np.any(pore_heights <= 0.0):
186 msg = "pore_heights must be positive"
187 raise ValueError(msg)
188 if not 0.0 < porosity <= 1.0:
189 msg = "porosity must be in (0, 1]"
190 raise ValueError(msg)
191 if well_radius <= 0.0:
192 msg = "well_radius must be positive"
193 raise ValueError(msg)
194 if longitudinal_dispersivity <= 0.0:
195 msg = "longitudinal_dispersivity must be positive (the dispersion kernel requires alpha_L > 0)"
196 raise ValueError(msg)
197 if molecular_diffusivity < 0.0:
198 msg = "molecular_diffusivity must be non-negative"
199 raise ValueError(msg)
200 _validate_retardation_factor(retardation_factor)
201 if weights is not None and len(weights) != len(pore_heights):
202 msg = "weights must have the same length as pore_heights"
203 raise ValueError(msg)
204 if not np.isfinite(regional_flux):
205 msg = "regional_flux must be finite"
206 raise ValueError(msg)
207 if n_modes is not None and n_modes < 1:
208 msg = "n_modes must be >= 1"
209 raise ValueError(msg)
212def _echo_operator(
213 *,
214 flow: npt.NDArray[np.floating],
215 tedges: pd.DatetimeIndex,
216 c_geos: npt.NDArray[np.floating],
217 well_radius: float,
218 longitudinal_dispersivity: float,
219 molecular_diffusivity: float,
220 retardation_factor: float,
221 weights: npt.NDArray[np.floating],
222 n_quad: int,
223) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.bool_], npt.NDArray[np.bool_]]:
224 """Weight-averaged single-cycle echo operator ``W`` (``cout' = W @ cin'_inj``) over the streamtubes.
226 Builds the closed-form echo matrix per streamtube (geometry constant ``c_geo = pi b n``) and
227 averages by ``weights``. Used by both the forward (``cout = W @ cin``) and the reverse (Tikhonov).
229 Returns
230 -------
231 w_ens : ndarray, shape (n_ext, n_inj)
232 Weight-averaged echo operator.
233 inj_mask, ext_mask : ndarray of bool
234 Injection (``flow > 0``) and extraction (``flow < 0``) bin masks.
235 """
236 inj_mask, ext_mask = flow > 0.0, flow < 0.0
237 dt = dt_to_days(tedges)
238 inj_vol = np.concatenate(([0.0], np.cumsum((flow * dt)[inj_mask]))) # 0 .. S_inj
239 ext_vol = np.concatenate(([0.0], np.cumsum((-flow * dt)[ext_mask]))) # 0 .. T_end
240 inj_flow_scale = float(np.mean(flow[inj_mask])) if np.any(inj_mask) else 1.0
241 ext_flow_scale = float(np.mean(-flow[ext_mask])) if np.any(ext_mask) else 1.0
242 w_ens = np.zeros((int(np.sum(ext_mask)), int(np.sum(inj_mask))))
243 for c_geo, w_i in zip(c_geos, weights, strict=True):
244 w_ens += w_i * single_cycle_echo_matrix(
245 inj_volume_edges=inj_vol,
246 ext_volume_edges=ext_vol,
247 c_geo=c_geo,
248 r_w=well_radius,
249 alpha_l=longitudinal_dispersivity,
250 inj_flow_scale=inj_flow_scale,
251 ext_flow_scale=ext_flow_scale,
252 retardation_factor=retardation_factor,
253 molecular_diffusivity=molecular_diffusivity,
254 n_quad=n_quad,
255 )
256 return w_ens / np.sum(weights), inj_mask, ext_mask
259def _reuse_ensemble(
260 cin_deviation: npt.NDArray[np.floating],
261 *,
262 flow: npt.NDArray[np.floating],
263 dt_days: npt.NDArray[np.floating],
264 c_geos: npt.NDArray[np.floating],
265 well_radius: float,
266 longitudinal_dispersivity: float,
267 molecular_diffusivity: float,
268 retardation_factor: float,
269 weights: npt.NDArray[np.floating],
270 n_quad: int,
271) -> npt.NDArray[np.floating]:
272 """Weight-averaged multi-cycle extracted-flux deviation over the streamtubes.
274 Runs the reused-propagator-matrix multi-cycle engine once per streamtube (geometry constant
275 ``c_geo = pi b n``) and averages by ``weights``. ``cin_deviation`` may be ``(n,)`` or ``(n, k)`` -- a
276 column batch is transported through one engine pass per streamtube (the per-phase propagator / source /
277 readout matrices are cin-independent, so they are built once and applied to every column). Used by the
278 forward (with ``cin``) and the reverse (with the unit-pulse column batch).
280 Returns
281 -------
282 ndarray, shape (n,) or (n, k)
283 Weight-averaged extracted-flux deviation on extraction bins (matching ``cin_deviation``), ``0``
284 elsewhere.
285 """
286 # cout_deviation returns NaN on injection / rest bins by contract (structural, expected). Only the
287 # extraction bins are accumulated, so those structural NaNs are dropped without a blanket nan_to_num --
288 # a genuine NaN on an EXTRACTION bin (a real numerical failure) then propagates and surfaces instead of
289 # silently reading as a physical zero.
290 ext_mask = flow < 0.0
291 acc = np.zeros(np.shape(cin_deviation))
292 for c_geo, w_i in zip(c_geos, weights, strict=True):
293 dev = cout_deviation(
294 cin_deviation=cin_deviation,
295 flow=flow,
296 dt_days=dt_days,
297 c_geo=c_geo,
298 r_w=well_radius,
299 alpha_l=longitudinal_dispersivity,
300 molecular_diffusivity=molecular_diffusivity,
301 retardation_factor=retardation_factor,
302 n_quad=n_quad,
303 )
304 acc[ext_mask] += w_i * dev[ext_mask]
305 return acc / np.sum(weights)
308def _auto_n_modes(
309 flow: npt.NDArray[np.floating],
310 dt_days: npt.NDArray[np.floating],
311 c_geo: float,
312 well_radius: float,
313 longitudinal_dispersivity: float,
314 v_d: float,
315) -> int:
316 r"""Azimuthal truncation ``M`` sized from the drift ratio ``eps`` and the rest-phase displacement.
318 The pumping-phase mode amplitudes decay geometrically, ``|c_m| ~ eps^|m|`` with
319 ``eps = v_d R_b / A_0``, so keeping modes ``-M .. M`` truncates the azimuthal field at
320 ``O(eps^{M+1})``; ``M`` is chosen so that tail is below ``~5e-3``. An interior rest phase
321 additionally translates the plume by ``delta = v_d t_rest`` (``R = 1``, conservative; idle bins
322 before the first or after the last pumping do not move the field), populating harmonics up
323 to ``~ delta / width`` (``width`` the radial breakthrough std) -- the second bound. The result is
324 clamped to ``[2, 8]`` (the slow-drift envelope -- beyond ``eps ~ 0.6`` the far-field escape this
325 engine does not model dominates anyway); the rest kernel's honest spectral-tail guard raises if a long
326 rest still outruns the clamp (pass ``n_modes`` explicitly then). ``A_0`` uses the **smallest** pumping
327 magnitude (the worst-case largest ``eps``, consistent with the stagnation-radius envelope guard),
328 ``R_b`` the peak net injected radius. This function is only reached for nonzero drift.
330 Returns
331 -------
332 int
333 Azimuthal truncation ``M``.
334 """
335 pumping = np.abs(flow[flow != 0.0])
336 if pumping.size == 0: # all-rest schedule: nothing pumps; the engine returns all-NaN downstream
337 return 2
338 a0 = float(np.min(pumping)) / (2.0 * c_geo)
339 net_volume = np.concatenate(([0.0], np.cumsum(flow * dt_days)))
340 peak_volume = max(float(net_volume.max()), 0.0)
341 r_b = np.sqrt(well_radius**2 + peak_volume / c_geo)
342 nz = np.flatnonzero(flow != 0.0)
343 interior = slice(nz[0], nz[-1] + 1) # leading/trailing idle bins do not move the field
344 delta = abs(v_d) * float(np.sum(dt_days[interior][flow[interior] == 0.0]))
345 eps = min(abs(v_d) * (r_b + delta) / abs(a0), _RS_FRAC)
346 m_eps = int(np.ceil(np.log(5e-3) / np.log(eps))) if eps > 0.0 else 2
347 width = np.sqrt(longitudinal_dispersivity * r_b + longitudinal_dispersivity**2)
348 m_shift = int(np.ceil(delta / width)) + 2 if delta > 0.0 else 2
349 return int(np.clip(max(m_eps, m_shift), 2, 8))
352def _block_ensemble(
353 cin_deviation: npt.NDArray[np.floating],
354 *,
355 flow: npt.NDArray[np.floating],
356 dt_days: npt.NDArray[np.floating],
357 c_geos: npt.NDArray[np.floating],
358 porosity: float,
359 well_radius: float,
360 longitudinal_dispersivity: float,
361 molecular_diffusivity: float,
362 retardation_factor: float,
363 regional_flux: float,
364 n_modes: int | None,
365 weights: npt.NDArray[np.floating],
366 n_quad: int,
367) -> npt.NDArray[np.floating]:
368 """Weight-averaged multi-cycle extracted-flux deviation with regional drift over the streamtubes.
370 Runs the azimuthal-mode block engine (:func:`gwtransport._radial_asr_drift_kernels.block_cout_deviation`)
371 once per streamtube (geometry constant ``c_geo = pi b n``) and averages by ``weights``. The drift
372 seepage ``v_d = U / n`` is the same for every streamtube (a regional Darcy flux through the porosity);
373 only the radial strength ``A_0 ~ 1/c_geo`` varies, so faster (thinner) streamtubes see a smaller drift
374 ratio. ``n_modes`` is auto-sized per streamtube from its drift ratio when not given. ``cin_deviation``
375 may be ``(n,)`` or ``(n, k)`` -- a column batch is transported through one engine pass per streamtube
376 (used by the reverse operator build).
378 Returns
379 -------
380 ndarray, shape (n,) or (n, k)
381 Weight-averaged extracted-flux deviation on extraction bins (matching ``cin_deviation``), ``0``
382 elsewhere.
383 """
384 v_d = regional_flux / porosity
385 acc = np.zeros(np.shape(cin_deviation))
386 for c_geo, w_i in zip(c_geos, weights, strict=True):
387 m = (
388 n_modes
389 if n_modes is not None
390 else _auto_n_modes(flow, dt_days, c_geo, well_radius, longitudinal_dispersivity, v_d)
391 )
392 acc += w_i * np.nan_to_num(
393 block_cout_deviation(
394 cin_deviation=cin_deviation,
395 flow=flow,
396 dt_days=dt_days,
397 c_geo=c_geo,
398 r_w=well_radius,
399 alpha_l=longitudinal_dispersivity,
400 v_d=v_d,
401 molecular_diffusivity=molecular_diffusivity,
402 retardation_factor=retardation_factor,
403 n_modes=m,
404 n_quad=n_quad,
405 )
406 )
407 return acc / np.sum(weights)
410def infiltration_to_extraction(
411 *,
412 cin: npt.ArrayLike,
413 flow: npt.ArrayLike,
414 tedges: pd.DatetimeIndex,
415 cout_tedges: pd.DatetimeIndex,
416 pore_heights: npt.ArrayLike,
417 porosity: float,
418 well_radius: float,
419 longitudinal_dispersivity: float,
420 molecular_diffusivity: float = 0.0,
421 retardation_factor: float = 1.0,
422 weights: npt.ArrayLike | None = None,
423 background: float = 0.0,
424 regional_flux: float = 0.0,
425 n_modes: int | None = None,
426 n_quad: int = 240,
427) -> npt.NDArray[np.floating]:
428 """Compute the extracted flux concentration at a radial well for a signed flow schedule.
430 Parameters
431 ----------
432 cin : array-like, shape (n,)
433 Injected concentration per time bin (used only on injection bins, ``flow > 0``).
434 flow : array-like, shape (n,)
435 Signed flow per time bin [m^3/day]: ``> 0`` injection, ``< 0`` extraction, ``0`` rest.
436 tedges : DatetimeIndex
437 Time bin edges (``n + 1`` for ``n`` bins).
438 cout_tedges : DatetimeIndex
439 Output time bin edges; must equal ``tedges``. Output is NaN on injection / rest bins.
440 pore_heights : array-like
441 Effective streamtube pore height(s) ``b`` [m] -- a scalar (one homogeneous screen) or an array
442 of streamtube heights for the velocity-heterogeneity macrodispersion ensemble (each streamtube
443 carries the full flow; smaller ``b`` = faster). See the module docstring and
444 :func:`gamma_infiltration_to_extraction`.
445 porosity : float
446 Porosity ``n`` [-].
447 well_radius : float
448 Well (screen) radius ``r_w`` [m].
449 longitudinal_dispersivity : float
450 Longitudinal dispersivity ``alpha_L`` [m].
451 molecular_diffusivity : float, optional
452 Molecular diffusivity ``D_m`` [m^2/day]. Default 0. ``D_m = 0`` uses the vectorized Airy branch;
453 ``D_m > 0`` uses the log-derivative Riccati kernel -- exact to the de Hoog floor at any ``A_0/D_m``
454 with no precision cap, reducing continuously to the Airy branch as ``D_m -> 0``.
455 retardation_factor : float, optional
456 Linear retardation ``R >= 1``. Default 1.
457 weights : array-like, optional
458 Per-streamtube averaging weights (same length as ``pore_heights``). Default equal weights.
459 background : float, optional
460 Ambient aquifer concentration ``c_bg``. The deviation ``cin - c_bg`` is transported and
461 ``c_bg`` is added back; constant ``cin = c_bg`` returns ``cout = c_bg``. Default 0.
462 regional_flux : float, optional
463 Steady uniform regional background Darcy flux ``U`` [m/day] in ``+x`` (drift seepage
464 ``v_d = U / n``). ``0`` (default) reproduces the radial-symmetric engine bit-for-bit. A nonzero
465 value engages the azimuthal-mode block engine, which captures the drift-induced recovery loss
466 (the down-gradient plume is partly swept past the well). The slow-drift envelope requires the
467 plume -- including its rest-phase drift displacement -- to stay well inside the stagnation radius
468 ``r_s = |A_0| / |v_d|`` (a ``ValueError`` is raised otherwise). Rest phases (``flow == 0``) are
469 propagated by the exact free-space drift kernel (translation ``v_d t / R`` plus anisotropic
470 Gaussian spread, with a Neumann-image closure at the shut well face). See
471 :ref:`concept-drift-envelope` for a worked multi-year feasibility table.
472 n_modes : int, optional
473 Azimuthal truncation ``M`` for the drift engine (keeps modes ``-M .. M``). Default ``None``
474 auto-sizes ``M`` from the plume-front drift ratio ``eps = v_d R_b / A_0`` and the rest-phase
475 displacement (clamped to ``[2, 8]``). Ignored when ``regional_flux == 0``.
476 n_quad : int, optional
477 Gauss-Legendre node count for the resident-profile superposition. Default 240.
479 Returns
480 -------
481 ndarray, shape (n,)
482 Extracted flux concentration; NaN on injection and rest bins.
483 """
484 cin = np.asarray(cin, dtype=float)
485 flow = np.asarray(flow, dtype=float)
486 pore_heights = np.atleast_1d(np.asarray(pore_heights, dtype=float))
487 weights_arr = np.ones(len(pore_heights)) if weights is None else np.atleast_1d(np.asarray(weights, dtype=float))
488 _validate(
489 cin_or_cout=cin,
490 flow=flow,
491 tedges=tedges,
492 cout_tedges=cout_tedges,
493 pore_heights=pore_heights,
494 porosity=porosity,
495 well_radius=well_radius,
496 longitudinal_dispersivity=longitudinal_dispersivity,
497 molecular_diffusivity=molecular_diffusivity,
498 retardation_factor=retardation_factor,
499 weights=None if weights is None else weights_arr,
500 regional_flux=regional_flux,
501 n_modes=n_modes,
502 )
503 c_geos = np.pi * pore_heights * porosity
504 cout = np.full(len(flow), np.nan)
505 if regional_flux != 0.0:
506 # Steady regional drift breaks radial symmetry: the azimuthal-mode block engine carries it.
507 ext_mask = flow < 0.0
508 cout_dev = _block_ensemble(
509 cin - background,
510 flow=flow,
511 dt_days=dt_to_days(tedges),
512 c_geos=c_geos,
513 porosity=porosity,
514 well_radius=well_radius,
515 longitudinal_dispersivity=longitudinal_dispersivity,
516 molecular_diffusivity=molecular_diffusivity,
517 retardation_factor=retardation_factor,
518 regional_flux=regional_flux,
519 n_modes=n_modes,
520 weights=weights_arr,
521 n_quad=n_quad,
522 )
523 cout[ext_mask] = background + cout_dev[ext_mask]
524 return cout
525 # A rest phase combined with molecular diffusion (seasonal storage / ATES) cannot use the
526 # flushed-volume echo operator, which is blind to a rest's wall-clock diffusion; route it to the
527 # reuse engine (which propagates the rest with the Bessel pure-diffusion kernel).
528 use_echo = _is_single_cycle(flow) and not (molecular_diffusivity > 0.0 and np.any(flow == 0.0))
529 if use_echo:
530 w_ens, inj_mask, ext_mask = _echo_operator(
531 flow=flow,
532 tedges=tedges,
533 c_geos=c_geos,
534 well_radius=well_radius,
535 longitudinal_dispersivity=longitudinal_dispersivity,
536 molecular_diffusivity=molecular_diffusivity,
537 retardation_factor=retardation_factor,
538 weights=weights_arr,
539 n_quad=n_quad,
540 )
541 cout[ext_mask] = background + w_ens @ (cin[inj_mask] - background)
542 return cout
543 ext_mask = flow < 0.0
544 cout_dev = _reuse_ensemble(
545 cin - background,
546 flow=flow,
547 dt_days=dt_to_days(tedges),
548 c_geos=c_geos,
549 well_radius=well_radius,
550 longitudinal_dispersivity=longitudinal_dispersivity,
551 molecular_diffusivity=molecular_diffusivity,
552 retardation_factor=retardation_factor,
553 weights=weights_arr,
554 n_quad=n_quad,
555 )
556 cout[ext_mask] = background + cout_dev[ext_mask]
557 return cout
560def extraction_to_infiltration(
561 *,
562 cout: npt.ArrayLike,
563 flow: npt.ArrayLike,
564 tedges: pd.DatetimeIndex,
565 cout_tedges: pd.DatetimeIndex,
566 pore_heights: npt.ArrayLike,
567 porosity: float,
568 well_radius: float,
569 longitudinal_dispersivity: float,
570 molecular_diffusivity: float = 0.0,
571 retardation_factor: float = 1.0,
572 weights: npt.ArrayLike | None = None,
573 background: float = 0.0,
574 regional_flux: float = 0.0,
575 n_modes: int | None = None,
576 regularization_strength: float = 1e-10,
577 n_quad: int = 240,
578) -> npt.NDArray[np.floating]:
579 """Recover the injected concentration from extracted-water measurements (Tikhonov inverse).
581 Inverts the forward operator built by :func:`infiltration_to_extraction`. Returns the injected
582 concentration on injection bins (NaN on extraction / rest bins).
584 Parameters
585 ----------
586 cout : array-like, shape (n,)
587 Measured extracted concentration (used on extraction bins, ``flow < 0``).
588 flow, tedges, cout_tedges, pore_heights, porosity, well_radius, longitudinal_dispersivity
589 As in :func:`infiltration_to_extraction`.
590 molecular_diffusivity, retardation_factor, weights, background, regional_flux, n_modes, n_quad
591 As in :func:`infiltration_to_extraction`.
592 regularization_strength : float, optional
593 Tikhonov parameter. Default ``1e-10``.
595 Returns
596 -------
597 ndarray, shape (n,)
598 Recovered injected concentration; NaN on extraction / rest bins.
600 Raises
601 ------
602 ValueError
603 If ``cout`` contains NaN on any extraction bin (``flow < 0``), which would poison the
604 least-squares solve. Structural NaN on injection / rest bins is allowed.
605 """
606 cout = np.asarray(cout, dtype=float)
607 flow = np.asarray(flow, dtype=float)
608 pore_heights = np.atleast_1d(np.asarray(pore_heights, dtype=float))
609 weights_arr = np.ones(len(pore_heights)) if weights is None else np.atleast_1d(np.asarray(weights, dtype=float))
610 _validate(
611 cin_or_cout=cout,
612 flow=flow,
613 tedges=tedges,
614 cout_tedges=cout_tedges,
615 pore_heights=pore_heights,
616 porosity=porosity,
617 well_radius=well_radius,
618 longitudinal_dispersivity=longitudinal_dispersivity,
619 molecular_diffusivity=molecular_diffusivity,
620 retardation_factor=retardation_factor,
621 weights=None if weights is None else weights_arr,
622 regional_flux=regional_flux,
623 n_modes=n_modes,
624 )
625 # A NaN measurement on an extraction bin (the only bins the inverse reads) would poison the whole
626 # least-squares solve into an all-NaN cin; raise instead, as the advection / diffusion inverses do.
627 # Structural NaN on injection / rest bins is allowed (those bins are ignored by the inverse).
628 if np.any(np.isnan(cout[flow < 0.0])):
629 msg = "cout contains NaN values on extraction bins, which are not allowed"
630 raise ValueError(msg)
631 c_geos = np.pi * pore_heights * porosity
632 # A rest phase with molecular diffusion routes to the reuse engine (see the forward function); regional
633 # drift never uses the radial echo operator (it breaks the azimuthal symmetry the echo relies on).
634 use_echo = (
635 regional_flux == 0.0 and _is_single_cycle(flow) and not (molecular_diffusivity > 0.0 and np.any(flow == 0.0))
636 )
637 if use_echo:
638 w_ens, inj_mask, ext_mask = _echo_operator(
639 flow=flow,
640 tedges=tedges,
641 c_geos=c_geos,
642 well_radius=well_radius,
643 longitudinal_dispersivity=longitudinal_dispersivity,
644 molecular_diffusivity=molecular_diffusivity,
645 retardation_factor=retardation_factor,
646 weights=weights_arr,
647 n_quad=n_quad,
648 )
649 else:
650 # Build the dense forward operator W whose columns are the unit-injection-pulse responses; the reverse
651 # cannot reuse the cheap single-solve forward path. The per-phase propagator / source / readout matrices
652 # are cin-independent (flow + geometry only), so the whole unit-pulse batch is transported in ONE engine
653 # pass -- the matrices are built once and applied to every column. The block (drift) engine carries
654 # the regional drift; the scalar reuse engine the drift-free case.
655 inj_mask, ext_mask = flow > 0.0, flow < 0.0
656 inj_idx = np.flatnonzero(inj_mask)
657 dt_days = dt_to_days(tedges)
658 pulses = np.zeros((len(flow), len(inj_idx)))
659 pulses[inj_idx, np.arange(len(inj_idx))] = 1.0
660 if regional_flux != 0.0:
661 cols = _block_ensemble(
662 pulses,
663 flow=flow,
664 dt_days=dt_days,
665 c_geos=c_geos,
666 porosity=porosity,
667 well_radius=well_radius,
668 longitudinal_dispersivity=longitudinal_dispersivity,
669 molecular_diffusivity=molecular_diffusivity,
670 retardation_factor=retardation_factor,
671 regional_flux=regional_flux,
672 n_modes=n_modes,
673 weights=weights_arr,
674 n_quad=n_quad,
675 )
676 else:
677 cols = _reuse_ensemble(
678 pulses,
679 flow=flow,
680 dt_days=dt_days,
681 c_geos=c_geos,
682 well_radius=well_radius,
683 longitudinal_dispersivity=longitudinal_dispersivity,
684 molecular_diffusivity=molecular_diffusivity,
685 retardation_factor=retardation_factor,
686 weights=weights_arr,
687 n_quad=n_quad,
688 )
689 w_ens = cols[ext_mask, :]
690 # Tikhonov least-squares min ||W x - (cout-bg)||^2 + lambda ||x||^2 via the stable augmented
691 # system [W; sqrt(lambda) I] x = [cout-bg; 0]. The echo / reuse operator has column sums ~1
692 # (mass conservation per injection bin) and overdetermined rows, so a direct Tikhonov fit is used.
693 n_inj = w_ens.shape[1]
694 augmented = np.vstack([w_ens, np.sqrt(regularization_strength) * np.eye(n_inj)])
695 rhs = np.concatenate([cout[ext_mask] - background, np.zeros(n_inj)])
696 cin_dev = np.linalg.lstsq(augmented, rhs, rcond=None)[0]
697 cin = np.full(len(flow), np.nan)
698 cin[inj_mask] = background + cin_dev
699 return cin
702def gamma_infiltration_to_extraction(
703 *,
704 cin: npt.ArrayLike,
705 flow: npt.ArrayLike,
706 tedges: pd.DatetimeIndex,
707 cout_tedges: pd.DatetimeIndex,
708 porosity: float,
709 well_radius: float,
710 longitudinal_dispersivity: float,
711 screen_height: float,
712 velocity_cv: float,
713 n_bins: int = 100,
714 molecular_diffusivity: float = 0.0,
715 retardation_factor: float = 1.0,
716 background: float = 0.0,
717 regional_flux: float = 0.0,
718 n_modes: int | None = None,
719 n_quad: int = 240,
720) -> npt.NDArray[np.floating]:
721 """Radial transport with gamma-distributed screen velocity (within-screen macrodispersion).
723 The well screen has a **known** height ``screen_height``; macrodispersion is the spread of arrival
724 times from velocity heterogeneity across that fixed height. The layer velocity is gamma-distributed
725 with mean equal to the homogeneous value (a streamtube at the mean velocity has effective pore
726 height ``screen_height``) and coefficient of variation ``velocity_cv``. A streamtube with velocity
727 ratio ``rho`` (gamma, mean 1) has effective pore height ``screen_height / rho`` -- faster layers are
728 thinner and break through sooner. The gamma is discretized into ``n_bins`` equal-probability bins
729 (:func:`gwtransport.gamma.bins`) and averaged by probability mass via
730 :func:`infiltration_to_extraction`.
732 Parameters
733 ----------
734 screen_height : float
735 Known well-screen height ``H`` [m] (the fixed total; the mean streamtube velocity is set by it).
736 velocity_cv : float
737 Coefficient of variation of the layer velocity (the macrodispersion strength). ``0`` is a
738 homogeneous screen (a single streamtube, sharp breakthrough); typically ``< 1`` -- larger values
739 give a heavy slow-velocity tail.
740 n_bins : int, optional
741 Number of equal-probability velocity bins. Default 100.
742 cin, flow, tedges, cout_tedges, porosity, well_radius, longitudinal_dispersivity
743 As in :func:`infiltration_to_extraction`.
744 molecular_diffusivity, retardation_factor, background, regional_flux, n_modes, n_quad
745 As in :func:`infiltration_to_extraction`.
747 Returns
748 -------
749 ndarray, shape (n,)
750 Extracted flux concentration; NaN on injection / rest bins.
751 """
752 pore_heights, weights = _velocity_gamma_streamtubes(screen_height, velocity_cv, n_bins)
753 return infiltration_to_extraction(
754 cin=cin,
755 flow=flow,
756 tedges=tedges,
757 cout_tedges=cout_tedges,
758 pore_heights=pore_heights,
759 porosity=porosity,
760 well_radius=well_radius,
761 longitudinal_dispersivity=longitudinal_dispersivity,
762 molecular_diffusivity=molecular_diffusivity,
763 retardation_factor=retardation_factor,
764 weights=weights,
765 background=background,
766 regional_flux=regional_flux,
767 n_modes=n_modes,
768 n_quad=n_quad,
769 )
772def gamma_extraction_to_infiltration(
773 *,
774 cout: npt.ArrayLike,
775 flow: npt.ArrayLike,
776 tedges: pd.DatetimeIndex,
777 cout_tedges: pd.DatetimeIndex,
778 porosity: float,
779 well_radius: float,
780 longitudinal_dispersivity: float,
781 screen_height: float,
782 velocity_cv: float,
783 n_bins: int = 100,
784 molecular_diffusivity: float = 0.0,
785 retardation_factor: float = 1.0,
786 background: float = 0.0,
787 regional_flux: float = 0.0,
788 n_modes: int | None = None,
789 regularization_strength: float = 1e-10,
790 n_quad: int = 240,
791) -> npt.NDArray[np.floating]:
792 """Inverse of :func:`gamma_infiltration_to_extraction` (gamma-distributed screen velocity).
794 Returns
795 -------
796 ndarray, shape (n,)
797 Recovered injected concentration; NaN on extraction / rest bins.
798 """
799 pore_heights, weights = _velocity_gamma_streamtubes(screen_height, velocity_cv, n_bins)
800 return extraction_to_infiltration(
801 cout=cout,
802 flow=flow,
803 tedges=tedges,
804 cout_tedges=cout_tedges,
805 pore_heights=pore_heights,
806 porosity=porosity,
807 well_radius=well_radius,
808 longitudinal_dispersivity=longitudinal_dispersivity,
809 molecular_diffusivity=molecular_diffusivity,
810 retardation_factor=retardation_factor,
811 weights=weights,
812 background=background,
813 regional_flux=regional_flux,
814 n_modes=n_modes,
815 regularization_strength=regularization_strength,
816 n_quad=n_quad,
817 )
820def _velocity_gamma_streamtubes(
821 screen_height: float, velocity_cv: float, n_bins: int
822) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
823 """Streamtube pore heights and weights for a gamma-distributed screen velocity (mean velocity <-> H).
825 The layer velocity ratio ``rho`` is gamma(mean 1, std ``velocity_cv``); the effective pore height is
826 ``screen_height / rho`` (velocity ~ 1/height), so the mean velocity corresponds to height ``H``.
828 Returns
829 -------
830 pore_heights : ndarray
831 Effective streamtube pore heights ``screen_height / rho`` per velocity bin.
832 weights : ndarray
833 Probability mass per velocity bin.
835 Raises
836 ------
837 ValueError
838 If ``screen_height`` is not positive or ``velocity_cv`` is negative.
839 """
840 if screen_height <= 0.0:
841 msg = "screen_height must be positive"
842 raise ValueError(msg)
843 if velocity_cv < 0.0:
844 msg = "velocity_cv must be non-negative"
845 raise ValueError(msg)
846 if velocity_cv == 0.0:
847 # A degenerate gamma (std 0) is not a valid distribution; velocity_cv = 0 is the homogeneous
848 # screen -- a single streamtube at the mean velocity (pore height H), matching the doc.
849 return np.array([screen_height]), np.array([1.0])
850 bins = gamma.bins(mean=1.0, std=velocity_cv, n_bins=n_bins)
851 return screen_height / bins["expected_values"], bins["probability_mass"]