Environmental Modelling & Software 51 (2014) 12e25

Contents lists available at ScienceDirect

Environmental Modelling & Software

journal homepage: www.elsevier.com/locate/envsoft

Drinking water distribution systems contamination management to

reduce public health impacts and system service interruptions
Amin Rasekh*, Kelly Brumbelow
Zachry Department of Civil Engineering, Texas A&M University, 3136 TAMU, College Station, TX 77843-3136, USA

a r t i c l e i n f o a b s t r a c t

Article history: Decisions on protecting public health against drinking water systems contamination threats should be
Received 2 September 2012 made with careful consideration of credibility of threat observations and adverse impacts of response on
Received in revised form system serviceability. Decision support models are developed in this study to prepare water utility op-
23 September 2013
erators for making these critical decisions during the intense course of an emergency. A pressure-
Accepted 24 September 2013
dependent demand model is developed to simulate the system hydraulics and contaminant propaga-
Available online 17 October 2013
tion under pressure-deficit conditions that emerge after the response actions are executed. Contrary to
conventional demand-driven models, this hydraulic analysis approach prevents potential occurrence of
Keywords:
Water distribution system
negative pressures during the simulation and may identify better response protocols through exploring a
Contamination larger search space. Response mechanisms of contaminant containment and discharge are optimized
Emergency management using evolutionary algorithms to achieve public health protection with minimum service interruption.
Evolutionary computation Sensitivity analyses are conducted to assess optimal response performance for varying response delay,
Multiobjective number of hydrants, and intrusion characteristics. Different methods for quantifying impacts on public
Pressure-dependant demand analysis health and system serviceability are explored and the sensitivity of the optimal response plan to these
different formulations is investigated. The simulation-optimization schemes are demonstrated and
discussed using a virtual water distribution system.
Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction risk of health impacts has returned to pre-event levels). An emer-

