This section describes the methods engaged to simulate earthquake hazards on this specific WDN. Figure 1 depicts the four sub-modules built for this task in the water network tool for resilience (WNTR)^{33}. In sub-module 1, an input file (i.e., has the .INP extension) that contains descriptive data of the WDN is compiled in EPANET^{34}. In this file, system components and their characteristics are entered. For example, size, location, elevation, age, building materials, inlet and outlet type data are entered to describe storage tanks. For pipes, typical data includes location, slope, diameter, upstream and downstream joints, material, etc. Then, system components are assigned fragility curves as described in “WDN components failure” section.

Sub-module 2, house the earthquake simulation process. First, a pre-determined hypothetical scenario is cast where an earthquake location, magnitude, and depth are assigned. The values for these characteristics could be chosen randomly or, as in this study, based on knowledge of the system. For example, to predict a worst-case scenario, a highly probable magnitude, a particularly sensitive location (e.g., near pumping station), and a wider depth are selected. Second, proper Ground Motion Prediction Equations (GMPE) are assigned to represent the seismic wave attenuation. Location, topography, and simulated earthquake characteristics are the main considerations in this step. Third, a formulation of peak ground acceleration (PGA), peak ground velocity (PGV), and repair rate (RR) are inserted and computed for each scenario.

In the third sub-module, simulation of the WDN hydraulics is performed using the PDD approach (“Hydraulic analysis” section). Then, using fragility curves and the computed values of PGA, PGV, and RR, the status of each WDN component is determined. The last step of this sub-module is to quantify system damages. For example, if the shaking had caused a power outage in one pump and multiple leaks in tanks, joints, and pipes; subsequent pressure drops, and leak demand amounts, and locations for this damage are quantified in this step.

WDN responsivity and serviceability under earthquake hazard is investigated through resilience metrics in the fourth sub-module (“Network resilience” section). In this sub-module, resilience metrics are defined, formulated, and calculated. This sub-module may be considered as a summary of the simulation results which is very important for analysts and decision-makers.

As this simulation is a scenario with multiple outcomes, Monte Carlo simulation is used to tackle the uncertainty associated with multiple realizations and to estimate the probabilistic seismic reliability. For meaningful results and as a validation of the simulation, it is recommended to consider a large number of iterations^{7}. In this work, fifty realizations of each earthquake scenario (i.e., magnitude, location, and depth) are considered. The number of realizations was adopted from Klise et al.^{4} and can vary depending on the simulating machine capacity and speed.

### Earthquake occurrence

The effects of an earthquake on a WDN vary depending on earthquake intensity, location, and depth; topological characteristics; and WDN layout and facilities. When constructing an earthquake generation module, it is crucial to choose the input earthquake characteristics based on these factors^{14}. Indubitably, the input earthquake will differ for each study area. According to Hazards United States Multi-Hazard (HAZUS-MH), six different options could be used to define an earthquake module; some of which are customized for the United States data^{35}. Table 2 presents these options.

In this paper, the input earthquake events are extracted from historical data of the target WDN and expanded to arbitrary events. Historical data show that for a radius of 100 km, M-City has experienced earthquakes in the past four centuries with a maximum intensity of 4.1 M and a maximum depth of 18 km^{36}. Historical data also shows that about 70%, 20%, 5%, and 5% of earthquakes had an intensity between 2 and 3 M, less than 2 M, between 3 and 4 M, and more than 4 M, respectively. Historically, WDN in the study area have not yet experienced devastating damages to their components caused by seismic hazards. However, the occurrence of a higher intensity earthquake will most likely create a new case of risks and damages to the WDN in the region.

### Earthquake attenuation law

Earthquakes Canada defines an earthquake as “The sudden release of stored elastic energy caused by the sudden fracture and movement of rocks along a fault”^{36}. Some of the energy is released in the form of seismic waves, that cause the ground to shake. This ground shaking causes several ground motions depending on the transmitted energy path and geological characteristics of earth surface^{37}. Ground Motion Prediction Equation (GMPE) formulates this process of releasing energy from the earthquake epicentre to the earth’s surface. Two physical representations of this motion are velocity and acceleration of the seismic waves. When investigating the risks of seismic hazard, one should consider the peak of such motions (i.e., PGV and PGA). According to HAZUS-MH, vulnerabilities of buried pipes are related to PGV, and vulnerabilities of other WDN components (e.g., tanks, pumps, and water treatment plants) are related to PGA^{38}. Attenuation models are developed for a certain earthquake and a particular region. Therefore, modelers usually use the average of several empirical models to conclude a general behaviour^{4,7,39}. In this study, the GMPE, Eq. 1, for the PGV in a soil characteristic topology proposed by Kawashima et al.^{40}; and the GMPE, Eq. 2, for the PGA proposed by Yu and Jin^{41} are adopted.

$$PGV = 10^{{ – 0.285 + 0.711M – 1.85\left( {R + 17} \right)}}$$

