Coverage for src/gwtransport/recharge.py: 100%
173 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1"""
2Recharge-Driven Transport for Aquifers with Areal Recharge.
4Concentration at extraction has two sources. 1) Water infiltrates and is
5transported through an aquifer with constant thickness to extraction. 2) During
6transport, rainfall is mixed instantaneously over the height of the aquifer. In
7an unbounded aquifer all extracted water originates as recharge. Transport is
8advective with linear sorption; there is no microdispersion, molecular diffusion,
9or macrodispersion. Only forward modeling is supported. No assumption is made
10about whether the flow is radial or orthogonal. Two conceptual models share one
11entry point:
13- **Unbounded aquifer** (``aquifer_pore_volume=None``): all extracted water
14 originates as recharge. The residence-time distribution is exponential with
15 mean ``retardation_factor * aquifer_pore_depth / N`` — independent of the
16 pumping rate, hydraulic conductivity, capture-zone size, and planform shape
17 (Haitjema, 1995). In the cumulative-recharge clock
18 ``u(t) = ∫ N dt / (retardation_factor * aquifer_pore_depth)`` (pore volumes
19 flushed) the model is the stationary unit filter ``dC/du = cin_recharge - C``,
20 which this module integrates in closed form per bin. No flow rate is needed.
22- **Bounded aquifer** (``aquifer_pore_volume`` set): the aquifer extent is
23 capped at pore volume ``aquifer_pore_volume`` (strip area
24 ``aquifer_pore_volume / aquifer_pore_depth``). Water with concentration
25 ``cin`` enters at the upstream side at rate ``q_b = flow - N * area``
26 whenever extraction exceeds the rainfall on the strip. When rainfall exceeds
27 extraction (``q_b < 0``) the surplus flows out across the upstream boundary
28 and is lost; the outside has no memory, so when extraction later dominates
29 again the inflow carries the current ``cin``. The exact solution is the
30 unbounded exponential kernel acting on ``cin_recharge``, truncated at the
31 boundary-entry time of the extracted water, with the residual tail weight
32 placed as an atom on ``cin`` at the entry time. With zero recharge this
33 reduces exactly to single-pore-volume piston flow
34 (:func:`gwtransport.advection.infiltration_to_extraction`); with the
35 boundary never feeding the well it reduces exactly to the unbounded model.
37Available functions:
39- :func:`recharge_to_extraction` - Compute extracted concentration from
40 recharge concentration (and, in the bounded model, upstream-boundary
41 concentration). Exact closed-form solution; output is a flow-weighted
42 (bounded) or recharge-weighted (unbounded) bin average.
44References
45----------
46Haitjema, H.M. (1995). On the residence time distribution in idealized
47groundwatersheds. Journal of Hydrology, 172(1-4), 127-146.
48https://doi.org/10.1016/0022-1694(95)02732-5
50This file is part of gwtransport which is released under AGPL-3.0 license.
51See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
52"""
54import numpy as np
55import numpy.typing as npt
56import pandas as pd
58from gwtransport._time import tedges_to_days
59from gwtransport._validation import (
60 _validate_no_nan,
61 _validate_non_negative_array,
62 _validate_positive_scalar,
63 _validate_retardation_factor,
64 _validate_tedges_parity,
65)
67# Kernel weights older than this many pore-volume flushes are below one ulp of the
68# row sum (e^-60 ~ 9e-27); truncating them keeps the gathers banded without
69# changing any double-precision result.
70_KERNEL_CUTOFF = 60.0
73def recharge_to_extraction(
74 *,
75 cin: npt.ArrayLike | None = None,
76 cin_recharge: npt.ArrayLike,
77 flow: npt.ArrayLike | None = None,
78 recharge: npt.ArrayLike,
79 tedges: pd.DatetimeIndex,
80 cout_tedges: pd.DatetimeIndex,
81 aquifer_pore_volume: float | None = None,
82 aquifer_pore_depth: float,
83 retardation_factor: float = 1.0,
84) -> npt.NDArray[np.floating]:
85 """Compute the concentration of extracted water under uniform areal recharge.
87 Unbounded model (``aquifer_pore_volume=None``): exponential residence-time
88 distribution with mean ``retardation_factor * aquifer_pore_depth / N``
89 (Haitjema, 1995), exact for bin-constant inputs. Bounded model
90 (``aquifer_pore_volume`` set, together with ``cin`` and ``flow``): the
91 exponential kernel is truncated at the upstream-boundary entry time and the
92 residual weight is an atom on ``cin``; water pushed out across the boundary
93 during rainfall surplus is lost.
95 Parameters
96 ----------
97 cin : array-like, optional
98 Concentration of the water entering at the upstream side of the
99 bounded aquifer [concentration units]. Required when
100 ``aquifer_pore_volume`` is set; must be None otherwise.
101 cin_recharge : array-like
102 Concentration of the recharge water entering via the surface
103 [concentration units]. Length must equal ``len(tedges) - 1``; constant
104 over each interval ``[tedges[i], tedges[i+1])``.
105 flow : array-like, optional
106 Extraction rate [m3/day]. Required when ``aquifer_pore_volume`` is
107 set; must be None otherwise, because the unbounded model is
108 independent of the pumping rate (see Notes). Must be non-negative and
109 NaN-free.
110 recharge : array-like
111 Areal recharge rate N [m/day; same length unit as
112 ``aquifer_pore_depth``]. Length must equal ``len(tedges) - 1``.
113 Must be non-negative and NaN-free.
114 tedges : pandas.DatetimeIndex
115 Time bin edges for the input series.
116 cout_tedges : pandas.DatetimeIndex
117 Time bin edges for the output series. Bins not fully inside the
118 ``tedges`` range return NaN.
119 aquifer_pore_volume : float, optional
120 Pore volume of the bounded aquifer [m3]. The strip area between the
121 upstream boundary and the well is
122 ``aquifer_pore_volume / aquifer_pore_depth``. Default None (unbounded).
123 aquifer_pore_depth : float
124 Pore volume per unit surface area: porosity times saturated thickness
125 [m]. The only static aquifer parameter of the unbounded model.
126 retardation_factor : float, optional
127 Compound retardation factor (>= 1.0), by default 1.0. Dilates the
128 solute clock; mixing fractions are unaffected.
130 Returns
131 -------
132 numpy.ndarray
133 Extracted concentration per ``cout_tedges`` bin, length
134 ``len(cout_tedges) - 1``. Flow-weighted bin average (bounded model) or
135 recharge-weighted bin average (unbounded model). NaN for bins outside
136 the input time range, for zero-recharge bins (unbounded), and for
137 zero-extraction bins (bounded).
139 Raises
140 ------
141 ValueError
142 If array lengths do not match the bin-edge pattern, inputs contain NaN
143 or negative values, physical parameters are out of range, or only part
144 of the bounded-model triple (``cin``, ``flow``,
145 ``aquifer_pore_volume``) is provided.
147 See Also
148 --------
149 gwtransport.advection.infiltration_to_extraction : Zero-recharge limit of the bounded model.
150 gwtransport.deposition.deposition_to_extraction : Distributed source along the flow path.
151 :ref:`concept-residence-time` : Background on residence times.
152 :ref:`concept-transport-equation` : Flow-weighted averaging approach.
154 Notes
155 -----
156 The unbounded model needs no flow rate because the capture zone
157 self-adjusts: the well always draws exactly its pumping rate from
158 recharge, over a capture area ``flow / N``. Pumping harder widens the
159 capture area proportionally, leaving the age composition of the extracted
160 water -- set by the ratio of pore storage per unit area
161 (``aquifer_pore_depth``) to recharge per unit area (``N``) -- unchanged,
162 so the flow rate cancels exactly (Haitjema, 1995). In the bounded model
163 the area is fixed by ``aquifer_pore_volume`` instead of adjusting to the
164 well, so the flow rate no longer cancels and must be given.
166 Spin-up follows the ``"constant"`` policy: all inputs are treated as
167 constant at their first values before ``tedges[0]``. For the bounded model
168 this is the steady concentration profile
169 ``C(V) = cr0 + (cin0 - cr0) * (V_R - apv) / (V_R - V)`` when the boundary
170 feeds the well (``q_b(0) > 0``, ``V_R = flow[0] * aquifer_pore_depth /
171 recharge[0]``), and the uniform profile ``cin_recharge[0]`` otherwise.
173 Under constant inputs with ``flow > N * area`` the extracted water is the
174 mass-balance mixture ``cin_recharge + (cin - cin_recharge) * q_b / flow``:
175 an exponential residence-time density carrying the recharge fraction plus
176 a piston atom of mass ``q_b / flow`` at the boundary-to-well travel time.
178 The exponential kernel lives on the dimensionless clock ``u`` and is
179 parameter-free; the pumping rate enters the bounded model only through the
180 boundary-entry times. All formulas are closed-form (exp/log of bin-local
181 quantities), exact to machine precision for bin-constant inputs.
183 References
184 ----------
185 Haitjema, H.M. (1995). On the residence time distribution in idealized
186 groundwatersheds. Journal of Hydrology, 172(1-4), 127-146.
187 https://doi.org/10.1016/0022-1694(95)02732-5
189 Examples
190 --------
191 >>> import numpy as np
192 >>> import pandas as pd
193 >>> from gwtransport.recharge import recharge_to_extraction
194 >>> tedges = pd.date_range("2020-01-01", periods=11, freq="D")
195 >>> cout = recharge_to_extraction(
196 ... cin_recharge=np.full(10, 2.5),
197 ... recharge=np.full(10, 0.002),
198 ... tedges=tedges,
199 ... cout_tedges=tedges[3:],
200 ... aquifer_pore_depth=3.0,
201 ... )
202 >>> np.allclose(cout, 2.5)
203 True
204 """
205 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
206 cr = np.asarray(cin_recharge, dtype=float)
207 rech = np.asarray(recharge, dtype=float)
208 bounded = aquifer_pore_volume is not None
210 if (cin is None) != (flow is None) or (cin is None) == bounded:
211 msg = "cin, flow, and aquifer_pore_volume must be provided together (bounded model) or all be None"
212 raise ValueError(msg)
213 _validate_tedges_parity(tedges, cr, tedges_name="tedges", values_name="cin_recharge")
214 _validate_tedges_parity(tedges, rech, tedges_name="tedges", values_name="recharge")
215 _validate_no_nan(cr, name="cin_recharge")
216 _validate_no_nan(rech, name="recharge")
217 _validate_non_negative_array(rech, name="recharge")
218 _validate_positive_scalar(aquifer_pore_depth, name="aquifer_pore_depth")
219 _validate_retardation_factor(retardation_factor)
221 t = tedges_to_days(tedges)
222 tq = tedges_to_days(cout_tedges, ref=tedges[0])
223 dt = np.diff(t)
224 k = rech / (retardation_factor * aquifer_pore_depth)
225 u = np.concatenate([[0.0], np.cumsum(k * dt)])
226 covered = (tq[:-1] >= t[0]) & (tq[1:] <= t[-1])
228 if not bounded:
229 # The unbounded aquifer is the no-boundary special case: an infinite
230 # pore volume puts the boundary beyond reach (pure kernel and spin-up
231 # terms), and a synthetic flow proportional to recharge makes the
232 # flow-weighted bin average the recharge-weighted one.
233 return _bounded_average(t=t, dt=dt, u=u, k=k, q=k, cr=cr, cb=cr, apv=np.inf, tq=tq, covered=covered)
235 cb = np.asarray(cin, dtype=float)
236 q = np.asarray(flow, dtype=float)
237 _validate_tedges_parity(tedges, cb, tedges_name="tedges", values_name="cin")
238 _validate_tedges_parity(tedges, q, tedges_name="tedges", values_name="flow")
239 _validate_no_nan(cb, name="cin")
240 _validate_no_nan(q, name="flow")
241 _validate_non_negative_array(q, name="flow")
242 _validate_positive_scalar(aquifer_pore_volume, name="aquifer_pore_volume")
243 # The solute clock divides flow by the retardation factor; weighting ratios are unaffected.
244 return _bounded_average(
245 t=t,
246 dt=dt,
247 u=u,
248 k=k,
249 q=q / retardation_factor,
250 cr=cr,
251 cb=cb,
252 apv=float(aquifer_pore_volume),
253 tq=tq,
254 covered=covered,
255 )
258def _arrival_times(*, t, dt, k, q, apv):
259 """Forward arrival time at the well of parcels released at ``(t[j], V=apv)``.
261 Trajectories obey ``dV/dt = -(q - k V)`` with the per-bin closed form
262 ``V(s) = V_R + (V - V_R) e^{k s}``. Returns NaN for parcels that are
263 expelled across the boundary (lost) or have not arrived by ``t[-1]``.
264 These arrivals are exactly the output times where the boundary-entry bin
265 of the extracted water changes, including the pre-record transition. The
266 bin loop over the live front is O(n^2) in the number of input bins.
268 Returns
269 -------
270 ndarray
271 Arrival time in days per release edge ``t[0..n-1]``; NaN if the parcel
272 never reaches the well within the record.
273 """
274 n = len(dt)
275 arrivals = np.full(n, np.nan)
276 pos = np.full(n, np.nan)
277 with np.errstate(over="ignore", invalid="ignore"):
278 for i in range(n):
279 pos[i] = apv
280 idx = np.nonzero(np.isfinite(pos[: i + 1]))[0]
281 vi = pos[idx]
282 if k[i] > 0:
283 v_r = q[i] / k[i]
284 v_end = v_r + (vi - v_r) * np.exp(k[i] * dt[i])
285 hit = v_end <= 0.0
286 s_hit = np.log(v_r / (v_r - vi[hit])) / k[i]
287 else:
288 v_end = vi - q[i] * dt[i]
289 hit = v_end <= 0.0
290 s_hit = vi[hit] / q[i] if q[i] > 0 else vi[hit][:0] # q == 0 cannot hit the well
291 arrivals[idx[hit]] = t[i] + s_hit
292 # Resolved (arrived) parcels and parcels at or beyond the boundary
293 # (expelled, or parked exactly on a stagnant boundary -- the next
294 # edge's release duplicates a parked parcel) drop out of the front.
295 pos[idx] = np.where(hit | (v_end >= apv), np.nan, v_end)
296 return arrivals
299def _backward_entries(*, queries, t, dt, k, q, apv, v_start=0.0, edge_side="right"):
300 """Trace backward characteristics from ``(query, V=v_start)`` to the boundary or to ``t[0]``.
302 Returns ``(s, v0)``: the boundary-entry time ``s`` (NaN if the water
303 predates the record) and the landing position ``v0`` at ``t[0]`` for
304 pre-record water (NaN otherwise). Uses the numerically local per-bin form
305 ``V_a = V_R + (V_b - V_R) e^{-k seg}``; differencing global accumulators
306 instead loses precision catastrophically once ``u >> 1``.
308 ``v_start=apv`` with ``edge_side="left"`` traces the grazing trajectory
309 that touches the boundary exactly at an (on-edge) query time backward into
310 the preceding bin, yielding the left-branch limit of the entry-time map at
311 that arrival (its release time); for queries not preceded by outflow the
312 walk exits immediately and returns the query time itself.
314 Returns
315 -------
316 tuple of ndarray
317 ``(s, v0)`` per query: boundary-entry time in days (NaN for pre-record
318 water) and landing position at ``t[0]`` (NaN for entered water).
319 """
320 n = len(dt)
321 nq = len(queries)
322 m = np.clip(np.searchsorted(t, queries, side=edge_side) - 1, 0, n - 1)
323 s_out = np.full(nq, np.nan)
324 v0 = np.full(nq, np.nan)
325 pos = np.full(nq, float(v_start))
326 open_ = np.ones(nq, dtype=bool)
327 with np.errstate(invalid="ignore", divide="ignore"): # stagnant boundary (vi == v_r == apv)
328 for i in range(n - 1, -1, -1):
329 sel = np.nonzero(open_ & (m >= i))[0]
330 if sel.size == 0:
331 continue
332 starts_here = m[sel] == i
333 seg = np.where(starts_here, queries[sel] - t[i], dt[i])
334 t_hi = np.where(starts_here, queries[sel], t[i + 1])
335 vi = pos[sel]
336 if k[i] > 0:
337 v_r = q[i] / k[i]
338 va = v_r + (vi - v_r) * np.exp(-k[i] * seg)
339 ent = va >= apv
340 back = -np.log((apv - v_r) / (vi[ent] - v_r)) / k[i]
341 back = np.where(np.isfinite(back), back, 0.0)
342 else:
343 va = vi + q[i] * seg
344 # The q > 0 guard keeps a parcel parked exactly on the boundary
345 # through a fully stagnant bin (no flow, no recharge) walking into
346 # earlier bins: nothing moves and nothing is lost there.
347 ent = (va >= apv) & (q[i] > 0.0)
348 back = (apv - vi[ent]) / q[i] if q[i] > 0 else vi[ent][:0]
349 s_ent = np.clip(t_hi[ent] - back, t[i], t_hi[ent])
350 s_out[sel[ent]] = s_ent
351 open_[sel[ent]] = False
352 pos[sel[~ent]] = va[~ent]
353 v0[open_] = pos[open_]
354 return s_out, v0
357def _bounded_average(*, t, dt, u, k, q, cr, cb, apv, tq, covered):
358 """Flow-weighted bin averages of the bounded model, exact via piece integration.
360 The output interval is split at: input edges, output edges, and the
361 arrival times of boundary parcels released at the input edges. Between
362 consecutive breakpoints the entry time stays within one source bin and
363 every term integrates in closed form. The boundary atom integrates to
364 ``cin_js * q_b,js * (s2 - s1)`` via the change of variables
365 ``q e^{-u(t)} dt = e^{-u(s)} q_b(s) ds`` along the entry map; the kernel
366 terms reduce to bin-local exponentials times
367 ``I2 = ∫ q e^{-(u(t)-u(t2))} dt``; pre-record water carries the steady
368 spin-up profile evaluated at the landing position ``v0 = G(t)``, whose
369 flow-weighted integral is ``cr0``-linear plus a closed-form logarithm.
371 Returns
372 -------
373 ndarray
374 Flow-weighted average per output bin; NaN where undefined.
375 """
376 n = len(dt)
377 n_out = len(tq) - 1
378 qb = q - k * np.where(k > 0, apv, 0.0) # the where avoids 0 * inf in the unbounded (apv = inf) routing
380 arrivals = _arrival_times(t=t, dt=dt, k=k, q=q, apv=apv)
381 bp = np.unique(np.concatenate([t, tq[(tq >= t[0]) & (tq <= t[-1])], arrivals[np.isfinite(arrivals)]]))
382 mids = 0.5 * (bp[:-1] + bp[1:])
383 s_bp, v0_bp = _backward_entries(queries=bp, t=t, dt=dt, k=k, q=q, apv=apv)
384 s_mid, _ = _backward_entries(queries=mids, t=t, dt=dt, k=k, q=q, apv=apv)
386 # The entry-time map s*(t) is discontinuous at the arrival of a parcel
387 # released at an edge where boundary inflow resumes after an expulsion or
388 # stagnation episode (the skipped span is the lost window), and the walk at
389 # exactly that arrival resolves to one branch by floating-point luck. Both
390 # one-sided limits are computed robustly instead: the right limit is the
391 # release edge t[g] itself; the left limit is the release time of the
392 # grazing trajectory, traced backward from (t[g], V=apv).
393 g_idx = np.nonzero(np.isfinite(arrivals))[0]
394 s_left = np.full(len(arrivals), np.nan)
395 v0_left = np.full(len(arrivals), np.nan)
396 g_pos = g_idx[g_idx >= 1] # the g = 0 arrival grazes the boundary exactly at t[0] (landing position apv)
397 if g_pos.size:
398 s_left[g_pos], v0_left[g_pos] = _backward_entries(
399 queries=t[g_pos], t=t, dt=dt, k=k, q=q, apv=apv, v_start=apv, edge_side="left"
400 )
401 av = arrivals[g_idx]
402 order = np.argsort(av)
403 av, ae = av[order], g_idx[order]
405 def arrival_edge(x):
406 if av.size == 0:
407 return np.full(len(x), -1)
408 pos = np.minimum(np.searchsorted(av, x), av.size - 1)
409 return np.where(av[pos] == x, ae[pos], -1)
411 t1, t2 = bp[:-1], bp[1:]
412 span = t2 - t1
413 m = np.clip(np.searchsorted(t, mids, side="right") - 1, 0, n - 1)
414 kc = np.searchsorted(tq, mids, side="right") - 1
415 in_out = (kc >= 0) & (kc < n_out)
416 kc = np.clip(kc, 0, n_out - 1)
418 pre = np.isnan(s_mid)
419 js = np.clip(np.searchsorted(t, np.where(pre, t[0], s_mid), side="right") - 1, 0, n - 1)
420 # Entry times at piece endpoints: one-sided limits at arrival breakpoints,
421 # the walked values elsewhere, NaN (grazing/pre-record endpoints) falling
422 # back to the entry-bin edge; the final clip into the piece's own entry bin
423 # [t[js], t[js+1]] is a roundoff guard.
424 e1, e2 = arrival_edge(t1), arrival_edge(t2)
425 s1_raw = np.where(e1 >= 0, t[np.maximum(e1, 0)], s_bp[:-1])
426 s2_raw = np.where(e2 >= 0, s_left[np.maximum(e2, 0)], s_bp[1:])
427 s_lo, s_hi = t[js], t[js + 1]
428 s1 = np.clip(np.where(np.isnan(s1_raw), s_lo, s1_raw), s_lo, s_hi)
429 s2 = np.clip(np.where(np.isnan(s2_raw), s_hi, s2_raw), s_lo, s_hi)
430 ds = np.maximum(s2 - s1, 0.0)
432 u1t = u[m] + k[m] * (t1 - t[m])
433 u2t = u[m] + k[m] * (t2 - t[m])
434 kpos = k[m] > 0
435 vol = q[m] * span
436 fac = np.where(kpos, q[m] / np.where(kpos, k[m], 1.0), vol)
438 def piece_integral(x):
439 """Integrate ``q e^{x - u(t)}`` over each piece in closed form (callers keep x <= u(t1) bounded).
441 Returns
442 -------
443 ndarray
444 One integral per piece.
445 """
446 return np.where(kpos, fac * (np.exp(x - u1t) - np.exp(x - u2t)), fac * np.exp(x - u2t))
448 # Kernel mass: full-bin sum for j in [lo, m), then the current-bin and
449 # entry-bin partial corrections (telescopes to vol when cr == cb == 1).
450 lo = np.where(pre, 0, js)
451 # The per-bin weight magnitude scales with exp(u[col+1] - u1t) (u1t <= u2t is
452 # the piece's smallest clock value), so the cutoff must key to u1t: keying it
453 # to u2t drops still-significant bins when one input bin flushes >_KERNEL_CUTOFF
454 # pore volumes.
455 lo_eff = np.maximum(lo, np.clip(np.searchsorted(u, u1t - _KERNEL_CUTOFF, side="right") - 1, 0, None))
456 w_max = int((m - lo_eff).max(initial=0))
457 ker = np.zeros(len(mids))
458 if w_max > 0:
459 cols = lo_eff[:, None] + np.arange(w_max)[None, :]
460 valid = cols < m[:, None]
461 colsc = np.clip(cols, 0, n - 1)
462 # Each bin weight is an adjacent difference of exp values at its two u
463 # edges; consecutive bins share an edge, so evaluating exp once per edge
464 # over the extended [lo_eff, m] edge array halves the exp count.
465 edges = np.clip(lo_eff[:, None] + np.arange(w_max + 1)[None, :], 0, n)
466 e_lo = np.exp(u[edges] - u1t[:, None])
467 e_hi = np.exp(u[edges] - u2t[:, None])
468 w_lo = e_lo[:, 1:] - e_lo[:, :-1]
469 w_hi = e_hi[:, 1:] - e_hi[:, :-1]
470 wgt = np.where(kpos[:, None], w_lo - w_hi, w_hi)
471 ker = np.einsum("pw,pw->p", np.where(valid, wgt, 0.0), cr[colsc])
472 mass = cr[m] * vol - cr[m] * piece_integral(u[m]) + ker * fac
474 # Boundary atom (entered) or steady-profile spin-up atom (pre-record). The
475 # where keeps -inf * 0 (unbounded routing, pre-record pieces) out of the
476 # discarded branch.
477 entered_atom = (cb[js] - cr[js]) * np.where(pre, 0.0, qb[js]) * ds + cr[js] * piece_integral(u[js])
478 qb0 = qb[0]
479 if qb0 > 0 and k[0] > 0:
480 # Landing positions v0 = G(t) at the piece endpoints. At the
481 # pre-record transition the one-sided limit is the landing position of
482 # the grazing continuation from (t[g], apv): apv itself only when the
483 # transition parcel is the t[0] release; when earlier releases were
484 # expelled (delayed transition) the grazing path lands INSIDE at t0.
485 v0_arr = np.where(np.isnan(v0_left), apv, v0_left)
486 v0_1 = np.where(e1 >= 0, v0_arr[np.maximum(e1, 0)], v0_bp[:-1])
487 v0_2 = np.where(e2 >= 0, v0_arr[np.maximum(e2, 0)], v0_bp[1:])
488 v0_1 = np.where(pre, np.where(np.isnan(v0_1), apv, v0_1), 0.0)
489 v0_2 = np.where(pre, np.where(np.isnan(v0_2), apv, v0_2), 0.0)
490 v_r0 = q[0] / k[0]
491 ic_log = np.log((v_r0 - v0_1) / (v_r0 - v0_2))
492 ic_atom = cr[0] * piece_integral(0.0) + (cb[0] - cr[0]) * (v_r0 - apv) * ic_log
493 elif qb0 > 0:
494 ic_atom = cb[0] * piece_integral(0.0) # piston pre-record: domain full of boundary water
495 else:
496 ic_atom = cr[0] * piece_integral(0.0) # boundary never fed the domain before t0
497 mass += np.where(pre, ic_atom, entered_atom)
499 take = in_out & (span > 0)
500 masses = np.zeros(n_out)
501 vols = np.zeros(n_out)
502 np.add.at(masses, kc[take], mass[take])
503 np.add.at(vols, kc[take], vol[take])
504 with np.errstate(invalid="ignore", divide="ignore"):
505 return np.where(covered & (vols > 0), masses / vols, np.nan)