Coverage for src/gwtransport/radial_asr.py: 84%
126 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +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:
11volume coordinate ``V(r) = pi b n (r^2 - r_w^2)``, Scheidegger velocity-dependent dispersion
12``D = alpha_L |u| + D_m`` (microdispersion ``alpha_L |u|`` plus molecular diffusion ``D_m``),
13Kreft-Zuber flux boundary conditions, and the exact per-phase kernels (Airy for ``D_m = 0``; the
14log-derivative Riccati ODE for ``D_m > 0``). Nothing 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. Every signed-flow schedule (single cycle, multi-cycle ASR, intervening rests) is
19composed by the propagator-matrix engine (``gwtransport._radial_asr_reuse``), which carries the
20resident field across each flow reversal with the exact per-phase interior Green's functions
21(Airy / Riccati / Bessel); cycles are expressed through the flow sign pattern, not an argument.
22Molecular diffusion during pumping (the ``D_m > 0`` Whittaker kernel) is evaluated through the
23log-derivative Riccati ODE -- exact to the de Hoog inversion floor at any ``A_0/D_m``, with no
24special-function precision cap, and reducing continuously to the Airy branch as ``D_m -> 0``. During
25a **rest** (``Q = 0``) advection and microdispersion vanish and molecular diffusion acts alone on the
26wall-clock clock; it is carried exactly by the order-0 modified Bessel pure-diffusion kernel, the
27dominant mixing for seasonal storage / ATES. The only numerical steps are Gauss-Legendre quadrature
28and de Hoog Laplace inversion of exact special-function kernels. An independent finite-volume solve
29of the same PDE (``tests/src/_radial_asr_fv_oracle.py``) is used only as a test oracle.
31The reported ``cout`` is the flow-weighted average over each output bin -- defined on extraction bins
32(``flow < 0``) and ``NaN`` on injection / rest bins (nothing is recovered there).
34Macrodispersion within the well screen
35--------------------------------------
36The well screen has a **known** height; macrodispersion is the spread of arrival times caused by
37*velocity heterogeneity across the screen*. It is modelled as parallel streamtubes (``pore_heights``):
38each streamtube is an independent radial cell carrying the full flow, with an effective pore height
39that sets its velocity, and the output is the weight-averaged breakthrough. A streamtube of effective
40height ``b`` has velocity ``proportional to 1/b`` (its pore volume to radius ``r`` is
41``pi b n (r^2 - r_w^2)``), so smaller ``b`` means faster breakthrough.
42:func:`gamma_infiltration_to_extraction` builds this ensemble from a gamma distribution of the layer
43**velocity** within the fixed screen height (see that function); the mean velocity is set by the screen
44height and the spread by a velocity coefficient of variation. The spread is a within-screen velocity
45distribution -- velocity heterogeneity across the well screen -- not an aquifer pore-volume distribution.
47Regional background flow (drift)
48--------------------------------
49With a steady uniform regional Darcy flux ``regional_flux`` (``U``, drift seepage ``v_d = U/n``) the well
50field is superimposed on a regional gradient, so the stored bubble drifts and recovery degrades. The
51radial symmetry is broken and the transport is solved by an **azimuthal Fourier-mode** expansion
52``c(r, theta) = sum_m c_m(r) e^{i m theta}`` (``m = 0`` is the radial engine; drift couples ``m`` to
53``m +- 1``), composed through the same per-phase interior Green's functions
54(``gwtransport._radial_asr_drift_kernels``). ``regional_flux = 0`` (default) dispatches to the radial path
55bit-for-bit. The engine is for the **slow-drift** envelope -- the plume (including its rest-phase drift
56displacement) must stay well inside the stagnation radius ``r_s = |A_0|/|v_d|`` (else a ``ValueError``).
57Rest phases (``flow == 0``) are propagated by the exact free-space drift kernel (translate + anisotropic
58spread). The drift-induced recovery loss is validated against an independent 2-D finite-volume oracle.
60Available functions:
62- :func:`infiltration_to_extraction` - Forward transport: compute the extracted flux concentration ``cout`` from an
63 injected concentration ``cin``, a signed flow schedule (``> 0`` injection, ``< 0`` extraction, ``0`` rest) and the
64 well geometry. Grid-free -- the exact per-phase kernels are composed across every flow reversal, and the phases
65 (and hence the number of push-pull / ASR cycles) follow from the flow sign pattern rather than from an argument.
66 ``cout`` is the flow-weighted average over each output bin, in the units of ``cin``, defined on the extraction bins
67 and ``NaN`` on injection and rest bins; ``cout_tedges`` must equal ``tedges``. ``pore_heights`` is a scalar
68 (one homogeneous screen) or an array of streamtube heights -- each streamtube carries the full flow, and the
69 breakthroughs are averaged with ``weights`` (equal by default). Sorption is linear through ``retardation_factor``;
70 ``background`` is the ambient aquifer concentration, so the deviation ``cin - background`` is transported and
71 ``background`` added back. Nonzero ``regional_flux`` engages the azimuthal-mode drift engine within its
72 slow-drift envelope, and raises ``ValueError`` when the plume leaves that envelope.
74- :func:`extraction_to_infiltration` - Reverse operation: recover the injected concentration from measured extracted
75 concentrations under the same flow schedule and geometry. The forward operator is assembled column-by-column from
76 unit injection-pulse responses and inverted by Tikhonov least squares (``regularization_strength``). Returns one
77 value per bin: the recovered ``cin`` on injection bins and ``NaN`` on extraction and rest bins. ``NaN`` in ``cout``
78 on an extraction bin raises ``ValueError``; structural ``NaN`` on injection and rest bins is ignored.
80- :func:`gamma_infiltration_to_extraction` - Forward transport as in :func:`infiltration_to_extraction`, with the
81 streamtube ensemble built from a gamma distribution of the layer velocity across a well screen of known height
82 ``screen_height``. The velocity ratio has mean 1 and coefficient of variation ``velocity_cv``, a streamtube of
83 velocity ratio ``rho`` gets effective pore height ``screen_height / rho``, and the distribution is discretized into
84 ``n_bins`` equal-probability bins averaged by probability mass. ``velocity_cv = 0`` is the homogeneous screen: a
85 single streamtube at pore height ``screen_height``.
87- :func:`gamma_extraction_to_infiltration` - Reverse operation as in :func:`extraction_to_infiltration`, over the same
88 gamma layer-velocity streamtube ensemble as :func:`gamma_infiltration_to_extraction`.
90References
91----------
92The references below give the published closed-form solutions for the **single-phase** radial *injection*
93problem (steady divergent flow from one well) -- the per-phase forward kernel this module composes. The
94convergent-extraction dual and the multi-cycle push-pull / ASR composition across flow reversals are
95built on top of those kernels here and are not in the single-injection references. All
96share the assumptions used here: a single fully-penetrating well in a homogeneous medium with steady
97divergent flow ``v = Q / (2 pi b n r)``, plus retardation.
99The ``D_m = 0`` kernel (velocity-proportional microdispersion ``D = alpha_L |u|``, Airy functions)
100is the classical radial-dispersion problem: Tang & Babu (1979) under a Dirichlet (resident-concentration)
101well boundary, and Chen (1987) under the Cauchy / third-type (flux) boundary used here -- explicitly the
102Kreft-Zuber flux concentration, with transfer function ``Ai(Y) / [Ai(Y0)/2 - p^(1/3) Ai'(Y0)]`` equal to
103the flux operator this module evaluates. The ``D_m > 0`` kernel (``D = alpha_L |u| + D_m``, Kummer /
104confluent-hypergeometric functions) under the same flux boundary, with retardation, is Aichi & Akitaya
105(2018) -- whose well operator ``U(a,b) + 2a U(a+1,b+1)`` is this module's Whittaker flux boundary; they
106record the ``D_m -> 0`` reduction to Chen (1987) as an open problem, which this module performs
107continuously -- the log-derivative Riccati kernel reduces smoothly to the Airy branch as ``D_m -> 0``.
108The ``alpha_L = 0`` limit (constant diffusion, drift-dominated radial transport, Whittaker equation) is
109Akanji & Falade (2019). Each is an injection-only solution; none treats extraction or multi-cycle push-pull.
111Kreft, A., & Zuber, A. (1978). On the physical meaning of the dispersion equation and its solutions
112for different initial and boundary conditions. Chemical Engineering Science, 33(11), 1471-1480.
114Tang, D. H., & Babu, D. K. (1979). Analytical solution of a velocity dependent dispersion problem.
115Water Resources Research, 15(6), 1471-1478.
117Chen, C.-S. (1987). Analytical solutions for radial dispersion with Cauchy boundary at injection well.
118Water Resources Research, 23(7), 1217-1224.
120Aichi, M., & Akitaya, K. (2018). Analytical solution for a radial advection-dispersion equation
121including both mechanical dispersion and molecular diffusion for a steady-state flow field in a
122horizontal aquifer caused by a constant rate injection from a well. Hydrological Research Letters,
12312(3), 23-27.
125Akanji, L. T., & Falade, G. K. (2019). Closed-form solution of radial transport of tracers in porous
126media influenced by linear drift. Energies, 12(1), 29.
128This file is part of gwtransport which is released under AGPL-3.0 license.
129See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
130"""
132import numpy as np
133import numpy.typing as npt
134import pandas as pd
136from gwtransport import gamma
137from gwtransport._radial_asr_drift_kernels import _RS_FRAC, block_cout_deviation
138from gwtransport._radial_asr_reuse import cout_deviation
139from gwtransport._time import dt_to_days
140from gwtransport._validation import _validate_retardation_factor
143def _validate(
144 *,
145 cin_or_cout: npt.NDArray[np.floating],
146 flow: npt.NDArray[np.floating],
147 tedges: pd.DatetimeIndex,
148 cout_tedges: pd.DatetimeIndex,
149 pore_heights: npt.NDArray[np.floating],
150 porosity: float,
151 well_radius: float,
152 longitudinal_dispersivity: float,
153 molecular_diffusivity: float,
154 retardation_factor: float,
155 weights: npt.NDArray[np.floating] | None,
156 regional_flux: float,
157 n_modes: int | None,
158) -> None:
159 """Validate inputs for the radial single-well transport functions (signed flow is allowed).
161 Raises
162 ------
163 ValueError
164 On inconsistent lengths, non-positive geometry, out-of-range porosity, negative dispersion,
165 ``retardation_factor < 1``, mismatched ``weights``, NaN in ``flow``, a non-finite ``regional_flux``,
166 an ``n_modes < 1``, or a ``cout_tedges`` that differs from ``tedges`` (a distinct output grid is not
167 yet supported).
168 """
169 if len(tedges) != len(flow) + 1:
170 msg = "tedges must have one more element than flow"
171 raise ValueError(msg)
172 if len(cin_or_cout) != len(flow):
173 msg = "cin/cout must have the same length as flow"
174 raise ValueError(msg)
175 if not tedges.equals(cout_tedges):
176 msg = "cout_tedges must equal tedges (a distinct output grid is not yet supported)"
177 raise ValueError(msg)
178 if np.any(np.isnan(flow)):
179 msg = "flow contains NaN values, which are not allowed"
180 raise ValueError(msg)
181 if np.any(pore_heights <= 0.0):
182 msg = "pore_heights must be positive"
183 raise ValueError(msg)
184 if not 0.0 < porosity <= 1.0:
185 msg = "porosity must be in (0, 1]"
186 raise ValueError(msg)
187 if well_radius <= 0.0:
188 msg = "well_radius must be positive"
189 raise ValueError(msg)
190 if longitudinal_dispersivity <= 0.0:
191 msg = "longitudinal_dispersivity must be positive (the dispersion kernel requires alpha_L > 0)"
192 raise ValueError(msg)
193 if molecular_diffusivity < 0.0:
194 msg = "molecular_diffusivity must be non-negative"
195 raise ValueError(msg)
196 _validate_retardation_factor(retardation_factor)
197 if weights is not None and len(weights) != len(pore_heights):
198 msg = "weights must have the same length as pore_heights"
199 raise ValueError(msg)
200 if not np.isfinite(regional_flux):
201 msg = "regional_flux must be finite"
202 raise ValueError(msg)
203 if n_modes is not None and n_modes < 1:
204 msg = "n_modes must be >= 1"
205 raise ValueError(msg)
208def _auto_n_modes(
209 flow: npt.NDArray[np.floating],
210 dt_days: npt.NDArray[np.floating],
211 c_geo: float,
212 well_radius: float,
213 longitudinal_dispersivity: float,
214 v_d: float,
215 retardation_factor: float,
216) -> int:
217 r"""Azimuthal truncation ``M`` sized from the drift ratio ``eps`` and the rest-phase displacement.
219 The pumping-phase mode amplitudes decay geometrically, ``|c_m| ~ eps^|m|`` with
220 ``eps = v_d (R_b + delta) / A_0`` (worst-case ratio at the rest-displaced plume edge, so the
221 rest displacement ``delta`` feeds this bound too), so keeping modes ``-M .. M`` truncates the azimuthal field at
222 ``O(eps^{M+1})``; ``M`` is chosen so that tail is below ``~5e-3``. An interior rest phase
223 additionally translates the plume by ``delta = v_d t_rest / R`` (the free-space rest kernel's
224 own translation; idle bins before the first or after the last pumping do not move the field),
225 populating harmonics up to ``~ delta / width`` (``width`` the radial breakthrough std) -- the
226 second bound. The result is clamped to ``[2, 8]`` (the slow-drift envelope -- beyond
227 ``eps ~ 0.6`` the far-field escape this engine does not model dominates anyway); the rest
228 kernel's honest spectral-tail guard raises if a long rest still outruns the clamp (pass
229 ``n_modes`` explicitly then). ``A_0`` uses the **smallest** pumping magnitude (the worst-case
230 largest ``eps``, consistent with the stagnation-radius envelope guard), ``R_b`` the peak
231 **retarded** solute-front radius. Sizing ``R_b`` and ``delta`` from the retarded front (the
232 ``/R`` factors) makes the selected ``M`` invariant under the operator rescale
233 ``(A_0, v_d, D_m, R) -> (A_0/R, v_d/R, D_m/R, 1)`` -- the two equivalent runs pick the same
234 truncation instead of differing (see ``test_retardation_rescale_selects_same_n_modes``) -- and
235 is tighter (fewer wasted modes) than the water-front radius. This function is only reached for
236 nonzero drift.
238 Returns
239 -------
240 int
241 Azimuthal truncation ``M``.
242 """
243 pumping = np.abs(flow[flow != 0.0])
244 if pumping.size == 0: # all-rest schedule: nothing pumps; the engine returns all-NaN downstream
245 return 2
246 a0 = float(np.min(pumping)) / (2.0 * c_geo)
247 net_volume = np.concatenate(([0.0], np.cumsum(flow * dt_days)))
248 peak_volume = max(float(net_volume.max()), 0.0)
249 # Retarded solute-front radius: the solute penetrates R times less area than the water front,
250 # r_b^2 - r_w^2 = peak_volume / (R c_geo). This is the radius at which the azimuthal plume
251 # structure actually sits, and it makes r_b invariant under the (flow, v_d, R) rescale.
252 r_b = np.sqrt(well_radius**2 + peak_volume / (retardation_factor * c_geo))
253 nz = np.flatnonzero(flow != 0.0)
254 interior = slice(nz[0], nz[-1] + 1) # leading/trailing idle bins do not move the field
255 # Rest translation is v_d t / R (the free-space rest kernel's exact drift), not v_d t.
256 delta = abs(v_d) * float(np.sum(dt_days[interior][flow[interior] == 0.0])) / retardation_factor
257 eps = min(abs(v_d) * (r_b + delta) / abs(a0), _RS_FRAC)
258 # eps > 0 always holds here: v_d != 0 (this function is only reached for nonzero drift),
259 # r_b >= r_w > 0, and a0 > 0.
260 m_eps = int(np.ceil(np.log(5e-3) / np.log(eps)))
261 width = np.sqrt(longitudinal_dispersivity * r_b + longitudinal_dispersivity**2)
262 m_shift = int(np.ceil(delta / width)) + 2 if delta > 0.0 else 2
263 return int(np.clip(max(m_eps, m_shift), 2, 8))
266def _streamtube_ensemble(
267 cin_deviation: npt.NDArray[np.floating],
268 *,
269 flow: npt.NDArray[np.floating],
270 dt_days: npt.NDArray[np.floating],
271 c_geos: npt.NDArray[np.floating],
272 weights: npt.NDArray[np.floating],
273 porosity: float,
274 well_radius: float,
275 longitudinal_dispersivity: float,
276 molecular_diffusivity: float,
277 retardation_factor: float,
278 regional_flux: float,
279 n_modes: int | None,
280 n_quad: int,
281) -> npt.NDArray[np.floating]:
282 """Weight-averaged extracted-flux deviation over the streamtube ensemble.
284 Runs the transport engine once per streamtube (geometry constant ``c_geo = pi b n``) and averages by
285 ``weights``: the reused-propagator-matrix engine
286 (:func:`gwtransport._radial_asr_reuse.cout_deviation`) without regional drift, the azimuthal-mode
287 block engine (:func:`gwtransport._radial_asr_drift_kernels.block_cout_deviation`) with it. The drift
288 seepage ``v_d = U / n`` is the same for every streamtube (a regional Darcy flux through the porosity);
289 only the radial strength ``A_0 ~ 1/c_geo`` varies, so faster (thinner) streamtubes see a smaller drift
290 ratio, and ``n_modes`` is auto-sized per streamtube from that ratio when not given. ``cin_deviation``
291 may be ``(n,)`` or ``(n, k)`` -- a column batch is transported through one engine pass per streamtube
292 (the per-phase propagator / source / readout matrices are cin-independent, so they are built once and
293 applied to every column). Used by the forward (with ``cin``) and the reverse (with the unit-pulse
294 column batch).
296 Returns
297 -------
298 ndarray, shape (n,) or (n, k)
299 Weight-averaged extracted-flux deviation on extraction bins (matching ``cin_deviation``), ``0``
300 elsewhere.
301 """
302 v_d = regional_flux / porosity
303 # Both engines return NaN on injection / rest bins by contract (structural, expected). Only the
304 # extraction bins are accumulated, so those structural NaNs are dropped without a blanket nan_to_num --
305 # a genuine NaN on an EXTRACTION bin (a real numerical failure) then propagates and surfaces instead of
306 # silently reading as a physical zero.
307 ext_mask = flow < 0.0
308 acc = np.zeros(np.shape(cin_deviation))
309 for c_geo, w_i in zip(c_geos, weights, strict=True):
310 if regional_flux == 0.0:
311 dev = cout_deviation(
312 cin_deviation=cin_deviation,
313 flow=flow,
314 dt_days=dt_days,
315 c_geo=c_geo,
316 r_w=well_radius,
317 alpha_l=longitudinal_dispersivity,
318 molecular_diffusivity=molecular_diffusivity,
319 retardation_factor=retardation_factor,
320 n_quad=n_quad,
321 )
322 else:
323 # Steady regional drift breaks radial symmetry: the azimuthal-mode block engine carries it.
324 m = (
325 n_modes
326 if n_modes is not None
327 else _auto_n_modes(
328 flow, dt_days, c_geo, well_radius, longitudinal_dispersivity, v_d, retardation_factor
329 )
330 )
331 dev = block_cout_deviation(
332 cin_deviation=cin_deviation,
333 flow=flow,
334 dt_days=dt_days,
335 c_geo=c_geo,
336 r_w=well_radius,
337 alpha_l=longitudinal_dispersivity,
338 molecular_diffusivity=molecular_diffusivity,
339 retardation_factor=retardation_factor,
340 v_d=v_d,
341 n_modes=m,
342 n_quad=n_quad,
343 )
344 acc[ext_mask] += w_i * dev[ext_mask]
345 return acc / np.sum(weights)
348def infiltration_to_extraction(
349 *,
350 cin: npt.ArrayLike,
351 flow: npt.ArrayLike,
352 tedges: pd.DatetimeIndex,
353 cout_tedges: pd.DatetimeIndex,
354 pore_heights: npt.ArrayLike,
355 porosity: float,
356 well_radius: float,
357 longitudinal_dispersivity: float,
358 molecular_diffusivity: float = 0.0,
359 retardation_factor: float = 1.0,
360 weights: npt.ArrayLike | None = None,
361 background: float = 0.0,
362 regional_flux: float = 0.0,
363 n_modes: int | None = None,
364 n_quad: int = 240,
365) -> npt.NDArray[np.floating]:
366 """Compute the extracted flux concentration at a radial well for a signed flow schedule.
368 Parameters
369 ----------
370 cin : array-like, shape (n,)
371 Injected concentration per time bin (used only on injection bins, ``flow > 0``).
372 flow : array-like, shape (n,)
373 Signed flow per time bin [m^3/day]: ``> 0`` injection, ``< 0`` extraction, ``0`` rest.
374 tedges : DatetimeIndex
375 Time bin edges (``n + 1`` for ``n`` bins).
376 cout_tedges : DatetimeIndex
377 Output time bin edges; must equal ``tedges``. Output is NaN on injection / rest bins.
378 pore_heights : array-like
379 Effective streamtube pore height(s) ``b`` [m] -- a scalar (one homogeneous screen) or an array
380 of streamtube heights for the velocity-heterogeneity macrodispersion ensemble (each streamtube
381 carries the full flow; smaller ``b`` = faster). See the module docstring and
382 :func:`gamma_infiltration_to_extraction`.
383 porosity : float
384 Porosity ``n`` [-].
385 well_radius : float
386 Well (screen) radius ``r_w`` [m].
387 longitudinal_dispersivity : float
388 Longitudinal dispersivity ``alpha_L`` [m].
389 molecular_diffusivity : float, optional
390 Molecular diffusivity ``D_m`` [m^2/day]. Default 0. ``D_m = 0`` uses the vectorized Airy branch;
391 ``D_m > 0`` uses the log-derivative Riccati kernel -- exact to the de Hoog floor at any ``A_0/D_m``
392 with no precision cap, reducing continuously to the Airy branch as ``D_m -> 0``. For ``D_m > 0`` each
393 one-signed pumping phase is advanced at its (constant) dt-weighted mean flow magnitude, so
394 within-phase variable flow is a constant-|Q| approximation whose error grows with the within-phase
395 flow variation; it is exact for constant-|Q| phases, and (at any flow schedule) for ``D_m = 0``.
396 retardation_factor : float, optional
397 Linear retardation ``R >= 1``. Default 1.
398 weights : array-like, optional
399 Per-streamtube averaging weights (same length as ``pore_heights``). Default equal weights.
400 background : float, optional
401 Ambient aquifer concentration ``c_bg``. The deviation ``cin - c_bg`` is transported and
402 ``c_bg`` is added back; constant ``cin = c_bg`` returns ``cout = c_bg``. Default 0.
403 regional_flux : float, optional
404 Steady uniform regional background Darcy flux ``U`` [m/day] in ``+x`` (drift seepage
405 ``v_d = U / n``). ``0`` (default) reproduces the radial-symmetric engine bit-for-bit. A nonzero
406 value engages the azimuthal-mode block engine, which captures the drift-induced recovery loss
407 (the down-gradient plume is partly swept past the well). The slow-drift envelope requires the
408 plume -- including its rest-phase drift displacement -- to stay well inside the stagnation radius
409 ``r_s = |A_0| / |v_d|`` (a ``ValueError`` is raised otherwise). Rest phases (``flow == 0``) are
410 propagated by the exact free-space drift kernel (translation ``v_d t / R`` plus anisotropic
411 Gaussian spread, with a Neumann-image closure at the shut well face). See
412 :ref:`concept-drift-envelope` for a worked multi-year feasibility table.
413 n_modes : int, optional
414 Azimuthal truncation ``M`` for the drift engine (keeps modes ``-M .. M``). Default ``None``
415 auto-sizes ``M`` from the plume-front drift ratio ``eps = v_d R_b / A_0`` and the rest-phase
416 displacement (clamped to ``[2, 8]``). Ignored when ``regional_flux == 0``.
417 n_quad : int, optional
418 Gauss-Legendre node count for the resident-profile superposition. Default 240.
420 Returns
421 -------
422 ndarray, shape (n,)
423 Extracted flux concentration; NaN on injection and rest bins.
424 """
425 cin = np.asarray(cin, dtype=float)
426 flow = np.asarray(flow, dtype=float)
427 pore_heights = np.atleast_1d(np.asarray(pore_heights, dtype=float))
428 weights_arr = np.ones(len(pore_heights)) if weights is None else np.atleast_1d(np.asarray(weights, dtype=float))
429 _validate(
430 cin_or_cout=cin,
431 flow=flow,
432 tedges=tedges,
433 cout_tedges=cout_tedges,
434 pore_heights=pore_heights,
435 porosity=porosity,
436 well_radius=well_radius,
437 longitudinal_dispersivity=longitudinal_dispersivity,
438 molecular_diffusivity=molecular_diffusivity,
439 retardation_factor=retardation_factor,
440 weights=None if weights is None else weights_arr,
441 regional_flux=regional_flux,
442 n_modes=n_modes,
443 )
444 ext_mask = flow < 0.0
445 cout_dev = _streamtube_ensemble(
446 cin - background,
447 flow=flow,
448 dt_days=dt_to_days(tedges),
449 c_geos=np.pi * pore_heights * porosity,
450 weights=weights_arr,
451 porosity=porosity,
452 well_radius=well_radius,
453 longitudinal_dispersivity=longitudinal_dispersivity,
454 molecular_diffusivity=molecular_diffusivity,
455 retardation_factor=retardation_factor,
456 regional_flux=regional_flux,
457 n_modes=n_modes,
458 n_quad=n_quad,
459 )
460 cout = np.full(len(flow), np.nan)
461 cout[ext_mask] = background + cout_dev[ext_mask]
462 return cout
465def extraction_to_infiltration(
466 *,
467 cout: npt.ArrayLike,
468 flow: npt.ArrayLike,
469 tedges: pd.DatetimeIndex,
470 cout_tedges: pd.DatetimeIndex,
471 pore_heights: npt.ArrayLike,
472 porosity: float,
473 well_radius: float,
474 longitudinal_dispersivity: float,
475 molecular_diffusivity: float = 0.0,
476 retardation_factor: float = 1.0,
477 weights: npt.ArrayLike | None = None,
478 background: float = 0.0,
479 regional_flux: float = 0.0,
480 n_modes: int | None = None,
481 regularization_strength: float = 1e-10,
482 n_quad: int = 240,
483) -> npt.NDArray[np.floating]:
484 """Recover the injected concentration from extracted-water measurements (Tikhonov inverse).
486 Inverts the forward operator built by :func:`infiltration_to_extraction`. Returns the injected
487 concentration on injection bins (NaN on extraction / rest bins).
489 Parameters
490 ----------
491 cout : array-like, shape (n,)
492 Measured extracted concentration (used on extraction bins, ``flow < 0``).
493 flow, tedges, cout_tedges, pore_heights, porosity, well_radius, longitudinal_dispersivity
494 As in :func:`infiltration_to_extraction`.
495 molecular_diffusivity, retardation_factor, weights, background, regional_flux, n_modes, n_quad
496 As in :func:`infiltration_to_extraction`.
497 regularization_strength : float, optional
498 Tikhonov parameter, must be non-negative. Default ``1e-10``.
500 Returns
501 -------
502 ndarray, shape (n,)
503 Recovered injected concentration; NaN on extraction / rest bins.
505 Raises
506 ------
507 ValueError
508 If ``cout`` contains NaN on any extraction bin (``flow < 0``), which would poison the
509 least-squares solve. Structural NaN on injection / rest bins is allowed.
510 If ``regularization_strength`` is negative: it would reach ``np.sqrt`` in the augmented
511 system and silently produce an all-NaN ``cin``.
512 """
513 cout = np.asarray(cout, dtype=float)
514 flow = np.asarray(flow, dtype=float)
515 pore_heights = np.atleast_1d(np.asarray(pore_heights, dtype=float))
516 weights_arr = np.ones(len(pore_heights)) if weights is None else np.atleast_1d(np.asarray(weights, dtype=float))
517 _validate(
518 cin_or_cout=cout,
519 flow=flow,
520 tedges=tedges,
521 cout_tedges=cout_tedges,
522 pore_heights=pore_heights,
523 porosity=porosity,
524 well_radius=well_radius,
525 longitudinal_dispersivity=longitudinal_dispersivity,
526 molecular_diffusivity=molecular_diffusivity,
527 retardation_factor=retardation_factor,
528 weights=None if weights is None else weights_arr,
529 regional_flux=regional_flux,
530 n_modes=n_modes,
531 )
532 # A NaN measurement on an extraction bin (the only bins the inverse reads) would poison the whole
533 # least-squares solve into an all-NaN cin; raise instead, as the advection / diffusion inverses do.
534 # Structural NaN on injection / rest bins is allowed (those bins are ignored by the inverse).
535 if np.any(np.isnan(cout[flow < 0.0])):
536 msg = "cout contains NaN values on extraction bins, which are not allowed"
537 raise ValueError(msg)
538 if regularization_strength < 0.0:
539 msg = f"regularization_strength must be >= 0, got {regularization_strength}"
540 raise ValueError(msg)
541 # Build the dense forward operator W whose columns are the unit-injection-pulse responses. The per-phase
542 # propagator / source / readout matrices are cin-independent (flow + geometry only), so the whole
543 # unit-pulse batch is transported in ONE engine pass -- built once and applied to every column.
544 inj_mask, ext_mask = flow > 0.0, flow < 0.0
545 inj_idx = np.flatnonzero(inj_mask)
546 pulses = np.zeros((len(flow), len(inj_idx)))
547 pulses[inj_idx, np.arange(len(inj_idx))] = 1.0
548 cols = _streamtube_ensemble(
549 pulses,
550 flow=flow,
551 dt_days=dt_to_days(tedges),
552 c_geos=np.pi * pore_heights * porosity,
553 weights=weights_arr,
554 porosity=porosity,
555 well_radius=well_radius,
556 longitudinal_dispersivity=longitudinal_dispersivity,
557 molecular_diffusivity=molecular_diffusivity,
558 retardation_factor=retardation_factor,
559 regional_flux=regional_flux,
560 n_modes=n_modes,
561 n_quad=n_quad,
562 )
563 w_ens = cols[ext_mask, :]
564 # Tikhonov least-squares min ||W x - (cout-bg)||^2 + lambda ||x||^2 via the stable augmented
565 # system [W; sqrt(lambda) I] x = [cout-bg; 0]. The forward operator has column sums ~1
566 # (mass conservation per injection bin) and overdetermined rows, so a direct Tikhonov fit is used.
567 n_inj = w_ens.shape[1]
568 augmented = np.vstack([w_ens, np.sqrt(regularization_strength) * np.eye(n_inj)])
569 rhs = np.concatenate([cout[ext_mask] - background, np.zeros(n_inj)])
570 cin_dev = np.linalg.lstsq(augmented, rhs, rcond=None)[0]
571 cin = np.full(len(flow), np.nan)
572 cin[inj_mask] = background + cin_dev
573 return cin
576def gamma_infiltration_to_extraction(
577 *,
578 cin: npt.ArrayLike,
579 flow: npt.ArrayLike,
580 tedges: pd.DatetimeIndex,
581 cout_tedges: pd.DatetimeIndex,
582 porosity: float,
583 well_radius: float,
584 longitudinal_dispersivity: float,
585 screen_height: float,
586 velocity_cv: float,
587 n_bins: int = 100,
588 molecular_diffusivity: float = 0.0,
589 retardation_factor: float = 1.0,
590 background: float = 0.0,
591 regional_flux: float = 0.0,
592 n_modes: int | None = None,
593 n_quad: int = 240,
594) -> npt.NDArray[np.floating]:
595 """Radial transport with gamma-distributed screen velocity (within-screen macrodispersion).
597 The well screen has a **known** height ``screen_height``; macrodispersion is the spread of arrival
598 times from velocity heterogeneity across that fixed height. The layer velocity is gamma-distributed
599 with mean equal to the homogeneous value (a streamtube at the mean velocity has effective pore
600 height ``screen_height``) and coefficient of variation ``velocity_cv``. A streamtube with velocity
601 ratio ``rho`` (gamma, mean 1) has effective pore height ``screen_height / rho`` -- faster layers are
602 thinner and break through sooner. The gamma is discretized into ``n_bins`` equal-probability bins
603 (:func:`gwtransport.gamma.bins`) and averaged by probability mass via
604 :func:`infiltration_to_extraction`.
606 Parameters
607 ----------
608 screen_height : float
609 Known well-screen height ``H`` [m] (the fixed total; the mean streamtube velocity is set by it).
610 velocity_cv : float
611 Coefficient of variation of the layer velocity (the macrodispersion strength). ``0`` is a
612 homogeneous screen (a single streamtube, sharp breakthrough); typically ``< 1`` -- larger values
613 give a heavy slow-velocity tail.
614 n_bins : int, optional
615 Number of equal-probability velocity bins. Default 100.
616 cin, flow, tedges, cout_tedges, porosity, well_radius, longitudinal_dispersivity
617 As in :func:`infiltration_to_extraction`.
618 molecular_diffusivity, retardation_factor, background, regional_flux, n_modes, n_quad
619 As in :func:`infiltration_to_extraction`.
621 Returns
622 -------
623 ndarray, shape (n,)
624 Extracted flux concentration; NaN on injection / rest bins.
625 """
626 pore_heights, weights = _velocity_gamma_streamtubes(screen_height, velocity_cv, n_bins)
627 return infiltration_to_extraction(
628 cin=cin,
629 flow=flow,
630 tedges=tedges,
631 cout_tedges=cout_tedges,
632 pore_heights=pore_heights,
633 porosity=porosity,
634 well_radius=well_radius,
635 longitudinal_dispersivity=longitudinal_dispersivity,
636 molecular_diffusivity=molecular_diffusivity,
637 retardation_factor=retardation_factor,
638 weights=weights,
639 background=background,
640 regional_flux=regional_flux,
641 n_modes=n_modes,
642 n_quad=n_quad,
643 )
646def gamma_extraction_to_infiltration(
647 *,
648 cout: npt.ArrayLike,
649 flow: npt.ArrayLike,
650 tedges: pd.DatetimeIndex,
651 cout_tedges: pd.DatetimeIndex,
652 porosity: float,
653 well_radius: float,
654 longitudinal_dispersivity: float,
655 screen_height: float,
656 velocity_cv: float,
657 n_bins: int = 100,
658 molecular_diffusivity: float = 0.0,
659 retardation_factor: float = 1.0,
660 background: float = 0.0,
661 regional_flux: float = 0.0,
662 n_modes: int | None = None,
663 regularization_strength: float = 1e-10,
664 n_quad: int = 240,
665) -> npt.NDArray[np.floating]:
666 """Inverse of :func:`gamma_infiltration_to_extraction` (gamma-distributed screen velocity).
668 Returns
669 -------
670 ndarray, shape (n,)
671 Recovered injected concentration; NaN on extraction / rest bins.
672 """
673 pore_heights, weights = _velocity_gamma_streamtubes(screen_height, velocity_cv, n_bins)
674 return extraction_to_infiltration(
675 cout=cout,
676 flow=flow,
677 tedges=tedges,
678 cout_tedges=cout_tedges,
679 pore_heights=pore_heights,
680 porosity=porosity,
681 well_radius=well_radius,
682 longitudinal_dispersivity=longitudinal_dispersivity,
683 molecular_diffusivity=molecular_diffusivity,
684 retardation_factor=retardation_factor,
685 weights=weights,
686 background=background,
687 regional_flux=regional_flux,
688 n_modes=n_modes,
689 regularization_strength=regularization_strength,
690 n_quad=n_quad,
691 )
694def _velocity_gamma_streamtubes(
695 screen_height: float, velocity_cv: float, n_bins: int
696) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
697 """Streamtube pore heights and weights for a gamma-distributed screen velocity (mean velocity <-> H).
699 The layer velocity ratio ``rho`` is gamma(mean 1, std ``velocity_cv``); the effective pore height is
700 ``screen_height / rho`` (velocity ~ 1/height), so the mean velocity corresponds to height ``H``.
702 Returns
703 -------
704 pore_heights : ndarray
705 Effective streamtube pore heights ``screen_height / rho`` per velocity bin.
706 weights : ndarray
707 Probability mass per velocity bin.
709 Raises
710 ------
711 ValueError
712 If ``screen_height`` is not positive or ``velocity_cv`` is negative.
713 """
714 if screen_height <= 0.0:
715 msg = "screen_height must be positive"
716 raise ValueError(msg)
717 if velocity_cv < 0.0:
718 msg = "velocity_cv must be non-negative"
719 raise ValueError(msg)
720 if velocity_cv == 0.0:
721 # A degenerate gamma (std 0) is not a valid distribution; velocity_cv = 0 is the homogeneous
722 # screen -- a single streamtube at the mean velocity (pore height H), matching the doc.
723 return np.array([screen_height]), np.array([1.0])
724 bins = gamma.bins(mean=1.0, std=velocity_cv, n_bins=n_bins)
725 return screen_height / bins["expected_values"], bins["probability_mass"]