(1)

$$PGA = 403.8\times 10^{0.265M} \left( {R + 30} \right)^{ – 1.218}$$

(2)

where M is the unitless earthquake magnitude, and R is the earthquake depth measured from the epicentre in kilometers (km). PGV is frequently used to estimate pipes repair rate (RR) which is defined as the number of repairs needed per one kilometer length of pipe (repairs/km). Equations 3 and 4 represent linear and power law RR models from American Lifelines Alliance^{39}, respectively. Attenuation models are developed for a specific location and earthquake. Therefore, correction factors are added to account for soil and pipe characteristics^{42}. Table 3 and Eq. 5 display the correction factor and their categories and weights adopted from Isoyama et al.^{42}.

$$RR = 0.00187\, \times \,PGV$$

(3)

$$RR = 0.00108\, \times \,PGV^{1.173}$$

(4)

$$C = C_{t} \, \times \,C_{d} { }\, \times \,{ }C_{l} { }\, \times \,C_{m}$$

(5)

### WDN components failure

Post-earthquake damages to the WDN are caused by ground motion. Fragility curves expressed in terms of ground motion functions are a common method used to predict damages to a WDN. Here, Fragility curves are statistical tools that predict the probability of a component reaching or exceeding a certain damage state under seismic excitation. WDN components such as tanks, pumps, and pipes can have different damage states. For example, pipe damage can be expressed in four states: breakage, major leak, minor leak, and no damage at all. A large database of earthquake characteristics and their damages to WDN components can be used to construct a type of fragility curve known as an empirical fragility curve. An example of this type of fragility curves is those reported in The American Lifelines Alliance reports^{39,43}. In these reports, PGV and RR are often used to estimate damages to buried pipes, where PGA is used in the case of pumps and tanks.

### Hydraulic analysis

In this study, WDN analyses are carried out in two compatible software environments, EPANET and Water Network Tool for Resilience (WNTR). EPANET has high reliability in hydraulic analysis of a steady-state WDN, and an approachable user interface. WNTR has a spectrum of capabilities that align well with this work: the ability to add disruptive incidents and response strategies, performing probabilistic simulations, and computing resilience metrics were essential to this study.

Throughout the WDN, hydraulic analyses are carried out for all nodes and links. Every node and link in the water network has a mass balance equation that quantifies inflow, outflow, and possible storage. Mass balance equations adopted in EPANET are shown in Eq. 6^{34}.

$$\mathop \sum \limits_{{p \in P_{n} }} q_{p,n} – { }D_{n}^{act} = 0\quad\, n{ } \in N{ }$$

(6)

where P_{n} is the set of pipes connected to node n, q_{p,n} is water flow rate (m^{3}/s) from pipe p into node n, D^{act}_{n} is the actual water demand leaving node n (m^{3}/s), and N is the set of all nodes. Here, q_{p,n} is positive unless water is outflowing node n into pipe p.

To account for head losses in links, conservation of energy formulae are used. In our hydraulic analyses, the Hazen-Williams head loss formula, Eq. 7^{34}, is selected in EPANET.

$$H_{nj} – { }H_{ni} { } = { }h_{l} = 10.667C^{ – 1.852} { }d^{ – 4.871} { }Lq^{1.852}$$

(7)

where h_{L} is the head loss in the pipe (m), C in the Hazen-Williams roughness coefficient (unitless), d is the pipe diameter (m), L is the pipe length (m), q is the water flow rate (m^{3}/s), H_{nj} and H_{ni} are are the heads at the starting and ending nodes (m), respectively.

During a disruptive incident like an earthquake, where pressure gradients in the WDN are unsteady, applying PDD methods in hydraulic simulations are more realistic and accurate than DDA. In WNTR, the PDD that is used in hydraulic simulations is that proposed by Wagner et al.^{44}. A formulation of this model is presented in Eq. 8.

$$D = { }\left[ {\begin{array}{*{20}c} { 0\quad p \le P_{0} } \\ {D_{f} \sqrt {\frac{{p{ } – { }P_{0} }}{{P_{f} – { }P_{0} }}{ }} \quad P_{0} < p < P_{f} } \\ {D_{f}\quad p \ge P_{f} } \\ \end{array} } \right]$$

(8)

where D is the actual demand delivered (m^{3}/s), p is the pressure (Pa), P_{0} is the pressure below which the customers cannot receive any water (Pa), D_{f} is the expected demand (m^{3}/s), and P_{f} is the pressure above which the expected demand should be received (Pa).

Another model that is crucial during a disruptive event is the leak model. Leaks are expected to occur in tanks, fittings, junctions, and pipes. Leaks negatively affect the hydraulics and serviceability of the WDN. In WNTR, leaks or leak demand can be modeled, and lost water can be quantified between the start of the disruptive event and the time of repair for those affected components. In cases where leaks are major, a breakage can be modelled by splitting a pipe into two sections and adding two new disconnected junctions at both ends. The equation, expressed here in general form as Eq. 9^{45} is used in WNTR to quantify leaks as a mass flowrate.

