Environ. Sci. Technol. 2004, 38, 5925-5931
Uncertainty Analysis in a Mechanistic Model of Bacterial Regrowth in Distribution Systems† FRANCIS A. DIGIANO* Department of Environmental Sciences and Engineering, University of North Carolina, 106 Rosenau Hall, CB 7400, Chapel Hill, North Carolina 27599-7431 WEIDONG ZHANG Hazen and Sawyer, P.C., 629 Green Valley Road, Suite 210, Greensboro, North Carolina 27408
The first generation of mechanistic models of bacterial regrowth in distribution systems (DS) provides insight into cause-and-effect relationships. However, the state of knowledge about the processes included in these models is insufficient to warrant deterministic predictions. Even if the process descriptions are reasonable, the uncertainty in values of key system constants limits predictions of bacterial growth. A new mechanistic model was developed to incorporate the accepted knowledge of physical, chemical, and biological processes with the hydraulic features in order to capture the unsteady state behavior of the DS. Sensitivity testing showed that the extent of bacterial regrowth was affected mainly by the rate constants for chlorine decay reactions in bulk water and on the pipe wall and by the maximum growth rate constant of attached bacteria. A simple hypothetical network was used to evaluate the effects of uncertainty in these three system constants by running 100 Monte Carlo simulations. Cumulative probability plots showed a wide range of predictions for concentrations of bacteria and chlorine in bulk water at various nodes in the DS. The magnitude of these concentrations and the range of values were greatly affected by water residence time to each node. Once the chlorine residual is depleted, bacterial growth is mainly influenced by the amount of substrate available. However, high values of the coefficients for the maximum growth rate of attached bacteria, the chlorine decay in bulk water, and the chlorine decay by wall reaction did not necessarily lead to the maximum bacterial growth at a given location.
Introduction Despite large investments in water treatment processes to meet new federal regulations, water quality has been shown to deteriorate in water distribution systems (DS) (1-7). The reactions of most concern to public health are those that involve chemical reduction of residual disinfectant, both in the water flowing within pipes and at pipe surfaces (8). These reactions lower the disinfectant residual and thus remove inhibition to bacterial regrowth if other conditions in the DS are suitable. Water residence time (water age), pipe materials, concentration of biodegradable organic matter, temperature, and pH are all important factors that determine the rate of †
This paper is part of the Walter J. Weber Jr. tribute issue. * Corresponding author phone: (919)966-2480; fax: (919)966-1171; e-mail:
[email protected]. 10.1021/es049745l CCC: $27.50 Published on Web 10/14/2004
2004 American Chemical Society
disinfectant loss and the extent of regrowth once the disinfectant residual has been lowered. The cause-and-effect relationship for bacterial regrowth defies simple description. A few field studies have produced statistically based models (1, 2, 4, 6, 7). However, these have failed to explain regrowth completely. The spatial and temporal limitations of typical sampling programs may understandably fail to capture the highly dynamic nature of the DS comprising thousands of pipe segments of different diameter and materials along with numerous storage tanks in which water quality can change due to diurnal as well as seasonal variations in water demand. Several mechanistic models have also been put forward (9-12), but testing has been limited. Each has a similar description of microbial processes adapted from Monod kinetics and biofilm modeling of reactors used in wastewater treatment (13). However, extension to DS modeling is not straightforward because of the need to account for microbial growth in the presence of an inhibitor (disinfectant), oliogtrophic growth conditions, inactivation by the disinfectant, and bacterial detachment from pipe walls. With regard to the last issue, recent modeling has shown that detachment accounts for most of the bacteria in the bulk water in contrast to growth within the bulk water, and thus, the mechanistic description of this process is critical (12). While current mechanistic models of regrowth and disinfectant residual may be over-specified in light of existing knowledge, they still may provide a framework by which to understand the interactions among physical, chemical, and biological processes. In recognition of imperfect knowledge, the objective of this paper is to incorporate uncertainty analysis into a mechanistic model using a Monte Carlo approach. Uncertainty has been incorporated into an increasing number of deterministic environmental process models using a Monte Carlo approach (14-18). The ultimate goal of this research would be to develop a frequency distribution of key variables (e.g., bacterial number concentration and residual chlorine concentration) at locations within the DS in recognition of uncertainty in the system constants of current mechanistic models.
Mechanistic Model A new mathematical model of bacterial regrowth was established to combine hydraulic calculations with a description of free and attached growth of bacteria, detachment of attached bacteria, endogenous respiration, growth rate suppression by chlorine, and decay of chlorine in the bulk water and at the pipe wall (12). This model improves upon the SANCHO model (9) by direct inclusion of hydraulic calculations, by allowance for dynamic variations in water velocity and water quality, and by reduction in specification of system constants from 19 to 12. While the PICCOBIO model (10) couples the regrowth and hydraulic network models, it does not capture the dynamics of bacterial regrowth for the typical unsteady-state conditions of a DS. Moreover, the mechanistic descriptions are far more detailed, which increases the number of parameters to be specified. Another drawback is the proprietary nature of PICCOBIO that prevents its modification. A set of partial differential equations described the concentrations of bacteria in bulk water (free bacteria), bacteria attached to pipe surfaces (biofilm), substrate, and chlorine as a function of position within the pipe network and of time, due both to the unsteady state in growth of bacteria and changes in water quality entering the DS. Advective-dispersive transport occurs in the axial direction VOL. 38, NO. 22, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
5925
of the pipe, with dispersion being constant throughout the DS. The bacterial growth rate is determined by the organic substrate concentration, the inhibition caused by chlorine and by temperature:
{
µb )
(
) (
)
(Cl2) - (Cl2,t) S exp S + Ks Cl2,c Topt - T 2 exp Topt - Ti Topt - T 2 S µmax,b exp S + Ks Topt - Ti
µmax,b
((
)) ( ) ((
))
when Cl2 > Cl2,t
when Cl2 e Cl2,t
}
(1) where µb is the specific growth rate of free bacteria in the bulk water (t-1); µmax,b is the maximum specific growth rate of these bacteria (t-1); S is the concentration of substrate (m/l3 as C), taken as the biodegradable dissolved organic carbon (BDOC); Ks is the half-saturation constant of substrate uptake (m/l3 as C); Topt is the optimal temperature for bacterial activity (10); Ti is a temperature-dependent shape parameter (°C); T is the in-situ temperature; Cl2,t is the threshold concentration (m/l3) above which chlorine suppresses bacterial activity; and Cl2,c is the characteristic chlorine concentration (m/l3) (9). The same form of equation holds for µa, the attached specific growth rate (t-1) for description of the biofilm. Cl2,t will be less for free bacteria than attached bacteria because biofilms are more resistant to growth rate suppression by chlorine. Free bacteria can deposit on the pipe wall to become attached bacteria. Detachment is possible to form free bacteria so that the net rate of change in free bacteria (cells l-3 t-1) is given by:
(rX,b) ) - kdepXb + kdetXaν/Rh
(2)
where kdep is a first-order kinetic constant for free bacteria that deposit onto the pipe wall (t-1), Xb is the number concentration of bacteria in the bulk water (cells/l3), kdet is a kinetic constant for detachment of bacteria (l-1), v is the water velocity (l/t), and Xa is the number of bacteria attached per unit area of pipe wall (cells/l2), and Rh is hydraulic radius (l). The attached biofilm is a uniform thin layer of biomass such that diffusion of substrate and chlorine into the biofilm are not assumed to be rate-limiting. The material balances for all constituents are listed below:
∂Xb ∂2 X b ∂Xb ) -v + Dd 2 + µbXb + kdetXav/Rh ∂t ∂x ∂x kdepXb - kdXb (3) ∂Xa ) µaXa - kdetXav - kdXa + kdepXbRh ∂t
( )
∂S ∂S ∂2S 1 ) -v + Dd 2 (µ X /R + µbXb) ∂t ∂x Ygβ a a h ∂x
(4)
9
Uncertainty in any model output variable should be the sum of products of sensitivities to input factors and their associated uncertainties (variances). This analysis of uncertainty, however, was restricted by choosing to include only the three most-sensitive input factors with respect to the extent of regrowth and not considering differences in uncertainty among them. Admittedly, input factors to eqs 1-6 with less sensitivity but with larger uncertainty could be important as well. However, the published data on these factors are too scant to rank their uncertainty. Similar to a study of the PICCOBIO model (21), a single pipe line (8000 m in length and 8 in. in diameter) with a water velocity of 0.13 m/s was used to measure sensitivity as defined by
S* ) (6)
Equation 3 includes transport of bacteria by advection and dispersion, growth of bacteria (as affected by substrate, chlorine residual, and temperature), deposition to the pipe wall and release from the pipe wall, and endogenous decay. Dispersion can become important near dead-end pipes where velocity is very low. The contribution of an inactivation rate term that follows the Chick-Watson rate law has been included in other models (10, 11) but was excluded in this model for both free (eq 3) and attached bacteria (eq 4). 5926
Methodology
(5)
2
∂Cl2 ∂Cl2 ∂ Cl2 ) -v + Dd 2 - kbCl2 - kw/Rh ∂t ∂x ∂x
Simulations showed that bacterial loss by inactivation was small for the mean value of the reported range of rate constants. On the basis of the values of Cl2,t and Cl2,c that were selected from the literature, the inactivation mechanism would not be important at chlorine concentrations above about 0.5 mg/L because µb is reduced by 90% by the suppression mechanism alone (see eq 1). Moreover, inactivation would not decrease bacterial numbers significantly as chlorine concentration falls below 0.5 mg/L because of the first-order dependency of inactivation rate on chlorine concentration and the fact that chlorine is lost simultaneously through chlorine demand reactions (eq 6). Thus, the addition of a chlorine inactivation mechanism may only cause additional loss of bacteria over a narrow band of residuals in the neighborhood of 0.5 mg/L. In eq 5, Yg is the growth yield coefficient of bacteria (cell uptake of C/substrate of C utilized), and β is the number of bacterial cells that are produced per unit mass of C uptake into cells (109 cells/mg of C uptake in this work), which depends on assumptions of percent C in cells and weight of a cell. In eq 6, kb is the first-order kinetic constant for chlorine decay in bulk water (t-1), and kw is the zero-order rate constant (m l-2 t-1) for the wall reaction (8, 19). Wall demand is not included in the SANCHO model (9) but is included in the PICCOBIO model (10). A series of simulations has been presented elsewhere for a small, illustrative pipe network using the alternating, splitoperator method to solve the above system of equations (12). The baseline values of all input parameters and their units are available in the Supporting Information. The simulations showed that in the absence of variation in water demand, a steady state in free chlorine residual concentration was achieved at all pipe network nodes after only 12 h from start up, where the initial condition was absence of free chlorine, bacteria, and substrate. In sharp contrast, the time to reach a steady-state concentration of substrate and bacteria at each node was much longer (40 d) because the associated microbial processes are much slower than chemical reduction of free chlorine.
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 22, 2004
∂X*b [Xb(R + ∆R) - Xb(R)]/Xb(R) ) ∂R* ∆R/R
(7)
where S* is the normalized, relative sensitivity coefficient, which is also referred to as the elasticity (20); Xb is the state variable of free bacteria; R is any model parameter of interest; and ∆R is the step size of the perturbation (∆R ) 1% × R) in this parameter in a repeated model simulation. The uncertainty in the three parameters (kw, µmax,a, and kb) was carried forward from the sensitivity test to the Monte Carlo simulations. Log-normal probability distributions (22) were developed for each parameter based on the mean (µ) and median values (m) of data from other studies (11, 19, 23, 24) and LOGNORM4.EXE (25) to calculate the standard devia-
TABLE 1. Statistics of the Log-Normal Distribution of Three Model Parameters Selected for Monte Carlo Simulations parameter (h-1)
kb kw (mg m-2 h-1) µmax,a (h-1)
m
µ
σL
0.030 26.5 0.12
0.054 38.8 0.32
0.626 0.873 1.401
FIGURE 1. Characteristics of the hypothetical pipe network to test the regrowth model. Nodes A-D are locations in the model with output and represent water ages at average water demand of 2.1 h (node A), 5.7 h (node B), 9.8 h (node C), and 25.9 h (node D). tion (σL)of the log-transformed values. The log-normal distribution statistics for kb, kw, and µmax,a are listed in Table 1. A pseudo-random number generator was used to sample the resulting log-normal distributions (26); such generators have been shown to work as well as a hardware random generator. The test DS for Monte Carlo simulations of the regrowth model is shown in Figure 1. This simple network consisted of five links and four nodes. A diurnal variation in water demand (from 0.2 to 3.2 times the daily average flow rate) was included with the morning peak being higher than the evening peak. Water quality calculations were made at four pipe nodes (A-D) representing water ages at the average water demand of 2.1 h (A), 5.7 h (B), 9.8 h (C), and 25.9 h (D). The pipe material was assumed to be heavily tuberculated cast iron in order to give a worst-case baseline in which the chlorine decay rate was maximized (i.e., kw for this pipe material is the largest). Each simulation was run for 60 d, a period that was found necessary to approximate steady-state conditions at each node (i.e., to reach less than 5% variation in concentrations of free bacteria, attached bacteria, and substrate). Chlorine residuals reached steady-state far sooner. Temperature was held constant at 20 °C because this would
not be an unreasonable assumption for a 60-d simulation. Each simulation gave output at a 2-h interval over a 24-h time period at four nodes to give 4800 data points; 100 of these simulations were made.
Results and Discussion Sensitivity Testing. S* varied for all model parameters with length along the test pipe (27). This is illustrated in Figure 2 for sensitivity of regrowth to µmax,a. Up to distance of 3.5 km, which corresponds to water age of 8 h, the model is insensitive to µmax,a because the chlorine residual is too high to allow regrowth (eq 1). This emphasizes the important role of chlorine residual in control of regrowth when growth would otherwise be able to occur owing to the presence BDOC. Field measurements have shown, however, that chlorine residual is depleted at long water age, thus allowing regrowth (6). Other related model inputs that govern the distance at which regrowth starts are kw, kb, Cl2,t,a, Cl2,0, and v. S* increases sharply once the chlorine residual declines but then decreases again. The increase is caused by more attached bacteria that subsequently produce more free bacteria by detachment. S* again declines because growth is less influenced by µmax,a once the inhibition effect of chlorine residual is eliminated (eq 1). S* becomes slightly negative near 8 km meaning that an increase in µmax,a produces less growth as a result of rapid depletion of the substrate concentration (eq 1). The same explanation holds for switches in the direction of sensitivity for µmax,b, kdep, Yg, Ks, kdep, and kdet. The baseline values of each model parameter are listed along with the maximum S* observed along the pipe length after 60 d of simulation in Table 2. The most sensitive parameter was kw followed by µmax,a. Both imply that surface reactions are more important than bulk reactions, consistent with other reported simulations (12). The maximum S* was somewhat higher for Cl2,t,a than for kb, but the latter was nonetheless selected for subsequent Monte Carlo simulations because kb has been widely linked to regrowth events, and many measurements were available to generate a credible log-normal distribution of values. Future research, however, should address confirmation of eq 1 and Cl2,t,a in particular. Effect of Diurnal Variations in Water Demand. The daily variations in concentrations of chlorine and free bacteria at nodes A-D for three successive days (days 58-60) are shown in Figures 3 and 4, respectively. The cyclical variation in chlorine at stations A-C is caused by the diurnal variation in water demand (a factor 16) that directly affects the time available for chlorine decay reactions to occur. Thus, the
FIGURE 2. Sensitivity of free bacteria concentration to maximum specific growth rate constant of attached bacteria (µmax) with length along 8-in. diameter pipeline. Relative sensitivity coefficient is calculated from eq 7. VOL. 38, NO. 22, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
5927
TABLE 2. Maximum Relative Sensitivity Coefficients of the Parameters and Variables in the Regrowth Model model parameter
kw µmax,a Cl2,t,a kb Cl2,0 v S0 kd Ks kdet Dd Cl2,c kdep Xb,0 µmax,b Yg Cl2,t,b
baseline value m-2 h-1
26 mg 0.16 h-1 0.20 mg/L 0.03 h-1 1.0 mg/L 0.13 m/s 0.40 mg/ L 0.06 h-1 0.4 mg of C/L 0.03 h-1(m/s)-1 0.005 m2/s 0.20 mg/L 0.20 h-1 10 cell/mL 0.20 h-1 0.15 mg/mg 0.05 mg/L
maximum S* 420 (+) 170 (+) 120 (+) 95 (+) 80 (-) 80 (-) 65 (+) 53 (-) 40 (-) 13 (-) 7 (-) 6 (+) 1.5 (+) 1 (+)
