Coverage for src / gwtransport / logremoval.py: 100%
37 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-04 18:51 +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).
48Available functions:
50- :func:`residence_time_to_log_removal` - Calculate log removal from residence times and
51 decay rate coefficient. Uses formula: Log Removal = log10_decay_rate * residence_time.
52 Handles single values, 1D arrays, or multi-dimensional arrays of residence times. Returns
53 log removal values with same shape as input.
55- :func:`decay_rate_to_log10_decay_rate` - Convert a natural-log decay rate constant
56 lambda [1/day] to a log10 decay rate mu [log10/day].
58- :func:`log10_decay_rate_to_decay_rate` - Convert a log10 decay rate mu [log10/day]
59 to a natural-log decay rate constant lambda [1/day].
61- :func:`parallel_mean` - Calculate weighted average log removal for parallel flow systems.
62 Computes overall efficiency when multiple treatment paths operate in parallel with different
63 log removal values and flow fractions. Uses formula: Total Log Removal = -log10(sum(F_i * 10^(-LR_i)))
64 where F_i is flow fraction and LR_i is log removal for path i. Supports multi-dimensional
65 arrays via axis parameter for batch processing. Assumes equal flow distribution if flow_fractions
66 not provided.
68- :func:`gamma_pdf` - Compute probability density function (PDF) of log removal given
69 gamma-distributed residence time. Since R = mu*T and T ~ Gamma(alpha, beta), R follows a
70 Gamma(alpha, mu*beta) distribution.
72- :func:`gamma_cdf` - Compute cumulative distribution function (CDF) of log removal given
73 gamma-distributed residence time. Returns probability that log removal is less than or equal
74 to specified values.
76- :func:`gamma_mean` - Compute effective (parallel) mean log removal for gamma-distributed
77 residence time. Uses the moment generating function of the gamma distribution to compute the
78 log-weighted average: LR_eff = alpha * log10(1 + beta * mu * ln(10)).
80- :func:`gamma_find_flow_for_target_mean` - Find flow rate that produces specified target
81 effective mean log removal given gamma-distributed aquifer pore volume. Solves inverse problem:
82 flow = apv_beta * mu * ln(10) / (10^(target_mean / apv_alpha) - 1).
84This file is part of gwtransport which is released under AGPL-3.0 license.
85See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
86"""
88import numpy as np
89import numpy.typing as npt
90from scipy import stats
93def residence_time_to_log_removal(
94 *, residence_times: npt.ArrayLike, log10_decay_rate: float
95) -> npt.NDArray[np.floating]:
96 """
97 Compute log removal given residence times and a log10 decay rate.
99 This function calculates the log removal based on residence times
100 and a log10 decay rate coefficient using first-order decay:
102 Log Removal = log10_decay_rate * residence_time
104 This corresponds to exponential decay of pathogen concentration:
105 C_out/C_in = 10^(-log10_decay_rate * residence_time).
107 Parameters
108 ----------
109 residence_times : array-like
110 Array of residence times in days. Must be positive values.
111 log10_decay_rate : float
112 Log10 decay rate coefficient (log10/day). Relates residence time
113 to log removal efficiency via first-order decay.
115 Returns
116 -------
117 log_removals : numpy.ndarray
118 Array of log removal values corresponding to the input residence times.
119 Same shape as input residence_times.
121 See Also
122 --------
123 decay_rate_to_log10_decay_rate : Convert natural-log decay rate to log10 decay rate
124 log10_decay_rate_to_decay_rate : Convert log10 decay rate to natural-log decay rate
125 gamma_mean : Compute mean log removal for gamma-distributed residence times
126 gamma_find_flow_for_target_mean : Find flow rate to achieve target log removal
127 parallel_mean : Calculate weighted average for parallel flow systems
128 gwtransport.residence_time.residence_time : Compute residence times from flow and pore volume
129 :ref:`concept-residence-time` : Time in aquifer determines pathogen contact time
131 Notes
132 -----
133 Log removal is a logarithmic measure of pathogen reduction:
134 - Log 1 = 90% reduction
135 - Log 2 = 99% reduction
136 - Log 3 = 99.9% reduction
138 The first-order decay model is mathematically identical to radioactive
139 decay used in tracer dating. To convert a published natural-log decay
140 rate lambda [1/day] to log10_decay_rate mu [log10/day], use
141 :func:`decay_rate_to_log10_decay_rate`.
143 Examples
144 --------
145 >>> import numpy as np
146 >>> from gwtransport.logremoval import residence_time_to_log_removal
147 >>> residence_times = np.array([10.0, 20.0, 50.0])
148 >>> log10_decay_rate = 0.2
149 >>> residence_time_to_log_removal(
150 ... residence_times=residence_times, log10_decay_rate=log10_decay_rate
151 ... ) # doctest: +NORMALIZE_WHITESPACE
152 array([ 2., 4., 10.])
154 >>> # Single residence time
155 >>> residence_time_to_log_removal(residence_times=5.0, log10_decay_rate=0.3)
156 np.float64(1.5)
158 >>> # 2D array of residence times
159 >>> residence_times_2d = np.array([[10.0, 20.0], [30.0, 40.0]])
160 >>> residence_time_to_log_removal(
161 ... residence_times=residence_times_2d, log10_decay_rate=0.1
162 ... )
163 array([[1., 2.],
164 [3., 4.]])
165 """
166 return log10_decay_rate * np.asarray(residence_times, dtype=float)
169def decay_rate_to_log10_decay_rate(decay_rate: float) -> float:
170 """
171 Convert a natural-log decay rate constant to a log10 decay rate.
173 Converts lambda [1/day] to mu [log10/day] using the relationship
174 mu = lambda / ln(10).
176 Parameters
177 ----------
178 decay_rate : float
179 Natural-log first-order decay rate constant lambda (1/day).
180 For example, from tracer dating: lambda = ln(2) / half_life.
182 Returns
183 -------
184 log10_decay_rate : float
185 Log10 decay rate mu (log10/day).
187 See Also
188 --------
189 log10_decay_rate_to_decay_rate : Inverse conversion
190 residence_time_to_log_removal : Apply the log10 decay rate
192 Examples
193 --------
194 >>> from gwtransport.logremoval import decay_rate_to_log10_decay_rate
195 >>> import numpy as np
196 >>> # Convert a decay rate of ln(2)/30 (half-life of 30 days)
197 >>> decay_rate = np.log(2) / 30
198 >>> decay_rate_to_log10_decay_rate(decay_rate) # doctest: +SKIP
199 0.01003...
200 """
201 return decay_rate / np.log(10)
204def log10_decay_rate_to_decay_rate(log10_decay_rate: float) -> float:
205 """
206 Convert a log10 decay rate to a natural-log decay rate constant.
208 Converts mu [log10/day] to lambda [1/day] using the relationship
209 lambda = mu * ln(10).
211 Parameters
212 ----------
213 log10_decay_rate : float
214 Log10 decay rate mu (log10/day).
216 Returns
217 -------
218 decay_rate : float
219 Natural-log first-order decay rate constant lambda (1/day).
221 See Also
222 --------
223 decay_rate_to_log10_decay_rate : Inverse conversion
225 Examples
226 --------
227 >>> from gwtransport.logremoval import log10_decay_rate_to_decay_rate
228 >>> log10_decay_rate_to_decay_rate(0.2) # doctest: +SKIP
229 0.4605...
230 """
231 return log10_decay_rate * np.log(10)
234def parallel_mean(
235 *, log_removals: npt.ArrayLike, flow_fractions: npt.ArrayLike | None = None, axis: int | None = None
236) -> np.floating | npt.NDArray[np.floating]:
237 """
238 Calculate the weighted average log removal for a system with parallel flows.
240 This function computes the overall log removal efficiency of a parallel
241 filtration system. If flow_fractions is not provided, it assumes equal
242 distribution of flow across all paths.
244 The calculation uses the formula:
246 Total Log Removal = -log10(sum(F_i * 10^(-LR_i)))
248 Where:
249 - F_i = fraction of flow through system i (decimal, sum to 1.0)
250 - LR_i = log removal of system i
252 Parameters
253 ----------
254 log_removals : array-like
255 Array of log removal values for each parallel flow.
256 Each value represents the log10 reduction of pathogens.
257 For multi-dimensional arrays, the parallel mean is computed along
258 the specified axis.
260 flow_fractions : array-like, optional
261 Array of flow fractions for each parallel flow.
262 Must sum to 1.0 along the specified axis and have compatible shape
263 with log_removals. If None, equal flow distribution is assumed
264 (default is None).
266 axis : int, optional
267 Axis along which to compute the parallel mean for multi-dimensional
268 arrays. If None, the array is treated as 1-dimensional
269 (default is None).
271 Returns
272 -------
273 float or array-like
274 The combined log removal value for the parallel system.
275 If log_removals is multi-dimensional and axis is specified,
276 returns an array with the specified axis removed.
278 Notes
279 -----
280 Log removal is a logarithmic measure of pathogen reduction:
282 - Log 1 = 90% reduction
283 - Log 2 = 99% reduction
284 - Log 3 = 99.9% reduction
286 For parallel flows, the combined removal is typically less effective
287 than the best individual removal but better than the worst.
288 For systems in series, log removals would be summed directly.
290 Examples
291 --------
292 >>> import numpy as np
293 >>> from gwtransport.logremoval import parallel_mean
294 >>> # Three parallel streams with equal flow and log removals of 3, 4, and 5
295 >>> log_removals = np.array([3, 4, 5])
296 >>> parallel_mean(log_removals=log_removals)
297 np.float64(3.431798275933005)
299 >>> # Two parallel streams with weighted flow
300 >>> log_removals = np.array([3, 5])
301 >>> flow_fractions = np.array([0.7, 0.3])
302 >>> parallel_mean(log_removals=log_removals, flow_fractions=flow_fractions)
303 np.float64(3.153044674980176)
305 >>> # Multi-dimensional array: parallel mean along axis 1
306 >>> log_removals_2d = np.array([[3, 4, 5], [2, 3, 4]])
307 >>> parallel_mean(log_removals=log_removals_2d, axis=1)
308 array([3.43179828, 2.43179828])
310 See Also
311 --------
312 residence_time_to_log_removal : Compute log removal from residence times
313 """
314 # Convert log_removals to numpy array if it isn't already
315 log_removals = np.asarray(log_removals, dtype=float)
317 # If flow_fractions is not provided, assume equal distribution
318 if flow_fractions is None:
319 if axis is None:
320 # 1D case: calculate the number of parallel flows
321 n = len(log_removals)
322 # Create equal flow fractions (avoid division by zero)
323 flow_fractions = np.full(n, 1.0 / n) if n > 0 else np.array([])
324 else:
325 # Multi-dimensional case: create equal flow fractions along the specified axis
326 n = log_removals.shape[axis]
327 shape = [1] * log_removals.ndim
328 shape[axis] = n
329 flow_fractions = np.full(shape, 1.0 / n)
330 else:
331 # Convert flow_fractions to numpy array
332 flow_fractions = np.asarray(flow_fractions, dtype=float)
334 fraction_sums = np.sum(flow_fractions, axis=axis)
335 if not np.all(np.isclose(fraction_sums, 1.0)):
336 msg = "flow_fractions must sum to 1.0 (along the specified axis)"
337 raise ValueError(msg)
339 # Convert log removal to decimal reduction values
340 decimal_reductions = 10 ** (-log_removals)
342 # Calculate weighted average decimal reduction
343 weighted_decimal_reduction = np.sum(flow_fractions * decimal_reductions, axis=axis)
345 # Convert back to log scale
346 return -np.log10(weighted_decimal_reduction)
349def gamma_pdf(
350 *, r: npt.ArrayLike, rt_alpha: float, rt_beta: float, log10_decay_rate: float
351) -> npt.NDArray[np.floating]:
352 """
353 Compute the PDF of log removal given gamma-distributed residence time.
355 Since log removal R = mu*T and T ~ Gamma(alpha, beta), the log removal R
356 follows a Gamma(alpha, mu*beta) distribution.
358 Parameters
359 ----------
360 r : array-like
361 Log removal values at which to compute the PDF.
362 rt_alpha : float
363 Shape parameter of the gamma distribution for residence time.
364 rt_beta : float
365 Scale parameter of the gamma distribution for residence time (days).
366 log10_decay_rate : float
367 Log10 decay rate mu (log10/day). Relates residence time to
368 log removal via R = mu * T.
370 Returns
371 -------
372 pdf : numpy.ndarray
373 PDF values corresponding to the input r values.
375 See Also
376 --------
377 gamma_cdf : Cumulative distribution function of log removal
378 gamma_mean : Mean of the log removal distribution
379 """
380 r = np.asarray(r)
381 return stats.gamma.pdf(r, a=rt_alpha, scale=log10_decay_rate * rt_beta)
384def gamma_cdf(
385 *, r: npt.ArrayLike, rt_alpha: float, rt_beta: float, log10_decay_rate: float
386) -> npt.NDArray[np.floating]:
387 """
388 Compute the CDF of log removal given gamma-distributed residence time.
390 Since log removal R = mu*T and T ~ Gamma(alpha, beta), the CDF is
391 P(R <= r) = P(T <= r/mu) = Gamma_CDF(r/mu; alpha, beta).
393 Parameters
394 ----------
395 r : array-like
396 Log removal values at which to compute the CDF.
397 rt_alpha : float
398 Shape parameter of the gamma distribution for residence time.
399 rt_beta : float
400 Scale parameter of the gamma distribution for residence time (days).
401 log10_decay_rate : float
402 Log10 decay rate mu (log10/day). Relates residence time to
403 log removal via R = mu * T.
405 Returns
406 -------
407 cdf : numpy.ndarray
408 CDF values corresponding to the input r values.
410 See Also
411 --------
412 gamma_pdf : Probability density function of log removal
413 gamma_mean : Mean of the log removal distribution
414 """
415 r = np.asarray(r)
416 return stats.gamma.cdf(r, a=rt_alpha, scale=log10_decay_rate * rt_beta)
419def gamma_mean(*, rt_alpha: float, rt_beta: float, log10_decay_rate: float) -> float:
420 """
421 Compute the effective (parallel) mean log removal for gamma-distributed residence time.
423 When water travels through multiple flow paths with gamma-distributed
424 residence times, the effective log removal is determined by mixing the
425 output concentrations (not by averaging individual log removals). This
426 uses the moment generating function of the gamma distribution:
428 LR_eff = -log10(E[10^(-mu*T)])
429 = alpha * log10(1 + beta * mu * ln(10))
431 This is always less than the arithmetic mean (mu * alpha * beta)
432 because short residence time paths contribute disproportionately
433 to the output concentration.
435 Parameters
436 ----------
437 rt_alpha : float
438 Shape parameter of the gamma distribution for residence time.
439 rt_beta : float
440 Scale parameter of the gamma distribution for residence time (days).
441 log10_decay_rate : float
442 Log10 decay rate mu (log10/day).
444 Returns
445 -------
446 mean : float
447 Effective (parallel) mean log removal value.
449 See Also
450 --------
451 gamma_find_flow_for_target_mean : Find flow for target mean log removal
452 parallel_mean : Discrete version of this calculation
453 gamma_pdf : PDF of the log removal distribution
454 gamma_cdf : CDF of the log removal distribution
455 :ref:`concept-pore-volume-distribution` : Why residence times are distributed
456 """
457 return rt_alpha * np.log10(1 + rt_beta * log10_decay_rate * np.log(10))
460def gamma_find_flow_for_target_mean(
461 *, target_mean: float, apv_alpha: float, apv_beta: float, log10_decay_rate: float
462) -> float:
463 """
464 Find the flow rate that produces a target effective mean log removal.
466 Given a gamma-distributed aquifer pore volume with parameters (apv_alpha, apv_beta),
467 the residence time distribution is Gamma(apv_alpha, apv_beta/flow). The effective
468 mean log removal (from :func:`gamma_mean`) is:
470 LR_eff = apv_alpha * log10(1 + (apv_beta/flow) * mu * ln(10))
472 Solving for flow:
474 flow = apv_beta * mu * ln(10) / (10^(target_mean / apv_alpha) - 1)
476 Parameters
477 ----------
478 target_mean : float
479 Target effective mean log removal value.
480 apv_alpha : float
481 Shape parameter of the gamma distribution for aquifer pore volume.
482 apv_beta : float
483 Scale parameter of the gamma distribution for aquifer pore volume.
484 log10_decay_rate : float
485 Log10 decay rate mu (log10/day).
487 Returns
488 -------
489 flow : float
490 Flow rate (same units as apv_beta per day) that produces the
491 target mean log removal.
493 See Also
494 --------
495 gamma_mean : Compute effective mean log removal for given parameters
496 """
497 return apv_beta * log10_decay_rate * np.log(10) / (10 ** (target_mean / apv_alpha) - 1)