**2. Methodology**

In this study, a seismic reliability-based WDS optimization approach is proposed. Figure 1 shows the flowchart of the proposed WDS optimal design model which consists of two sub-models: seismic reliability estimation and optimization model. Seismic reliability estimation model (SREM) first generates stochastic earthquake events by using earthquake generation module. Here, the locations and magnitude of stochastic earthquakes are determined and peak ground acceleration (PGA) reached at each pipe is calculated (described in Section 2.1). Note that no hydraulic calculation is performed from this module. This module also determines the failure mode of pipe where each pipe is classified as normal, leakage, or breakage (Section 2.2).

**Figure 1.** Flowchart of the proposed water distribution system (WDS) optimal design model.

The failure modes are transferred to hydraulic calculation module where nodal pressures under the failure conditions (*i.e.*, the boundary conditions to the WDS system equation) are calculated (Figure 1 and Section 2.3). Pipe leakage and breakage are modeled by the emitter in EPANET. In order to better simulate the conditions, a semi pressure-driven analysis (semi-PDA) is proposed. After an initial run of EPANET (a hydraulic simulator for DDA), the negative-pressure nodes' demand is set to zero, and the second run is made (Section 2.4).

The resulting hydraulics are provided to seismic reliability calculation module to quantify seismic reliability which is defined as the ratio of available supply to the required demand under stochastic earthquake events. The proposed WDS optimal design model maximizes the system's seismic reliability with constraints on economic cost and pressure requirements (Figure 1 and Section 2.5). Harmony search algorithm (HSA) is used for optimization (Figure 1) [38,39]. In the last section of the Methodology, Todini [40]'s resilience indicator is defined. This indicator is not used in the proposed model but used for a postoptimization analysis.

The following subsections describe the details of the methodology, such as earthquake simulation, objective function, optimization approach, and semi-PDA method. Note that the following subsections are in the order that Figure 1 represents.

### *2.1. Earthquake Intensity Attenuation*

The proposed model considers the physical characteristics of the WDS's components to determine the failure modes given earthquake intensity. Therefore, supplemental information such as the types of pipes and surrounding soil characteristics should be provided. Network information such as the network layout, nodal demands and coordinates, pipe diameters and lengths, and sizes of the pump and tanks are entered from the EPANET input file. The nodal coordinates are used to calculate the location of a pipe that is used for calculating the distance from the epicenter.

Earthquake intensity is attenuated as it travels from the epicenter. In the model, three attenuation equations were adopted to quantify the damped energy of seismic waves. Kawashima *et al.* [41] relates the distance from the epicenter and earthquake magnitude to PGA (cm/s2) based on the 90 earthquakes that have occurred in Japan as

$$\text{PGA} = 403.8 \times 10^{0.265 \text{M}} \times (\text{R} + \text{30})^{-1.218} \tag{1}$$

where M = earthquake magnitude (e.g., M7 stands for earthquake magnitude 7), and R = the distance from the epicenter (km).

Lee and Cho [42] presented a similar PGA function that was developed using the historical data of South Korea as

$$
\log \text{PGA} = -1.83 + 0.386 \text{M} - \log \text{R} - 0.0015 \text{R} \tag{2}
$$

Instead of using the horizontal distance from the epicenter, Baag *et al.* [43] relate the Euclidean distance from focus ( Δ) and earthquake magnitude to PGA as

$$\ln \text{PGA} = 0.40 + 1.2 \text{M} - 0.76 \ln \Delta - 0.0094 \Delta \tag{3}$$

The average of the three calculated PGAs (from Equation (1) to Equation (3)) was used in our model.

### *2.2. Determination of Pipe Failure Mode*

The fragility of the pipes was determined based on the pipe repair rate (RR, number of pipe repairs/km) and the PGA. In this study, the RR equation suggested by Isoyama *et al.* [44] and adopted by ALA [19] was used. Isoyama *et al.* related the RR to the PGA by multiplying some correction factors denoting the pipe types, diameters and soil characteristics (Table 1) as

$$\text{RR} = \text{C1} \times \text{C2} \times \text{C3} \times \text{C4} \times 0.00187 \times \text{PGA} \tag{4}$$

where C1, C2, C3, and C4 represent the correction factors according to the pipe diameter, pipe material, topography, and soil liquefaction, respectively. Because of the lack of information, this study only considered the correction factor of the pipe diameter, C1, in Equation (4), while assuming other correction factors are equal to one. As seen in Table 1, for the smaller pipe, C1 increases and results in the higher probability of damage given the same PGA.

**Table 1.** Correction factors (Isoyama *et al.* [44]).


The post-earthquake pipe status was classified into normal, leakage, or breakage conditions. Pipe leakage is defined as a small rupture on the pipe wall or at the joint through which continuous minor water loss occurs. Pipe breakage is defined as a complete separation of the pipe into two pieces that causes complete loss of transportation ability. ALA [19] suggested that the probability of pipe breakage (Pfbreakage,*i*) is a function of the RR and the length of the pipe (Li, km) expressed as

$$P\_{\mathbf{f}\_{\text{bruakago},i}} = 1 - \mathbf{e}^{-\mathsf{RR} \times \mathsf{L}\_i} \tag{5}$$

As a rule of thumb, the probability of pipe leakage (Pfleakage,*i*) is assumed to be five times higher than the probability of pipe breakage [19]. If the cumulative probability of the two conditions is less than 1.0 (e.g., Pfbreakage,*i* = 0.1 and Pfleakage,*i* = 0.5), the complementary probability is assigned for the normal condition (Pnormal = 0.4). However, if the cumulative probability of the leakage and breakage is greater than 1 (e.g., Pfbreakage,*i* = 0.18 and Pfleakage,*i* = 0.90), the model assumes that the normal condition (without any damage) is not available, and the two failure probabilities are normalized to make the sum of 1.0 (Pfbreakage,*i* = 0.18/(0.18 + 0.90) = 0.17 and Pfleakage,*i* = 0.90/(0.18 + 0.90) = 0.83).
