Coverage for src/gwtransport/logremoval.py: 96%
46 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1"""
2Log Removal Calculations for First-Order Decay Processes.
4This module provides utilities to calculate log removal values from first-order decay
5processes, including pathogen inactivation and radioactive decay. The module supports
6basic log removal calculations and parallel flow arrangements where multiple flow paths
7operate simultaneously.
9First-Order Decay Model
10-----------------------
11The log removal from any first-order decay process is:
13 Log Removal = log10_decay_rate * residence_time
15where ``log10_decay_rate`` has units [log10/day] and ``residence_time`` has units [days].
16This is equivalent to exponential decay ``C_out/C_in = 10^(-mu * t)``, where mu is the
17log10 decay rate and t is residence time. The natural-log decay rate constant lambda [1/day]
18is related to mu by ``lambda = mu * ln(10)``.
20This model applies to any process that follows first-order kinetics:
22- **Pathogen inactivation**: viruses, bacteria, and protozoa lose infectivity over time
23- **Radioactive decay**: isotopes used for groundwater dating (tritium, CFC, SF6)
24- **Chemical degradation**: first-order breakdown of contaminants
26Pathogen Removal in Bank Filtration
27------------------------------------
28For pathogen removal during soil passage, total removal consists of two distinct mechanisms
29(Schijven and Hassanizadeh, 2000):
311. **Inactivation (time-dependent)**: Pathogens lose infectivity over time through biological
32 decay. This follows first-order kinetics and is modeled by this module as
33 ``LR_decay = log10_decay_rate * residence_time``. The inactivation rate depends strongly
34 on temperature and pathogen type.
362. **Attachment (geometry-dependent)**: Pathogens are physically removed by adsorption to soil
37 grains and straining. This depends on aquifer geometry, distance, soil properties, and pH,
38 and is NOT modeled by this module. Users should add this component separately based on
39 site-specific data.
41Total log removal = LR_decay (this module) + LR_attachment (user-specified).
43At the Castricum dune recharge site, Schijven et al. (1999) found that attachment contributed
44approximately 97% of total MS2 removal, with inactivation contributing only 3%. Inactivation
45rates for common model viruses at 10 degrees C are typically 0.02-0.11 log10/day (Schijven and
46Hassanizadeh, 2000, Table 7).
48Gamma-distribution parameter notation
49-------------------------------------
50Several functions are parameterized by a gamma distribution. The parameter prefix marks
51*which* physical quantity is gamma-distributed, because two distinct quantities appear here:
53- ``rt_alpha`` / ``rt_beta`` / ``rt_loc`` (or the equivalent ``rt_mean`` / ``rt_std`` /
54 ``rt_loc``) parameterize the gamma distribution of the **residence time** (used by
55 :func:`gamma_pdf`, :func:`gamma_cdf`, :func:`gamma_mean`).
56- ``apv_alpha`` / ``apv_beta`` / ``apv_loc`` (or the equivalent ``apv_mean`` / ``apv_std`` /
57 ``apv_loc``) parameterize the gamma distribution of the **aquifer pore volume** (used by
58 :func:`gamma_find_flow_for_target_mean`).
60These prefixes are intentional and load-bearing: residence time and pore volume are different
61quantities, so a bare ``alpha`` / ``beta`` / ``loc`` would be ambiguous in this module. Both the
62shape/scale and the mean/std pairs are validated through :func:`gwtransport.gamma.parse_parameters`,
63so invalid parameters (e.g. a negative shape) raise ``ValueError`` rather than silently returning
64an unphysical result.
66Available functions:
68- :func:`residence_time_to_log_removal` - Calculate log removal from residence times and
69 decay rate coefficient. Uses formula: Log Removal = log10_decay_rate * residence_time.
70 Handles single values, 1D arrays, or multi-dimensional arrays of residence times. Returns
71 log removal values with same shape as input.
73- :func:`decay_rate_to_log10_decay_rate` - Convert a natural-log decay rate constant
74 lambda [1/day] to a log10 decay rate mu [log10/day].
76- :func:`log10_decay_rate_to_decay_rate` - Convert a log10 decay rate mu [log10/day]
77 to a natural-log decay rate constant lambda [1/day].
79- :func:`parallel_mean` - Calculate weighted average log removal for parallel flow systems.
80 Computes overall efficiency when multiple treatment paths operate in parallel with different
81 log removal values and flow fractions. Uses formula: Total Log Removal = -log10(sum(F_i * 10^(-LR_i)))
82 where F_i is flow fraction and LR_i is log removal for path i. Supports multi-dimensional
83 arrays via axis parameter for batch processing. Assumes equal flow distribution if flow_fractions
84 not provided.
86- :func:`gamma_pdf` - Compute probability density function (PDF) of log removal given
87 gamma-distributed residence time. Since R = mu*T and T ~ Gamma(alpha, beta), R follows a
88 Gamma(alpha, mu*beta) distribution.
90- :func:`gamma_cdf` - Compute cumulative distribution function (CDF) of log removal given
91 gamma-distributed residence time. Returns probability that log removal is less than or equal
92 to specified values.
94- :func:`gamma_mean` - Compute effective (parallel) mean log removal for gamma-distributed
95 residence time. Uses the moment generating function of the gamma distribution to compute the
96 log-weighted average: LR_eff = mu * loc + alpha * log10(1 + beta * mu * ln(10)).
98- :func:`gamma_find_flow_for_target_mean` - Find flow rate that produces specified target
99 effective mean log removal given gamma-distributed aquifer pore volume. For ``loc == 0`` this
100 is the closed-form inverse: flow = beta * mu * ln(10) / (10^(target_mean / alpha) - 1); for
101 ``loc > 0`` the transcendental equation is solved numerically.
103This file is part of gwtransport which is released under AGPL-3.0 license.
104See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
105"""
107import numpy as np
108import numpy.typing as npt
109from scipy import optimize, stats
111from gwtransport.gamma import parse_parameters
114def residence_time_to_log_removal(
115 *, residence_times: npt.ArrayLike, log10_decay_rate: float
116) -> npt.NDArray[np.floating]:
117 """
118 Compute log removal given residence times and a log10 decay rate.
120 This function calculates the log removal based on residence times
121 and a log10 decay rate coefficient using first-order decay:
123 Log Removal = log10_decay_rate * residence_time
125 This corresponds to exponential decay of pathogen concentration:
126 C_out/C_in = 10^(-log10_decay_rate * residence_time).
128 Parameters
129 ----------
130 residence_times : array-like
131 Residence times in days. The formula evaluates ``log10_decay_rate *
132 residence_times`` for any real input; negative values produce negative
133 log removal (mathematical amplification) and the caller is responsible
134 for sign interpretation.
135 log10_decay_rate : float
136 Log10 decay rate coefficient (log10/day). Relates residence time
137 to log removal efficiency via first-order decay. Negative values
138 correspond to first-order production rather than decay.
140 Returns
141 -------
142 log_removals : ndarray
143 Array of log removal values corresponding to the input residence times.
144 Same shape as input residence_times.
146 See Also
147 --------
148 decay_rate_to_log10_decay_rate : Convert natural-log decay rate to log10 decay rate
149 log10_decay_rate_to_decay_rate : Convert log10 decay rate to natural-log decay rate
150 gamma_mean : Compute mean log removal for gamma-distributed residence times
151 gamma_find_flow_for_target_mean : Find flow rate to achieve target log removal
152 parallel_mean : Calculate weighted average for parallel flow systems
153 gwtransport.residence_time.full : Compute residence times from flow and pore volume
154 :ref:`concept-residence-time` : Time in aquifer determines pathogen contact time
156 Notes
157 -----
158 Log removal is a logarithmic measure of pathogen reduction:
159 - Log 1 = 90% reduction
160 - Log 2 = 99% reduction
161 - Log 3 = 99.9% reduction
163 The first-order decay model is mathematically identical to radioactive
164 decay used in tracer dating. To convert a published natural-log decay
165 rate lambda [1/day] to log10_decay_rate mu [log10/day], use
166 :func:`decay_rate_to_log10_decay_rate`.
168 Examples
169 --------
170 >>> import numpy as np
171 >>> from gwtransport.logremoval import residence_time_to_log_removal
172 >>> residence_times = np.array([10.0, 20.0, 50.0])
173 >>> log10_decay_rate = 0.2
174 >>> residence_time_to_log_removal(
175 ... residence_times=residence_times, log10_decay_rate=log10_decay_rate
176 ... ) # doctest: +NORMALIZE_WHITESPACE
177 array([ 2., 4., 10.])
179 >>> # Single residence time
180 >>> residence_time_to_log_removal(residence_times=5.0, log10_decay_rate=0.3)
181 np.float64(1.5)
183 >>> # 2D array of residence times
184 >>> residence_times_2d = np.array([[10.0, 20.0], [30.0, 40.0]])
185 >>> residence_time_to_log_removal(
186 ... residence_times=residence_times_2d, log10_decay_rate=0.1
187 ... )
188 array([[1., 2.],
189 [3., 4.]])
190 """
191 return log10_decay_rate * np.asarray(residence_times, dtype=float)
194def decay_rate_to_log10_decay_rate(decay_rate: float) -> float:
195 """
196 Convert a natural-log decay rate constant to a log10 decay rate.
198 Converts lambda [1/day] to mu [log10/day] using the relationship
199 mu = lambda / ln(10).
201 Parameters
202 ----------
203 decay_rate : float
204 Natural-log first-order decay rate constant lambda (1/day).
205 For example, from tracer dating: lambda = ln(2) / half_life.
207 Returns
208 -------
209 log10_decay_rate : float
210 Log10 decay rate mu (log10/day).
212 See Also
213 --------
214 log10_decay_rate_to_decay_rate : Inverse conversion
215 residence_time_to_log_removal : Apply the log10 decay rate
217 Examples
218 --------
219 >>> from gwtransport.logremoval import decay_rate_to_log10_decay_rate
220 >>> import numpy as np
221 >>> # Convert a decay rate of ln(2)/30 (half-life of 30 days)
222 >>> decay_rate = np.log(2) / 30
223 >>> decay_rate_to_log10_decay_rate(decay_rate) # doctest: +ELLIPSIS
224 np.float64(0.01003...)
225 """
226 return decay_rate / np.log(10)
229def log10_decay_rate_to_decay_rate(log10_decay_rate: float) -> float:
230 """
231 Convert a log10 decay rate to a natural-log decay rate constant.
233 Converts mu [log10/day] to lambda [1/day] using the relationship
234 lambda = mu * ln(10).
236 Parameters
237 ----------
238 log10_decay_rate : float
239 Log10 decay rate mu (log10/day).
241 Returns
242 -------
243 decay_rate : float
244 Natural-log first-order decay rate constant lambda (1/day).
246 See Also
247 --------
248 decay_rate_to_log10_decay_rate : Inverse conversion
250 Examples
251 --------
252 >>> from gwtransport.logremoval import log10_decay_rate_to_decay_rate
253 >>> log10_decay_rate_to_decay_rate(0.2) # doctest: +ELLIPSIS
254 np.float64(0.4605...)
255 """
256 return log10_decay_rate * np.log(10)
259def parallel_mean(
260 *, log_removals: npt.ArrayLike, flow_fractions: npt.ArrayLike | None = None, axis: int | None = None
261) -> np.floating | npt.NDArray[np.floating]:
262 """
263 Calculate the weighted average log removal for a system with parallel flows.
265 This function computes the overall log removal efficiency of a parallel
266 filtration system. If flow_fractions is not provided, it assumes equal
267 distribution of flow across all paths.
269 The calculation uses the formula:
271 Total Log Removal = -log10(sum(F_i * 10^(-LR_i)))
273 Where:
274 - F_i = fraction of flow through system i (decimal, sum to 1.0)
275 - LR_i = log removal of system i
277 Parameters
278 ----------
279 log_removals : array-like
280 Array of log removal values for each parallel flow.
281 Each value represents the log10 reduction of pathogens.
282 For multi-dimensional arrays, the parallel mean is computed along
283 the specified axis.
285 flow_fractions : array-like, optional
286 Array of flow fractions for each parallel flow.
287 Must sum to 1.0 along the specified axis and have compatible shape
288 with log_removals. If None, equal flow distribution is assumed
289 (default is None).
291 axis : int, optional
292 Axis along which to compute the parallel mean for multi-dimensional
293 arrays. If None, the reduction matches the way ``np.mean`` / ``np.sum``
294 treat ``axis=None``: the parallel mean is computed over the flattened
295 input (default is None).
297 Returns
298 -------
299 np.floating or ndarray
300 The combined log removal value for the parallel system. Returns a
301 scalar when axis=None, otherwise an array with the specified axis
302 removed.
304 Raises
305 ------
306 ValueError
307 If ``flow_fractions`` does not sum to 1.0 along the specified axis.
309 See Also
310 --------
311 residence_time_to_log_removal : Compute log removal from residence times
313 Notes
314 -----
315 Log removal is a logarithmic measure of pathogen reduction:
317 - Log 1 = 90% reduction
318 - Log 2 = 99% reduction
319 - Log 3 = 99.9% reduction
321 For parallel flows, the combined removal is typically less effective
322 than the best individual removal but better than the worst.
323 For systems in series, log removals would be summed directly.
325 Examples
326 --------
327 >>> import numpy as np
328 >>> from gwtransport.logremoval import parallel_mean
329 >>> # Three parallel streams with equal flow and log removals of 3, 4, and 5
330 >>> log_removals = np.array([3, 4, 5])
331 >>> parallel_mean(log_removals=log_removals)
332 np.float64(3.431798275933005)
334 >>> # Two parallel streams with weighted flow
335 >>> log_removals = np.array([3, 5])
336 >>> flow_fractions = np.array([0.7, 0.3])
337 >>> parallel_mean(log_removals=log_removals, flow_fractions=flow_fractions)
338 np.float64(3.153044674980176)
340 >>> # Multi-dimensional array: parallel mean along axis 1
341 >>> log_removals_2d = np.array([[3, 4, 5], [2, 3, 4]])
342 >>> parallel_mean(log_removals=log_removals_2d, axis=1)
343 array([3.43179828, 2.43179828])
344 """
345 log_removals = np.asarray(log_removals, dtype=float)
346 decimal_reductions = 10.0 ** (-log_removals)
347 if flow_fractions is None:
348 return -np.log10(np.mean(decimal_reductions, axis=axis))
349 flow_fractions = np.asarray(flow_fractions, dtype=float)
350 if not np.all(np.isclose(np.sum(flow_fractions, axis=axis), 1.0)):
351 msg = "flow_fractions must sum to 1.0 (along the specified axis)"
352 raise ValueError(msg)
353 return -np.log10(np.sum(flow_fractions * decimal_reductions, axis=axis))
356def gamma_pdf(
357 *,
358 r: npt.ArrayLike,
359 rt_alpha: float | None = None,
360 rt_beta: float | None = None,
361 rt_loc: float = 0.0,
362 rt_mean: float | None = None,
363 rt_std: float | None = None,
364 log10_decay_rate: float,
365) -> npt.NDArray[np.floating]:
366 """
367 Compute the PDF of log removal given (shifted) gamma-distributed residence time.
369 With residence time ``T = T0 + rt_loc`` where ``T0 ~ Gamma(rt_alpha, rt_beta)``,
370 the log removal ``R = mu * T`` follows a shifted gamma distribution with shape
371 ``rt_alpha``, scale ``mu * rt_beta``, and location ``mu * rt_loc``.
373 The residence-time distribution is specified with either ``(rt_alpha, rt_beta)`` or
374 ``(rt_mean, rt_std)`` (optionally shifted by ``rt_loc``); both are routed through
375 :func:`gwtransport.gamma.parse_parameters`.
377 Parameters
378 ----------
379 r : array-like
380 Log removal values at which to compute the PDF.
381 rt_alpha : float, optional
382 Shape parameter of the gamma distribution for residence time. Must be positive.
383 rt_beta : float, optional
384 Scale parameter of the gamma distribution for residence time (days). Must be positive.
385 rt_loc : float, optional
386 Location (minimum residence time, days) of the residence time distribution.
387 Must be non-negative. Default is ``0.0``.
388 rt_mean : float, optional
389 Mean residence time (days). Alternative to ``rt_alpha``; supply with ``rt_std``.
390 Must be strictly greater than ``rt_loc``.
391 rt_std : float, optional
392 Standard deviation of the residence time (days). Alternative to ``rt_beta``;
393 supply with ``rt_mean``. Must be positive.
394 log10_decay_rate : float
395 Log10 decay rate mu (log10/day). Relates residence time to
396 log removal via R = mu * T.
398 Returns
399 -------
400 pdf : ndarray
401 PDF values corresponding to the input r values.
403 Raises
404 ------
405 ValueError
406 If parameter validation in :func:`gwtransport.gamma.parse_parameters` fails
407 (e.g. ``rt_loc`` negative, non-positive shape/scale, or neither/both
408 parameter pairs supplied).
410 See Also
411 --------
412 gamma_cdf : Cumulative distribution function of log removal
413 gamma_mean : Mean of the log removal distribution
414 """
415 rt_alpha, rt_beta, rt_loc = parse_parameters(mean=rt_mean, std=rt_std, loc=rt_loc, alpha=rt_alpha, beta=rt_beta)
416 return stats.gamma.pdf(r, a=rt_alpha, loc=log10_decay_rate * rt_loc, scale=log10_decay_rate * rt_beta)
419def gamma_cdf(
420 *,
421 r: npt.ArrayLike,
422 rt_alpha: float | None = None,
423 rt_beta: float | None = None,
424 rt_loc: float = 0.0,
425 rt_mean: float | None = None,
426 rt_std: float | None = None,
427 log10_decay_rate: float,
428) -> npt.NDArray[np.floating]:
429 """
430 Compute the CDF of log removal given (shifted) gamma-distributed residence time.
432 With residence time ``T = T0 + rt_loc`` where ``T0 ~ Gamma(rt_alpha, rt_beta)``,
433 the CDF is ``P(R <= r) = P(mu*(T0 + rt_loc) <= r) =
434 P(T0 <= (r - mu*rt_loc)/mu)`` which is the CDF of a shifted gamma distribution
435 with location ``mu * rt_loc``.
437 The residence-time distribution is specified with either ``(rt_alpha, rt_beta)`` or
438 ``(rt_mean, rt_std)`` (optionally shifted by ``rt_loc``); both are routed through
439 :func:`gwtransport.gamma.parse_parameters`.
441 Parameters
442 ----------
443 r : array-like
444 Log removal values at which to compute the CDF.
445 rt_alpha : float, optional
446 Shape parameter of the gamma distribution for residence time. Must be positive.
447 rt_beta : float, optional
448 Scale parameter of the gamma distribution for residence time (days). Must be positive.
449 rt_loc : float, optional
450 Location (minimum residence time, days) of the residence time distribution.
451 Must be non-negative. Default is ``0.0``.
452 rt_mean : float, optional
453 Mean residence time (days). Alternative to ``rt_alpha``; supply with ``rt_std``.
454 Must be strictly greater than ``rt_loc``.
455 rt_std : float, optional
456 Standard deviation of the residence time (days). Alternative to ``rt_beta``;
457 supply with ``rt_mean``. Must be positive.
458 log10_decay_rate : float
459 Log10 decay rate mu (log10/day). Relates residence time to
460 log removal via R = mu * T.
462 Returns
463 -------
464 cdf : ndarray
465 CDF values corresponding to the input r values.
467 Raises
468 ------
469 ValueError
470 If parameter validation in :func:`gwtransport.gamma.parse_parameters` fails
471 (e.g. ``rt_loc`` negative, non-positive shape/scale, or neither/both
472 parameter pairs supplied).
474 See Also
475 --------
476 gamma_pdf : Probability density function of log removal
477 gamma_mean : Mean of the log removal distribution
478 """
479 rt_alpha, rt_beta, rt_loc = parse_parameters(mean=rt_mean, std=rt_std, loc=rt_loc, alpha=rt_alpha, beta=rt_beta)
480 return stats.gamma.cdf(r, a=rt_alpha, loc=log10_decay_rate * rt_loc, scale=log10_decay_rate * rt_beta)
483def gamma_mean(
484 *,
485 rt_alpha: float | None = None,
486 rt_beta: float | None = None,
487 rt_loc: float = 0.0,
488 rt_mean: float | None = None,
489 rt_std: float | None = None,
490 log10_decay_rate: float,
491) -> float:
492 """
493 Compute the effective (parallel) mean log removal for (shifted) gamma-distributed residence time.
495 When water travels through multiple flow paths with gamma-distributed
496 residence times, the effective log removal is determined by mixing the
497 output concentrations (not by averaging individual log removals). For a
498 shifted gamma distribution ``T = T0 + rt_loc`` with ``T0 ~ Gamma(alpha, beta)``,
499 factoring the moment generating function gives:
501 LR_eff = -log10(E[10^(-mu*T)])
502 = -log10(10^(-mu*rt_loc) * E[10^(-mu*T0)])
503 = mu * rt_loc + alpha * log10(1 + beta * mu * ln(10))
505 The ``rt_loc`` term shifts the whole log-removal distribution by a constant
506 ``mu * rt_loc``; the alpha/beta term is unchanged. This is always less than
507 the arithmetic mean ``mu * (alpha * beta + rt_loc)`` because short residence
508 time paths contribute disproportionately to the output concentration.
510 The residence-time distribution is specified with either ``(rt_alpha, rt_beta)`` or
511 ``(rt_mean, rt_std)`` (optionally shifted by ``rt_loc``); both are routed through
512 :func:`gwtransport.gamma.parse_parameters`.
514 Parameters
515 ----------
516 rt_alpha : float, optional
517 Shape parameter of the gamma distribution for residence time. Must be positive.
518 rt_beta : float, optional
519 Scale parameter of the gamma distribution for residence time (days). Must be positive.
520 rt_loc : float, optional
521 Location (minimum residence time, days) of the residence time distribution.
522 Must be non-negative. Default is ``0.0``.
523 rt_mean : float, optional
524 Mean residence time (days). Alternative to ``rt_alpha``; supply with ``rt_std``.
525 Must be strictly greater than ``rt_loc``.
526 rt_std : float, optional
527 Standard deviation of the residence time (days). Alternative to ``rt_beta``;
528 supply with ``rt_mean``. Must be positive.
529 log10_decay_rate : float
530 Log10 decay rate mu (log10/day).
532 Returns
533 -------
534 mean : float
535 Effective (parallel) mean log removal value.
537 Raises
538 ------
539 ValueError
540 If parameter validation in :func:`gwtransport.gamma.parse_parameters` fails
541 (e.g. ``rt_loc`` negative, non-positive shape/scale, or neither/both
542 parameter pairs supplied).
544 See Also
545 --------
546 gamma_find_flow_for_target_mean : Find flow for target mean log removal
547 parallel_mean : Discrete version of this calculation
548 gamma_pdf : PDF of the log removal distribution
549 gamma_cdf : CDF of the log removal distribution
550 :ref:`concept-pore-volume-distribution` : Why residence times are distributed
551 """
552 rt_alpha, rt_beta, rt_loc = parse_parameters(mean=rt_mean, std=rt_std, loc=rt_loc, alpha=rt_alpha, beta=rt_beta)
553 return log10_decay_rate * rt_loc + rt_alpha * np.log1p(rt_beta * log10_decay_rate * np.log(10)) / np.log(10)
556def gamma_find_flow_for_target_mean(
557 *,
558 target_mean: float,
559 apv_alpha: float | None = None,
560 apv_beta: float | None = None,
561 apv_loc: float = 0.0,
562 apv_mean: float | None = None,
563 apv_std: float | None = None,
564 log10_decay_rate: float,
565) -> float:
566 """
567 Find the flow rate that produces a target effective mean log removal.
569 Given a (shifted) gamma-distributed aquifer pore volume with parameters
570 ``(apv_alpha, apv_beta, apv_loc)``, the residence time distribution at flow
571 ``Q`` is a shifted gamma with shape ``apv_alpha``, scale ``apv_beta/Q``, and
572 location ``apv_loc/Q``. From :func:`gamma_mean`:
574 LR_eff = mu * apv_loc / Q + apv_alpha * log10(1 + (apv_beta/Q) * mu * ln(10))
576 For ``apv_loc == 0`` this is closed-form:
578 Q = apv_beta * mu * ln(10) / (10^(target_mean / apv_alpha) - 1)
580 For ``apv_loc > 0`` the equation is transcendental and solved numerically
581 with :func:`scipy.optimize.brentq` by bracketing the root in ``1/Q``.
583 The pore-volume distribution is specified with either ``(apv_alpha, apv_beta)`` or
584 ``(apv_mean, apv_std)`` (optionally shifted by ``apv_loc``); both are routed through
585 :func:`gwtransport.gamma.parse_parameters`.
587 Parameters
588 ----------
589 target_mean : float
590 Target effective mean log removal value. Must be positive.
591 apv_alpha : float, optional
592 Shape parameter of the gamma distribution for aquifer pore volume. Must be positive.
593 apv_beta : float, optional
594 Scale parameter of the gamma distribution for aquifer pore volume. Must be positive.
595 apv_loc : float, optional
596 Location (minimum aquifer pore volume) of the gamma distribution.
597 Must be non-negative. Default is ``0.0``.
598 apv_mean : float, optional
599 Mean aquifer pore volume. Alternative to ``apv_alpha``; supply with ``apv_std``.
600 Must be strictly greater than ``apv_loc``.
601 apv_std : float, optional
602 Standard deviation of the aquifer pore volume. Alternative to ``apv_beta``;
603 supply with ``apv_mean``. Must be positive.
604 log10_decay_rate : float
605 Log10 decay rate mu (log10/day).
607 Returns
608 -------
609 flow : float
610 Flow rate (same units as apv_beta per day) that produces the
611 target mean log removal.
613 Raises
614 ------
615 ValueError
616 If ``target_mean`` is not positive, if ``log10_decay_rate`` is not positive
617 (no decay can never produce a positive target log removal), or if parameter
618 validation in :func:`gwtransport.gamma.parse_parameters` fails (e.g. ``apv_loc``
619 negative, non-positive shape/scale, or neither/both parameter pairs supplied).
621 See Also
622 --------
623 gamma_mean : Compute effective mean log removal for given parameters
624 """
625 apv_alpha, apv_beta, apv_loc = parse_parameters(
626 mean=apv_mean, std=apv_std, loc=apv_loc, alpha=apv_alpha, beta=apv_beta
627 )
628 if target_mean <= 0:
629 msg = "target_mean must be positive"
630 raise ValueError(msg)
631 if log10_decay_rate <= 0:
632 # Without decay, the effective mean log removal is identically zero
633 # regardless of flow, so no finite flow can attain a positive target.
634 msg = "log10_decay_rate must be positive to attain a positive target_mean"
635 raise ValueError(msg)
637 ln10 = np.log(10)
638 flow_closed_form = apv_beta * log10_decay_rate * ln10 / (10 ** (target_mean / apv_alpha) - 1)
640 if apv_loc == 0.0:
641 return float(flow_closed_form)
643 # Solve target = mu*apv_loc*u + apv_alpha*log10(1 + apv_beta*mu*ln(10)*u) for u = 1/flow.
644 # Both terms are monotonically increasing in u, so f(u) - target is monotonic with a
645 # unique positive root. Bracket: at u = 1/flow_closed_form the alpha/beta term alone
646 # equals target_mean, so the full f overshoots by exactly mu*apv_loc*u_upper > 0.
647 u_upper = 1.0 / flow_closed_form
649 def residual(u: float) -> float:
650 return float(
651 log10_decay_rate * apv_loc * u
652 + apv_alpha * np.log1p(apv_beta * log10_decay_rate * ln10 * u) / ln10
653 - target_mean
654 )
656 u_root = optimize.brentq(residual, 0.0, u_upper)
657 return 1.0 / float(u_root) # type: ignore[arg-type]