$$d^{leak} = C_{d} A \sqrt {2\rho P^{\alpha } }$$

(9)

where d^{leak} is leak demand (m^{3}/s), C_{d} is the discharge coefficient, A is leak hole surface area (m^{2}), ρ is the water density (kg/m^{3}), P is the in-pipe gauge pressure (Pa), and α is a correction factor. Assuming a turbulent flow and large leaks out of pipes, values for C_{d} and α are set to 0.75 and 0.5, respectively^{46}.

### Network resilience

In engineering, resilience of a system is defined as the time required to restore an equilibrium state after a disruptive event occurs^{47}. The WDN resilience cycle begins in the design stage, then, is shaped through operations and maintenance, and finally tested through repair and mitigation processes. Each stage aims to reduce the effect of disruptive events when they occur. Quantifying WDN resilience is crucial to predict the possible system response to a variety of disruptive events. The United States Environmental Protection Agency (USEPA) classifies WDN resilience metrics into five categories: economic, hydraulic, topographic, water quality, and water security^{48}. Selecting the most representative metric is important and depends on the disruptive event scenario, and available input data. In this study, hydraulic metrics are employed to measure the WDN resilience. This includes, pressure, demand, water serviceability, and population impacted. Earthquake expected damages to the WDN (e.g., leaks, breakages, power loss) and available hydraulic model data were the reason hydraulic metrics were chosen. Equations 10 and 11 compute water serviceability and population impacted, respectively.

$$WSA_{t} = \left( {\mathop \sum \limits_{n = 1}^{n = N} V_{nt} } \right)/\left( {\mathop \sum \limits_{n = 1}^{n = N} \tilde{V}_{nt} } \right)$$

(10)

$$Pl_{t} = \mathop \sum \limits_{n = 1}^{n = N} pop_{n} *d_{nt}$$

(11)

$$pop_{n} = q_{n} /R$$

(12)

$$d_{nt} = { }\left\{ {\begin{array}{*{20}l} 1 \hfill & { if\,\, V_{nt} /\tilde{V}_{nt} < threshold} \hfill \\ 0 \hfill &\quad\quad { otherwise} \hfill \\ \end{array} } \right\}$$

(13)

where WSA_{t} is the water serviceability of the WDN at time t, n is the node number of the total N network nodes, V_{nt} and Ṽ_{nt} are the actual and expected water volume received at node n at time t. For example, if we have a simple network of four nodes (n_{1}, n_{2}, n_{3}, and n_{4}) that at a time t have an expected demand of (Ṽ_{1t} = Ṽ_{2t} = Ṽ_{3t} = Ṽ_{4t} = 0.1 m^{3}), and an actual demand of (V_{1t} = 0.06 m^{3}, V_{2t} = 0.09 m^{3}, V_{3t} = 0.09 m^{3}, and V_{4t} = 0.08 m^{3}), then WSA_{t} = 0.8 or 80%. During a disruptive event, if the actual water volumes received dropped to (V_{1t} = 0.04 m^{3}, V_{2t} = 0.06 m^{3}, V_{3t} = 0.09 m^{3}, and V_{4t} = 0.05 m^{3}), then WSA_{t} would correspondingly drop to 0.6 or 60%.

In Eq. 11, Pl_{t} is the population impacted at time t, pop_{n} is the population at node n, d_{nt} is a binary variable, q_{n} is the average water volume consumed per day at node n under normal conditions, and R is the average water volume consumed per capita per day. Based on Klise et al.^{4}, and industrial experience, R is set to 0.75 m^{3}/capita/day and the threshold is 75%.

For the above given example, if q_{1} = q_{2} = q_{3} = q_{4} = 900 m^{3}/day, then pop_{1} = pop_{2} = pop_{3} = pop_{4} = 1200 capita. During the given disruptive event, V_{1t}/ Ṽ_{1t} = 0.4, V_{2t}/ Ṽ_{2t} = 0.6, V_{3t}/ Ṽ_{3t} = 0.9, V_{4t}/ Ṽ_{4t} = 0.5, and d_{1t} = 1, d_{2t} = 1, d_{3t} = 0, d_{4t} = 1. As a result, the total population impacted at time t, Pl_{t} would equal to 1200 + 1200 + 1200 = 3600 capita.

Support Lumiserver & Cynesys on Tipeee

**Visit
our sponsors**

Wise (formerly TransferWise) is the cheaper, easier way to send money abroad. It helps people move money quickly and easily between bank accounts in different countries. Convert 60+ currencies with ridiculously low fees - on average 7x cheaper than a bank. No hidden fees, no markup on the exchange rate, ever.

**
Now you
can get a free first transfer up to 500£
with your ESNcard. You can access this offer here.**

Source link