gency response protocol explains actions that managers may take
Aging drinking water infrastructure and increased risks of in response to the perceived state of the system after the emer-
terrorism have intensified concerns for vulnerability of water dis- gency begins, and it considers how best to achieve managers’
tribution systems (WDS) to accidental and intentional contamina- multiple objectives. These response actions can be classified as
tion. A contamination event may cause health and sociopolitical “assessment,” “preventive,” or “protective” actions, depending on
impacts, erode public trust, and interrupt system operation. To whether they collect information about the state of the system,
effectively cope with these threats, there is a need to prepare operate on the system to decrease impacts, or require action by the
contamination emergency management plans that describe the public to reduce exposure, respectively (Perry and Lindell, 2007).
actions a drinking water utility should take in preparation for and in Title IV of the Public Health Security and Bioterrorism Prepared-
response to a contamination threat or incident. An emergency ness and Response Act of 2002 (United States Congress, 2002) re-
management plan should be based upon careful risk assessments quires all community water systems serving a population greater
and cover the four phases of hazard mitigation, emergency pre- than 3300 in the United States to prepare or revise emergency
paredness, emergency response, and disaster recovery (Lindell response plans. The Response Protocol Toolbox (RPT) has been pre-
et al., 2006). pared by the United States Environmental Protection Agency (USEPA,
A contamination emergency response phase is initiated with an 2003) to help water utilities meet this requirement. It provides
actual (or potential) release of contaminant that spreads across a general guidelines on how response decisions should be made at the
WDS, and it extends until the situation is stabilized (i.e., when the various stages of a contamination event as more information is
collected. Because this toolbox is essentially a qualitative document,
however, it does not provide specific guidance on how appropriate
* Corresponding author. Tel.: þ1 979 739 9652; fax: þ1 979 862 1542.
response strategies should be devised for a particular WDS. This
E-mail addresses: [email protected] (A. Rasekh), kbrumbelow@ study develops quantitative simulation-optimization models to
civil.tamu.edu (K. Brumbelow). prepare emergency response protocols that specify functional

1364-8152/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.envsoft.2013.09.019
 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 13

contaminant containment and flushing operation rules for the contamination scenario which is treated as a model input here are
achievement of conflicting response objectives. explained and a series of existing contaminant source identification
methods are introduced. This is followed by a detailed description
2. Problem statement and solution approach of the pressure-dependant demand and exposure models for the
simulation of contamination events under pressure-deficit condi-
Emergency response is a progressive, interactive, and adaptive tions. Optimization objective functions and decision variables are
process that includes parallel activities of assessing unusual then explained along with an overview of the evolutionary-
contamination observations and making appropriate emergency computation-based optimization algorithms used in this study.
response decisions. As more information is obtained about
contamination, emergency management progresses through three 3.1. Contamination scenario
threat stages of “possible,” “credible,” and “confirmatory” (as
described in RPT) accompanied by an increase in seriousness of the Optimization of emergency response protocols needs to be
threat impacts and magnitude of response decisions. While public performed for a given WDS contamination scenario, which is
health protection is the primary response focus, emergency man- defined by a set of attributes including: (1) site(s) of contaminant
agement should carefully consider other potential consequences on intrusion, (2) contaminant type, (3) contaminant mass, (4) time of
infrastructure serviceability due to response implementation, spe- year, (5) the time of day the contamination event is started, and (6)
cifically in the early stages of the process where the attack credibility the intrusion duration (Rasekh and Brumbelow, 2013). The
level is relatively low. At this stage, a multiobjective response plan is simulation-optimization models developed here treat the
useful to identify the balance between actions taken to protect contamination scenario as input information. This scenario could
public health against a potential threat and limiting overaction that be a potential critical and base scenario for which the response
adversely impacts the ability of the system to meet multiple aspects protocols need to be optimized before an emergency occurs
of its overall mission. Nevertheless, if streaming threat information (Perelman and Ostfeld, 2010; Rasekh and Brumbelow, 2013).
and observations corroborate occurrence of a contamination event, Alternatively, it can be a scenario that is occurring and being
minimizing potential health impacts becomes the sole primary characterized by applying a contaminant source identification
objective that should be sought for. Under these circumstances, the model (Preis and Ostfeld, 2006; Zechman and Ranjithan, 2009;
large size of multiobjective optimization results may not be deci- Hart et al., 2009; Liu et al., 2011a; Rasekh and Brumbelow, 2012;
pherable for making timely emergency decisions and use of a single- Gugat, 2012; Shen and McBean, 2013) linked with the sensor
objective model may become preferable. network (Ostfeld et al., 2008; Janke et al., 2009). Therefore, the
To date, limited research has addressed this multicriteria nature simulation-optimization models developed here can be employed
of the contamination emergency response problem (Preis and after the potential critical or design basis scenarios have been
Ostfeld, 2008; Alfonso et al., 2010). Multiobjective frameworks characterized during the emergency preparedness phase or a real
proposed so far have only considered hydrant and valve locations as contamination scenario has been identified during a contamination
decision variables and have not optimized the operation timing. emergency.
While these studies have considered the number of operational
actions as an emergency response criterion, the important criterion 3.2. Hydraulic simulation under pressure-deficit conditions
of system service interruption has not been explicitly addressed.
Moreover, previous single and multiobjective studies have used Behavior of a WDS under normal operating conditions is most
demand-driven analysis (DDA) to simulate WDS behavior, and this commonly simulated using standard DDA models like EPANET
assumption inevitably limits the optimization search space to (Rossman, 2000) for design, operation, and rehabilitation purposes.
response protocols that do not cause excessively low pressure in DDA models are formulated on the premise that nodal water de-
the WDS. This may unfavorably filter out many possible response mands are known and completely met during the simulation
protocols with a high potential for reducing the health impacts. period so that nodal pressure and pipe flows can be calculated by
In the light of these needs, this study develops and integrates a solving a system of quasi-linear equations (Wu et al., 2006). Pro-
hydraulic pressure-dependant demand model (PDDM) and evolu- jections of network behavior that are based on DDA are reasonably
tionary optimization schemes to find the optimal emergency accurate under normal conditions when pressures are sufficiently
response protocols with explicit consideration of two important high. Under abnormal pressure-deficit conditions caused by
response criteria: (1) public health impacts, and (2) system service emergency response actions, however, DDA may illustrate a dis-
interruption. Emergency response is treated as both single and torted image of the true system behavior.
multiobjective optimization problems to address utility managers’ A standard DDA solves a system of energy and continuity
needs under different conditions. Operational rules for contami- equations to calculate unknown nodal heads and pipe flow rates
nant containment and discharge locations and timing are explicitly (Todini and Pilati, 1988; Rossman, 2000). This system of equations
treated as optimization decision variables. Sensitivity analyses are may be extended to include pressure-dependant demand functions
performed to provide insight into effective response protocols and (PDDF) as well to relate pressurized water availability to existing
assess sensitivity of optimal protocols to different parameters such nodal head under pressure-deficit conditions (Laucelli et al., 2012).
as response delay. Different formulations for quantifying impacts Full pressure-driven analysis approaches solve this extended sys-
on public health and service availability are examined with the help tem of equations simultaneously to determine unknown nodal
of the PDDM and an exposure model. Performance of the proposed demands, nodal heads, and pipe flow rates (Laucelli et al., 2012).
schemes is investigated using the WDS of Mesopolis virtual city, Alternatively, the standard DDA model may be iteratively run and
which resembles the interdependency and interconnectedness of nodal demands are updated sequentially after each iteration using
realeworld complex water distribution networks. the PDDFs until a satisfactory convergence is achieved (Liu et al.,
2011b; Jun and Guoping, 2013; Kanta and Brumbelow, in press).
3. Model development This modeling approach is used in this paper. A major advantage of
this approach is that it allows benefiting from the computational
Contamination scenario and the simulation and optimization efficiency and robustness of hydraulic and quality simulators of the
models will be described in this section. Different attributes of a well-established EPANET software. The fact that many existing
14 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25

Fig. 1. Flowcharts of simulation and optimization models.

WDS models are developed in this modeling environment and

8
most water utilities are familiar with EPANET further highlights this > 0 Hi < Himin
>
>
advantage. >
<  1=2
Several PDDFs have been proposed in previous studies that may Qi ¼ nom Hi Himin (1)
> Qi Hiser Himin
Himin < Hi < Hiser
be used for a pressure-driven hydraulic simulation (Shirzad et al. in >
>
>
:
press). This study uses the PDDF proposed by Wagner et al. (1988), Qinom Hi > Hiser
which is among the most popular formulations (Laucelli et al.,
2012) and is the best suited for describing water demand sup- A general framework of the simulator is shown in Fig. 1. Every
plied to customer (Giustolisi and Walski, 2011). Shirzad et al. (in demand node in the WDS belongs to one of the three categories
press) recently conducted a comprehensive laboratory study to defined by Eq. (1) at every iteration and its category may switch
compare different pressure-discharge relations with application in depending on their nodal head at the next iteration. Interested
pressure-driven WDS simulation. They concluded that Wagner’s readers may refer to Appendix A of Kanta and Brumbelow (in press)
PDDF best matched experimental data since it closely resembles for a detailed description of the iterative PDDM used here and step-
the orifice relation. The Wagner’s model is based on the rationale by-step implementation instructions.
that hydraulic laws for flows through devices indicate that the flow The results of PDDM are saved as an output hydraulic file,
is proportional to the square root of head. It considers a parabolic which is then used by the EPANET quality analysis model to
relationship between the service head Hiser for node i and the simulate the propagation of contaminant throughout the WDS.
minimum head Himin needed to satisfy the nominal demand Qinom : The contaminant transport is simulated as a perfect tracer: density
 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 15

effects, decay, and reaction with pipe wall materials and other X
Np
existing species are not considered. EPANET-MSX (Multi-Species fH2 ¼ mi (3)
eXtension) (Shang et al., 2008) may be used to enable EPANET to i¼1
model complex reaction schemes between multiple chemical and
biological species. The obtained WDS quality analysis results, where Np ¼ the population size served by the WDS; ai ¼ a binary
represented as spatial and temporal concentration of contaminant index; md ¼ known toxic dose; mi ¼ total contaminant mass
in the system, are exported to a human exposure model for eval- ingested by individual i during the whole course of event that is
uation of public health impacts. calculated by the PDDM and exposure model:

X
NI
3.3. Human exposure model mi ¼ Vi;j  Ci;j (4)
j¼1
Adverse health impact, defined here as either the number of
sicknesses or the total ingested contaminant mass resulting from a where NI ¼ number of water ingestion events for consumer i,
contamination event, is estimated using the PDDM and water Vi,j ¼ volume of water ingested by individual i at ingestion event j,
quality simulation model coupled with an exposure model. The and Ci,j ¼ concentration of contaminant in water volume ingested
quantity of contaminant ingested by individuals during a by consumer i at ingestion event j. Vi,j ¼ 0.186 liter (calculated
contamination event depends on water ingestion pattern, time- through dividing daily water ingestion rate of 0.93 by 5 daily in-
varying concentration of contaminant, and availability of drink- gestions times) if water is available for individual i at ingestion
ing water under pressure-deficit conditions. An individual is event j, and zero, otherwise. Water availability is affected by the
assumed to become sick if the cumulative amount of contaminant flushing and isolation actions taken by the managers during the
ingested during a contamination event exceeds a known toxic emergency.
dose. The timing-of-ingestion model selected for this study as- Function fH1 explicitly represents the health impacts in terms of
sumes that tap water is ingested at the common starting times for morbidity while fH2 projects the health consequences in the more
the three major meals (7:00, 12:00, and 18:00) and times halfway implicit form of ingested mass. However, fH2 does not require that
between these meals (9:30 and 15:00), given that there is water toxic dose be known in advance, which greatly simplifies the
available for drinking at a consumer’s demand node. The tap water analysis since this value is often difficult to determine and subject
intake rate used in the exposure model is central to accuracy of to high variance among individuals. When information is available
estimated impacts and is set to 0.93 L/day here based on USEPA to construct a probabilistic doseeresponse model (Chick et al.,
(2004). 2001), the population fH1 affected by a particular scenario may be
estimated probabilistically with consideration of the dosee
response uncertainties (Ostfeld et al., 2008).
3.4. Optimization objective functions
3.4.2. WDS service interruption
Quantification of emergency response criteria is the first step for
Interruption to system serviceability is also mathematically
the preparation of the optimization model. Different measures may
expressed using two different measures. The first measure is the
be formulated to quantify two important objectives of reducing
total number of hours the volume of water supplied to consumers is
public health impacts and system service interruptions. This study
below a certain percentage of their demand. Alternatively, the
investigates two distinct formulations for each of these criteria and
second measure sums up the difference between water demand
analyses the sensitivity of optimal response plan to each quantifi-
and supply for all consumers during the whole course of the
cation method.
emergency. Mathematically,

3.4.1. Public health impacts X

tend X
Nc
Some studies have quantified the health impacts in terms of fS1 ¼ bti ; bti ¼ 1 if ðSti < lDti Þ; bti ¼ 0 otherwise
contaminant concentration in system nodes either as total t¼0 i¼1
contaminant concentration (Baranowski and LeBoeuf, 2008) or (5)
total number of nodes with concentration above a specified
threshold (Alfonso et al., 2010). This approach does not account for
X
tend X
Nc
the ingestion timing and rate which can significantly influence fS2 ¼ ðDti  Sti Þ (6)
estimated impacts (Davis and Janke, 2008) and does not consider t¼0 i¼1
the critical fact that the connections serving higher populations are
comparatively more vulnerable. In general, health impacts metrics where tend ¼ duration of emergency response phase in hours;
may be differentiated depending upon whether or not they are Nc ¼ total number of consumers; bti ¼ a binary index;
based on a threshold (i.e., number of people/nodes that experience Sti ¼ volume of water supplied to consumer i at time step t
a concentration or contaminant ingested mass above a threshold) calculated by the PDDM model; Dti ¼ water demand of consumer
(Perelman and Ostfeld, 2010; Zechman, 2011; Zechman et al., 2011) i at time step t; l ¼ supply deficit threshold that ranges from 0 to
or sum of concentration or consumed mass for all nodes/people 100%. The function fS1 treats all consumers equally while fS2
(Baranowski and LeBoeuf, 2008; Preis and Ostfeld, 2008). Accord- places more emphasis on consumers with higher water demand.
ingly, using the developed exposure model, the health impact cri- Both these metrics are always unrealistically projected as zero if
terion is quantified here as either the total number of sicknesses fH1 DDA is used.
or the total contaminant mass fH2 ingested by all people during the
whole course of the event: 3.5. Optimization decision variables

X
Np Fig. 2 shows a schematic arbitrary timeline of contamination
fH1 ¼ ai ; ai ¼ 1 if mi > md ; ai ¼ 0 otherwise (2) emergency period and contaminant containment and discharge
i¼1 operations. The optimal order of response actions and their timing
16 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25

Response Flushing be performed to identify the most reasonable configuration of these

delay duration CVSs. In this study, it is based upon WDS layout, extensive hydraulic
CVS closure
simulations, WDS designers’ understanding of the system behavior,
duration
and computation intensity considerations. Future research may
enhance and automate the configuration process using systematic
Time optimization (Giustolisi and Savic, 2010), network clustering
(Perelman and Ostfeld, 2011), and minimum-cut set approaches
Hydrant CVS closure (Yannopoulos and Spiliotis, 2013), among others.
open time time Contaminant flushing is the other preventive response mecha-
Intrusion duration nism that is executed through opening hydrants to discharge a
large volume of contaminated water. System flushing should be
Fig. 2. A schematic timeline of a contamination emergency period. The variables that
are optimization decision variables (i.e., elements of a response protocol) are shown in
planned and implemented carefully so that it is performed at the
italics. sections of the system where contaminant concentration is higher.
Otherwise, it will worsen the situation by further spread of
contaminated water to uncontaminated areas as it can considerably
are identified by the single and multiobjective optimization alter flow regime (USEPA, 2003; Haxton et al., 2012). The decisions
models, which will be described in detail later in Section 3.6. include selecting the subset hydrants that should be opened (h¼
Contaminant containment through isolation valve operations is {h1,h2,.,hnh}) among the complete set of hydrants to flush
a preventive response mechanism implemented to prevent contaminated water, time when hydrants should be opened after
contaminant spread to uncontaminated regions of a network and to the response delay (th¼{th,1,th,2,...,th,nh}), and time period when
preclude consumers from withdrawing contaminated water hydrants should be remained opened (Dth¼{Dth,1,Dth,2,.,Dth,nh}). A
(USEPA, 2003). Decisions on combination of valves to be closed and maximum number of hydrants,nh, may be opened during the
timing must be made carefully and implemented quickly to be simulation depending on personnel and equipment availability.
effective and minimize accompanying impacts on non- Sensitivity analyses are performed on nh to assess how this
consumptive uses in isolated regions. The decisions include the parameter affects the response performance.
time at which closure valve sets (CVS) should be closed after the Implementation of response actions should account for the
response delay (ts¼{ts,1, ts,2,., ts,ns}) and the duration that the CVSs response delay times between when intrusion is initiated and
should remain closed (Dts¼{Dts,1,Dts,2,.,Dts,ns}), where ns is the when injection of contaminant in the network is identified as likely
total number of CVSs. The optimization model selects a subset of or confirmed through multiple contamination trigger events. It is
the whole set of ns CVSs to be operated. When it returns Dts ¼ 0 for a also during this time period that the contamination scenario is
particular CVS, it means that CVS should not be operated. identified using a source characterization model as explained in
The set of CVSs should be configured before the optimization is Section 3.1. Sensitivity analyses are conducted on response delay in
performed. While every single pipeline may be theoretically consid- this study to understand how delays influence the effectiveness of
ered as a possible closure location, this would result in a tremen- response actions.
dously large decision space that may not be handled by the optimizer
practically. A large set would also lead to elaborate response protocols 3.6. Optimization algorithms
that may be burdensome to execute in practice. CVSs that include
fewer valves will minimize the number of operational actions needed The objective functions described in Eqs. (2)e(6) are minimized
to isolate an area and are thus more preferable. Assessments need to through optimization of preventive response mechanisms of

Fig. 3. Configuration of closure valve sets (CVS) in Mesopolis.

 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 17

Fig. 4. Water consumption patterns for three representative types of consumers in

Mesopolis. Fig. 5. Accumulation of total ingested contaminant mass over time. The ultimate value
at the end of the set emergency period (i.e., 6 days) denotes fH2 (see Equation (3) and
Table 1).

contaminant containment and flushing. Emergency decisions that

should be optimized in response to a contamination threat are
et al., 2002), an elitist evolutionary algorithm that benefits from a
diverse and discrete, and the relationship between decisions and
fast non-dominated sorting strategy and does not need any user-
performance is complex. Evolutionary algorithms have been
defined parameter for preserving diversity in Pareto optimal
demonstrated as flexible and powerful tools for solving such
surface.
complicated optimization problems in the field of water resources
Performance of NSGA-II is evaluated using the hypervolume
management (Reed et al., 2013; Afshar et al., 2009; Rasekh et al.,
measure (Zitzler and Thiele, 1998). Hypervolume is the non-
2010; Savic et al., 2011; Laucelli and Giustolisi, 2011; Takbiri and
overlapping volume in objective space covered by members of a
Afshar, 2012; Labadie et al., 2012). A genetic algorithm (GA) is an
non-dominated ensemble of solutions with respect to a reference
adaptive search heuristic that mimics the process of genetic and
point, set as the worst value for each objective. The hypervolume is
natural evolutions (Holland, 1975). It uses strings of decisions var-
applied to evaluate the convergence of the non-dominated front as
iables (also referred to as chromosomes) to represent a possible
it evolves. Increasing values of the hypervolume necessarily indi-
solution to the problem. A GA starts with a population of chro-
cate an improvement in the uniformity and proximity to a true
mosomes and performs recombination and random perturbation to
Pareto optimality, which are characteristics that are indicative of
successively generate better solutions. Solutions with higher fitness
the quality of a non-dominated set of solutions.
values have a higher chance of survival and being transmitted to the
next generation. Fitness progressively improves through successive
generations and the algorithm converges to an optimal solution. 4. Application example
Single and multiobjective genetic algorithms are developed here
and coupled with the PDDM and exposure model for identification Mesopolis virtual city (Johnston and Brumbelow, 2008) is used
of emergency response protocols (Fig. 1). here to demonstrate and discuss performance of the proposed
The single-objective model uses an elitist real-coded genetic simulation-optimization schemes for the identification of optimal
algorithm with roulette wheel selection, simulated binary cross- response protocols. In what follows, a general description of the
over (SBX) (Deb and Agrawal, 1995), and polynomial mutation Mesopolis WDS, two considered contamination scenarios, and
(Deb, 2001). In the roulette wheel selection, the probability that a model settings are provided. The single-objective optimization
solution will be selected is given by the ratio of its fitness to the model is then applied to find the minimum-health-impacts hy-
total fitness of other members of the current population. The drant operation rules and sensitivity analyses are performed on the
offspring hydrants are generated through performing crossover and number of hydrants and response delay. The multiobjective opti-
mutation operations on longitude and latitude coordinates of the mization is finally performed to optimize simultaneous contami-
parent hydrants and picking the hydrant nearest to the resulting nant containment and discharge practices for the identification of
coordinates. the non-dominated response protocols with varying trade-offs
In contrast to the single-objective optimization approach that between the response objectives of public health protection and
searches for a single solution with the best scalar fitness value, service interruption prevention.
multiobjective optimization seeks a set of trade-off solutions which
together define the best multiobjective alternatives surface called 4.1. WDS, scenarios, and settings
the Pareto-optimal front (Deb, 2001). Since GAs are population-
based algorithms, they offer a means of finding the complete The WDS of Mesopolis is comprised of one reservoir, two water
Pareto-optimal front in a single optimization run. This study em- treatment plants (WTP), 13 tanks, 65 pumps, 876 hydrants, and
ploys Non-dominated Sorting Genetic Algorithm II (NSGA-II) (Deb 2062 water mains (Fig. 3). Demands are exerted at 706 residential,

Table 1
Contamination scenarios.

Scenario Location Loading (kilograms) Demand multiplier Start time Duration (hour) fH1 (people) fH2 (grams)

West West WTP 300 1.00 18:00 6 33,944 304.08

East East WTP 300 1.00 19:00 5 54,638 397.28
18 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25

eastern area through closing only two mains. Closure of CVS 5 and 6
shuts down the East and West WTPs, respectively. CVS 7 isolates
the highly populated residential area in the western part of the city
from both WTPs by closing only two mains. CVS 8 allows the
western area to be supplied by only West WTP and is thus helpful if
the East WTP is contaminated. Finally, CVS 9 includes three mains
and completely disconnects the eastern network from central and
western areas.
Model emitter discharge coefficients for hydrants are set to
166.5 gpm/psi0.5 (associated with a 3-inch diameter connection fire
hydrant) to calculate pressure-dependant outflow in the PDDM.
Values of Hiser and Himin for the PDDM are based on engineering
design standard of the cities of Bryan and College Station, Texas
(Cities of Bryan and College Station, 2005). Under normal condi-
tions, a design head of 35 psi should be maintained throughout the
Fig. 6. Single-objective optimization convergence history for five GA runs. It is asso- system. Water outflow at a connection is assumed to be unavailable
ciated with the East Scenario, a 6-h response delay, and 5 hydrants.
if its pressure drops below 20 psi, the minimum allowed during fire
flow conditions. Because the PDDM developed here involves an
iterative process, a convergence criterion should be devised. It is
industrial, and commercial/institutional nodes, representing a mid- considered to be when the percentage of demand nodes whose
size city of nearly 147,000 residents. The diurnal demand patterns current category has been switched falls below a preset value. It is
for three representative consumer types are illustrated in Fig. 4. set to 1% in this study. Under the very rare occasions when the
This study accounts for the spatial and temporal variability of de- iteration does not converge after a preset maximum time period,
mands but consideration of demand uncertainties requires appli- which is two minutes here, the solution is discarded. The optimi-
cation of optimization algorithms beyond those used here (Jin and zation results show that it occurs to only nearly 0.2% of all solutions.
Branke, 2005; Kasprzyk et al., 2012).
Response plans are optimized here for two particular contami- 4.2. Single-objective optimization
nation scenarios that are listed in Table 1. The demand multiplier
associated with each scenario is representative of aggregate water The single-objective model is applied first to identify optimal
demand for a WDS that typically varies throughout the year. The subset of hydrants and their operation timing that minimize the
contaminant agent is arsenic with a toxic dose of md ¼ 3.5 mg for a health impacts formulated as fH1 and fH2. Appropriate genetic al-
body weight of 70 kg as reported by White (1999). The contami- gorithm settings are determined based on sensitivity analyses.
nation zone associated with West and East Scenarios is their major Crossover and mutation rates are 0.80 and 0.05, respectively. SBX
distinction. Contamination of West WTP leads to contamination of crossover distribution and polynomial mutation indices are set to
a limited western region, while the entire Mesopolis network is 15 and 5, respectively. Population size is set to 40, and the model is
contaminated if the contamination occurs at the East WTP. Both run for 200 generations. Fig. 6 shows optimization convergence
contamination scenarios occur in the third day of simulation, which during the evolution process of five random GA runs for the East
is when the system has reached dynamic equilibrium. Temporal Scenario, nh ¼ 5, a response delay of 6 h. The observation that all the
accumulation of total ingested contaminant mass as the emergency randomly seeded runs lead to a similar reduction in health impacts
proceeds is illustrated in Fig. 5 for both contamination scenarios. It may be experimentally interpreted as the model’s acceptable
is observed that the health impacts practically reach a maximum robustness and efficiency in finding optimal response protocols.
intensity after 6 days and this duration is thus set as the model Complete demonstration of the model’s performance, however, is
simulation period. subject to availability of the global optimal response protocols.
Nine CVSs are considered here as shown in Fig. 3 based upon Optimal reduction in health impacts with increasing number of
WDS connectivity, hydraulic simulations, WDS developers’ under- hydrants is illustrated in Fig. 7 for both contamination scenarios
standing of the system behavior, and computation intensity limi- and health impact formulations. Increasing number of flushing
tations. CVS 1, 2, and 3 are located on long mains that transmit locations leads into greater reduction in health consequences,
water to far sections of WDS. CVS 4 allows isolation of a large which is not surprising. The consistently decreasing slope of the

Fig. 7. Percentage reduction in health impacts for varying number of hydrants for a 6-h response delay using two formulations of (a) total number of sicknesses fH1, and (b) total
ingested mass of contaminant fH2.
 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 19

Fig. 8. Percentage reduction in health impacts for varying response delay for 5 hydrants using two formulations of (a) total number of sicknesses fH1, and (b) total ingested mass of
contaminant fH2.

curves, however, suggests that opening an additional hydrant is not the population sectors that are at high exposure risk (i.e.,
as effective when a larger number of hydrants are opened. ingested contaminant mass above the toxic dose). Fig. 9(b) shows
Degradation of optimal response protocols performance with that the protocols that are based on fH1 may even magnify the
increasing response delay is observed in Fig. 8. Prolonged public exposure severity for the population who are at a very high risk
exposure and extended contamination area are majorly responsible (ingested contaminant mass above nearly 10 mg). These obser-
for this decreasing efficiency of emergency response practices. The vations may be explained by the fact that every person in the
response mechanism of contaminant flushing is reasonably effec- population that ingest a mass of contaminant greater than the
tive for short response delays but its performance rapidly di- toxic dose contribute equally to fH1 measure, no matter how
minishes with increasing delay. Figs. 7 and 8 indicate that much greater. This is, however, not true for the fH2 measure.
contamination would be noticeably better controlled when it Considering the population at low risk (ingested mass below the
happens at the West WTP than the East WTP. This observation may toxic dose), it is similarly observed that the response protocols
be attributed to the smaller impact area associated with the West based on fH2 are superior when the East Scenario occurs. When
Scenario. the West WTP is contaminated, however, this observation does
Public exposure to contaminated water varies for different not hold.
population segments across the city once a contamination event The optimal response protocols associated with the two
occurs. The understanding of this variance and how it is affected by contamination scenarios, fH2 health impacts formulation, and
implementing response actions can be improved through analyzing response delays of 6 and 12 h are graphically illustrated in Fig. 10. It
the cumulative distribution curves of ingested contaminant mass is observed that the optimal hydrant locations are concentrated in a
presented in Fig. 9. The vertical axis indicates the percentage of limited western area when the West WTP is contaminated whereas
total population that has ingested a mass of contaminant below the the optimal flushing locations are scattered over a larger eastern
corresponding value on the horizontal axis. fH1 is the population on region if the East WTP is the contamination site. This is linked to the
the vertical axis above the cross point of toxic dose line with each impact areas associated with the West and East contamination
cumulative curve, while fH2 is the area between the curve and scenarios. Provision of these optimal response maps for finer dis-
vertical axis after it is transformed to absolute population values cretization of response delay and other vulnerable nodes would
(i.e., multiplying the percentage values by the total population of serve as a supportive visual aid for utility managers to respond to
Mesopolis). the contamination in a more timely and efficient manner.
Fig. 9 shows the response distribution curves that are based Fig. 10 also shows the optimal timing associated with the
on the health impact formulation of fH2 are above those based on flushing locations for different scenarios and response delays. It
fH1 when ingested contaminant mass is larger than the toxic dose is clearly observed that the optimal flushing locations move
of 3.5 mg. This observation indicates that the optimal response farther from the contamination site as the response delay in-
protocols that are obtained through minimizing fH1 better protect creases and the contaminant plume moves further throughout

Fig. 9. Distribution of health impacts for minimum-health-impact protocols for (a) West Scenario and (b) East Scenario. Number of hydrants and response delay are 5 and 6 h,
respectively, for both scenarios.
20 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25

Fig. 10. The optimal discharge locations and timing associated with minimum-health-impact protocols using fH2 formulation.

the system. Moreover, the hydrants would be better opened 4.3. Multiobjective optimization
immediately after the response delay to minimize the health
impacts, which is not surprising. The optimal values of flushing Response mechanisms of contaminant containment and
durations, however, might not be as predictable; the optimal discharge are next optimized simultaneously to minimize the im-
durations vary significantly over different discharge locations. An pacts on public health and system service. The multiobjective
analysis inspired by this observation shows that setting the du- optimization algorithm coupled with the PDDM identifies the most
rations to maximum possible values leads to degradation of efficient subset of hydrants and CVSs and their optimal operation
obtained optimal response protocols performance and this is timing for minimizing the health impacts and service interruptions.
consistent for all protocols. This rather unexpected observation The response protocols are optimized for the East Scenario, a
might be linked to the fact that while hydrant operation can response delay of 6 h, and nh ¼ 5. The total number of decision
discharge contaminated water, it also changes contaminant variables is thus 33 (15 for contaminant discharge and 18 for CVSs
plume and spread across the system. A prolonged flushing may operation). NSGA-II algorithm population size is 100, and the
not thus necessarily lead to an increased performance as it might optimization evolution process is continued for 200 generations.
divert the contaminant plume to originally safe regions across Crossover and mutation rates are set to 0.80 and 0.08 and SBX
the systems. crossover distribution and polynomial mutation indices are 15 and
 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 21

Fig. 11. Multiobjective optimization convergence history: (a) hypervolume evolution for three NSGA-II runs, (b) non-dominated front evolution corresponding to the hypervolume
curve plotted as red solid dash type. It is associated with the East Scenario, a 6-h response delay, and 5 hydrants.

10, respectively. The tournament size for the NSGA-II selection Eqs. (2)e(6), are illustrated in Fig. 12. Illustrated Pareto fronts
operator is set to 3. Additionally, the reference point for calculating demonstrate that there are significant trade-offs between the two
the hypervolume convergence measure is fH1 ¼ 54,638 people and response criteria of public health protection and system service-
fS1 ¼ 10,000 block-hours. ability. In light of the fact that impacts on public health are
Convergence of the hypervolume measure during the evolution considered more crucial, one may pose the question whether we
process for three randomly-seeded NSGA-II runs is shown in should be concerned about the trade-offs at all. In the trade-off
Fig. 11(a). This is associated with the East Scenario, response curves found here, each individual response protocol is associated
objective formulations of fH1 and fS1, five hydrants, and a 6-h with a level of reduction in health impacts and accompanied sys-
response delay. Increasing value of hypervolume as the number of tem interruption. A central point to remember is that, while
generations increase is an indicative of convergence to better non- implementation of such a protocol would lead to system inter-
dominated sets of response protocols as illustrated in Fig. 11(b). ruption for certain, reduction in health impacts is conditioned on
Obtained frontier migrates away from the worst point towards the the credibility of threat observations (i.e., the RPT sequencing of
coordinates reference point and spreads uniformly across both possible, credible, and confirmatory threat stages). For the trade-off
axes. The search practically converges and stabilizes after approx- analysis to be rational, the horizontal axis should be multiplied by
imately 100 generations. the probability that the contamination has actually occurred in
The obtained Pareto-optimal fronts for four possible combina- order for it to be consistent with the vertical axis representing
tions of objective functions formulations, which are represented by unconditioned system interruption. In reality, however, such

Fig. 12. Pareto-optimal fronts for minimizing public health impacts (fH1 and fH2) and system service interruption (fS1 and fS2) for the East Scenario, a 6-h response delay, and 5
hydrants.
22 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25

Fig. 13. The minimum-health-impact protocols using fH1 and fH2 formulations for East Scenario and a 6-h response delay.

probability values are not known exactly and must be inferred from behavior in the decision space and the effect of health impacts
an ensemble of uncertain sensor triggers and unusual observations. formulation used. Both protocols are highlighted in Fig. 12 and are
Thus, the trade-off curves must be understood through this filter of associated with the fS2 service interruption formulation and. Fig. 13
threat uncertainty. illustrates the optimal decision variables values. Spatial spread and
The maximum reduction in health impacts when the hydrants cumulative distribution curves of the health impacts associated
are operated solely is marked in Fig. 12. This helps to understand with each of the two response protocols are also shown in Figs. 14
how much the response performance is improved when the CVSs and 15, respectively.
are operated as well. It is observed that this leads to a significant Optimal hydrants and CVSs operation instructions are depicted
improvement in public health protection. For the fH1efS1 combi- in Fig. 13. The hydrants should be opened immediately after the
nation, for instance, the percentage reduction in health impacts is response delay and flushing duration varies for different hydrants,
almost doubled (from 40% to 79%) when the CVSs are operated which is in agreement with the single-objective optimization
together with the hydrants. findings. Immediate closure of CVS4 for both protocols blocks
Fig. 12 illustrates the non-dominated response protocols in the further spread of contaminant to the eastern regions of the city.
objective space. Further analysis has been performed on two Closure of CVS6 is the other common action for the two response
selected minimum-health-impact protocols to understand their protocols. This creates a pressure gradient, which may guide the
 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 23

the contaminant plume from reaching the western population.

Opening the hydrants in the isolated central region will flush out
the contaminant and protects the population in this area from
contaminated water as well. As Fig. 15 also indicates, this isolation
and flushing essentially reduces ingested contaminant mass to zero
for nearly 40% of the total population.
For the min-fH1 protocol, however, similar mechanics of isola-
tion and discharge are not clearly observed; while closure of CVS4
would protect the eastern region, no effective action is taken to
block further propagation of contamination across the western
region as Fig. 13 indicates. Identification of this particular
minimum-health-impact response protocol by the optimization
model is motivated by the expression of health impacts based upon
a preset and fixed threshold. Instead of attempting to protect the
whole population, the model tries to only decrease the number of
people who have ingested a mass of contaminant above that
threshold. This is partially achieved through distributing the
injected mass across a larger area to lower the concentration.
Fig. 14 indicates that the western population does not become
sick even when no response action is taken. Therefore, this popu-
lation does not contribute to the health impacts objective function
when fH1 formulation is used. This clearly explains why no attempt
is made to protect people in the western region for the min-fH1
protocol and health impacts practically remain unchanged. A
similar conclusion may be made for the people residing right on the
northeast downstream of the East WTP, where the health impacts
are noticeably more reduced when the min-fH2 protocol is
executed.

5. Conclusions

This study showed the formulation used for the quantification of

health impacts may have significant effect on the structure and
performance of optimal response protocols. A metric that is based
on a pre-specified exposure threshold has the advantage of being
understandable and explicit as it communicates the impacts
sensibly in terms of number of sicknesses or casualties when
Fig. 14. Spatial distribution of health impacts for Min-fH1 and Min-fH2 protocols for thresholds are known accurately. However, the results of this study
East Scenario and a 6-h response delay. Pipes are removed from the map for better
readability.
indicated that using the alternative metric of total ingested
contaminant mass may lead to more effective emergency response
recommendations even when exposure thresholds can be esti-
contaminant plume to the central area that is then partly flushed
mated accurately. In this sense, quantification of health impacts as
out through opening a number of hydrants there.
total mass of contaminant ingested by whole population may be
When it comes to protection of the western Mesopolis from the
preferable. This formulation also has the advantage that the
contaminant spread, however, the two response protocols based on
optimal response plan found for a specific contamination scenario
the number of sicknesses and total ingested contaminant mass
is still optimal for any other similar attack scenario with different
behave very differently. The min-fH2 protocol requires closing three
contaminant type and loading. This is because this metric is a linear
CVSs in the western area for a relatively long duration to prevent
function of injected mass, whereas the number of sicknesses and
deaths vary nonlinearly as injected mass changes. Ideally, dosee
response curves can be utilized after optimization is performed
based upon total ingested contaminant mass metric to provide a
more explicit and understandable presentation of adverse health
consequences.
It was shown that integrated response of contaminant
containment and discharge may effectively reduce the contami-
nation health impacts. The CVSs were located and configured
manually in this study. Future research may apply systematic ap-
proaches to automate and enhance this process for preparation of
more efficient optimal response protocols. Future work should
examine the incremental effectiveness of valve and hydrants
operation when they are executed in combination with other
response mechanisms such as public warnings disseminated
through news media. Public health can be best protected through a
Fig. 15. Distribution of health impacts for the minimum-health-impact protocols using comprehensive response plan that includes multiple response
formulations fH1 and fH2 for East Scenario, a 6-h response delay, and 5 hydrants. strategies.
24 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25

Although simultaneous contaminant containment and https://ceprofs.tamu.edu/kbrumbelow/Mesopolis/Developing%20Mesopolis.

pdf.
discharge operations may significantly reduce the health conse-
Jun, L., Guoping, Y., 2013. Iterative methodology of pressure-dependent demand
quences, it is observed that their performance exponentially de- based on EPANET for pressure-deficient water distribution analysis. J. Water
grades with delay time. This limits the direct application of the Res. Plan. Manag. 139 (1), 34e44.
computationally intensive simulation-optimization schemes dur- Kanta, L., Brumbelow, K., 2013. Vulnerability, risk, and mitigation assessment of
water distribution systems for insufficient fire flows. J. Water Res. Plan. Manag..
ing an emergency. Instead, these schemes may more appropriately http://dx.doi.org/10.1061/(ASCE)WR.1943-5452.0000281 (in press).
used to find optimal response protocols for an ensemble of repre- Kasprzyk, J.R., Reed, P.M., Characklis, G.W., Kirsch, B.R., 2012. Many-objective de
sentative contamination scenarios called the design basis threats. Novo water supply portfolio planning under deep uncertainty. Environ. Model.
Software 34, 87e104.
Clustering methods based on pattern analysis approaches are Labadie, J.W., Zheng, F., Wan, Y., 2012. Optimal integrated operation of reservoir-
currently being developed for the identification of such represen- assisted stormwater treatment areas for estuarine habitat restoration. Envi-
tative scenario. The set of threat-response information e also called ron. Model. Software 38, 271e282.
Laucelli, D., Berardi, L., Giustolisi, O., 2012. Assessing climate change and asset
the Response Planning Matrix in RPT (USEPA, 2003) e may be also deterioration impacts on water distribution networks: demand-driven or
used directly after an event happens in real-time or being fed as pressure-driven network modeling? Environ. Model. Software 37, 206e216.
input to dynamic simulation-optimization models for adaptation to Laucelli, D., Giustolisi, O., 2011. Scour depth modelling by a multi-objective evolu-
tionary paradigm. Environ. Model. Software 26 (4), 498e509.
the emerging conditions as an emergency proceeds. Efforts are Lindell, M.K., Prater, C.S., Perry, R.W., 2006. Fundamentals of Emergency Manage-
underway to develop these dynamic decision support models. ment. Federal Emergency Management Agency, Hyattsville, MD.
Liu, L., Ranjithan, S.R., Mahinthakumar, G., 2011a. Contamination source identifi-
cation in water distribution systems using an adaptive dynamic optimization
Acknowledgements procedure. J. Water Res. Plan. Manag. 137 (2), 183e192.
Liu, J., Yu, G., Savic, D., 2011b. Deficient-network simulation considering pressure-
dependent demand. In: International Conference on Pipelines and Trenchless
This material is based upon work supported by the National Technology (ICPTT), Beijing, China.
Science Foundation under Grant Number CMMI-0927739. Any Ostfeld, A., et al., 2008. The battle of the water sensor networks (BWSN): a design
opinions, findings, or recommendations expressed in this paper are challenge for engineers and algorithms. J. Water Res. Plan. Manag. 134 (6), 556e
568.
those of the authors and do not necessarily reflect the views of the Perelman, L., Ostfeld, A., 2010. Extreme impact contamination events sampling
National Science Foundation. We also gratefully acknowledge the for water distribution systems security. J. Water Res. Plan. Manag. 136 (1),
constructive and helpful comments of three anonymous reviewers 80e87.
Perelman, L., Ostfeld, A., 2011. Topological clustering for water distribution systems
that significantly improved this paper. analysis. Environ. Model. Software 26 (7), 969e972.
Perry, R.W., Lindell, M.K., 2007. Emergency Planning. John Wiley, Hoboken NJ.
Preis, A., Ostfeld, A., 2006. Contamination source identification in water systems: a
References hybrid model treeselinear programming scheme. J. Water Res. Plan. Manag. 132
(4), 263e273.
Afshar, A., Rasekh, A., Afshar, M.H., 2009. Risk-based optimization of large Preis, A., Ostfeld, A., 2008. Multiobjective contaminant response modelling for
flood-diversion systems using genetic algorithms. Eng. Optim. 41 (3), 259e water distributions systems security. J. Hydroinf. 10 (4), 267e274.
273. Rasekh, A., Afshar, A., Afshar, M.H., 2010. Risk-cost optimization of hydraulic
Alfonso, L., Jonoski, A., Solomatine, D., 2010. Multiobjective optimization of oper- structures: methodology and case study. Water Res. Manag. 24 (11), 2833e
ational responses for contaminant flushing in water distribution networks. 2851.
J. Water Res. Plan. Manag. 136 (1), 48e58. Rasekh, A., Brumbelow, K., 2012. Machine Learning Approach for Contamination
Baranowski, T.M., LeBoeuf, E.J., 2008. Consequence management utilizing optimi- Source Identification in Water Distribution Systems. World Environmental and
zation. J. Water Res. Plan. Manag. 134 (4), 386e394. Water Resources Congress, Palm Springs, CA.
Chick, S.E., Koopman, J.S., Soorapanth, S., Brown, M.E., 2001. Infection transmission Rasekh, A., Brumbelow, K., 2013. Probabilistic analysis and optimization to char-
system models for microbial risk assessment. Sci. Total Environ. 274 (1), 197e acterize critical water distribution system contamination scenarios. J. Water
207. Res. Plan. Manag. 139 (2), 191e199.
Cities of Bryan and College Station, 2005. B/CS Unified Design Guideline Manual, Reed, P.M., Hadka, D., Herman, J.D., Kasprzyk, J.R., Kollat, J.B., 2013. Evolutionary
Water, Sewer, and Streets. Bryan and College Station, Available www.bcsunited. multiobjective optimization in water resources: the past, present, and future.
net/2005Files/DesignManual2005.pdf. Adv. Water Res. 51, 438e456.
Davis, M.J., Janke, R., 2008. Importance of exposure model in estimating impacts Rossman, L.A., 2000. EPANET 2.0, User’s Manual. National Risk Management
when a water distribution system is contaminated. J. Water Res. Plan. Manag. Research Laboratory, U.S. EPA, Cincinnati, OH.
134 (5), 449e456. Savic, D.A., Bicik, J., Morley, M.S., 2011. A DSS generator for multiobjective optimi-
Deb, K., 2001. Multiobjective Optimization Using Evolutionary Algorithms. Wiley, sation of spreadsheet-based models. Environ. Model. Software 26, 551e561.
Chichester. Shang, F., Uber, J.G., Rossman, L., 2008. EPANET Multi-species Extension Software.
Deb, K., Agrawal, R.B., 1995. Simulated binary crossover for continuous search space. U.S. Environmental Protection Agency, Washington, DC. EPA/600/C-10/002.
Complex Sys. 9, 115e148. Shen, H., McBean, E.A., 2013. Application of parallel computing in data mining for
Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multi-objective contaminant source identification in water distribution systems. Can. Water
genetic algorithm NSGA-II. IEEE Trans. Evol. Comput. 6 (2), 182e197. Res. J. 38 (1), 34e39.
Giustolisi, O., Savic, D., 2010. Identification of segments and optimal isolation valve Shirzad, A., Tabesh, M., Farmani, R., Mohammadi, M., 2012. Pressure-discharge re-
system design in water distribution networks. Urban Water J. 7 (1), 1e15. lations with application in head driven simulation of water distribution net-
Giustolisi, O., Walski, T.M., 2011. Demand components in water distribution works. J. Water Res. Plan. Manag.. http://dx.doi.org/10.1061/(ASCE)WR.1943-
network analysis. J. Water Res. Plan. Manag. 138 (4), 356e367. 5452.0000305 (in press).
Gugat, M., 2012. Contamination source determination in water distribution net- Takbiri, Z., Afshar, A., 2012. Multi-objective optimization of fusegates system under
works. SIAM J. Appl. Math. 72 (6), 1772e1791. hydrologic uncertainties. Water Res. Manag. 26 (8), 2323e2345.
Hart, D.B., McKenna, S.A., Klise, K.A., Wilson, M.P., Murray, R., 2009. CANARY User’s Todini, E., Pilati, S., 1988. A gradient algorithm for the analysis of pipe networks. In:
Manual and Software Upgrades. U.S. Environmental Protection Agency, Wash- Coulbeck, B., Orr, Chun-Hou (Eds.), Computer Applications in Water Supply.
ington, D.C.. EPA/600/R-08/040A. Wiley: Research Studies Press, pp. 1e20.
Haxton, T., Murray, R., Klise, K., 2012. Examining the Application of Modeling Tools United States Congress, 2002. Public Health Security and Bioterrorism Preparedness
to Identify Effective Flushing Locations. World Environmental and Water Re- and Response act of 2002. Public Law No. 107-188, H.R. 3448, 107th Congress.
sources Congress, Palm Springs, CA. Government Printing Office, Washington, D.C.
Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. University of United States Environmental Protection Agency (USEPA), 2003. Response Protocol
Michigan Press, Ann Arbor, Mich. Toolbox: Planning for and Responding to Drinking Water Contamination
Janke, R.J., Murray, R., Haxton, T.M., Taxon, T., Bahadur, R., Samuels, W., Berry, J., Threats and Incidents. USEPA, Washington, D.C.
Boman, E., Hart, W., Riesen, L., Uber, J., Morley, K., 2009. Threat Ensemble United States Environmental Protection Agency (USEPA), 2004. Estimated per capita
Vulnerability Assessment-Sensor Placement Optimization Tool (TEVA-spot) water ingestion in the United States. USEPA, Office of Water, Washington, D.C.
Graphical User Interface. U.S. Environmental Protection Agency, Washington, Wagner, J.M., Shamir, U., Marks, D.H., 1988. Water distribution reliability: simulation
DC. EPA/600/C-08/010, 2009. methods. J. Water Res. Plan. Manag. 114 (3), 276e294.
Jin, Y., Branke, J., 2005. Evolutionary optimization in uncertain environments-a White, J.W., 1999. Hazards of Short-term Exposure to Arsenic-contaminated Soil.
survey. IEEE Trans. Evol. Comput. 9 (3), 303e317. Office of Environmental Health Assessment Services, Olympia, WA.
Johnston, G., Brumbelow, K., 2008. Developing Mesopolis: A ‘virtual City’ for Wu, Z.Y., Wang, R.H., Walski, T.M., Yang, S.Y., Bowdler, D., Baggett, C.C., 2006. Effi-
Research in Water Distribution Systems and Interdependent Infrastructures. cient pressure dependent demand model for large water distribution system
 A. Rasekh, K. Brumbelow / Environmental Modelling & Software 51 (2014) 12e25 25

analysis. In: 8th Annual Water Distribution Systems Analysis Symposium, Zechman, E.M., 2011. Agent-based modeling to simulate contamination events and
Cincinnati, Ohio. evaluate threat management strategies in water distribution systems. Risk Anal.
Yannopoulos, S., Spiliotis, M., 2013. Water distribution system reliability based on 31 (5), 758e772.
minimum cuteset approach and the hydraulic availability. Water Res. Manag. Zechman, E.M., Ranjithan, S., 2009. Evolutionary computation-based methods for
27 (6), 1821e1836. characterizing contaminant sources in a water distribution system. J. Water Res.
Zechman, E.M., Brumbelow, K., Lindell, M., Mumpower, J., Rasekh, A., Shafiee, M., Plan. Manag. 135 (5), 334e343.
2011. Including consumer behavior and system dynamics in vulnerability Zitzler, E., Thiele, L., 1998. Multiobjective Optimization Using Evolutionary Algo-
analysis for water distribution system contamination events. In: NSF CMMI rithms e a Comparative Case Study. Parallel Problem Solving from
Engineering Research and Innovation Conference, Atlanta, GA. NaturedPPSN V. Springer, Berlin/Heidelberg.