Urban water supply systems’ resilience under earthquake scenario

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)³³. In sub-module 1, an input file (i.e., has the .INP extension) that contains descriptive data of the WDN is compiled in EPANET³⁴. 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.

Flowchart of earthquake simulation in WDN.

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⁷. 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.⁴ 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¹⁴. 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³⁵. 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³⁶. 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”³⁶. 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³⁷. 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³⁸. 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.⁴⁰; and the GMPE, Eq. 2, for the PGA proposed by Yu and Jin⁴¹ are adopted.

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³⁹, respectively. Attenuation models are developed for a specific location and earthquake. Therefore, correction factors are added to account for soil and pipe characteristics⁴². Table 3 and Eq. 5 display the correction factor and their categories and weights adopted from Isoyama et al.⁴².

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³⁴.

where P_n is the set of pipes connected to node n, q_p,n is water flow rate (m³/s) from pipe p into node n, D^act_n is the actual water demand leaving node n (m³/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³⁴, is selected in EPANET.

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³/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.⁴⁴. A formulation of this model is presented in Eq. 8.

where D is the actual demand delivered (m³/s), p is the pressure (Pa), P₀ is the pressure below which the customers cannot receive any water (Pa), D_f is the expected demand (m³/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⁴⁵ is used in WNTR to quantify leaks as a mass flowrate.

where d^leak is leak demand (m³/s), C_d is the discharge coefficient, A is leak hole surface area (m²), ρ is the water density (kg/m³), 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⁴⁶.

Network resilience

In engineering, resilience of a system is defined as the time required to restore an equilibrium state after a disruptive event occurs⁴⁷. 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⁴⁸. 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.

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₁, n₂, n₃, and n₄) that at a time t have an expected demand of (Ṽ_1t = Ṽ_2t = Ṽ_3t = Ṽ_4t = 0.1 m³), and an actual demand of (V_1t = 0.06 m³, V_2t = 0.09 m³, V_3t = 0.09 m³, and V_4t = 0.08 m³), then WSA_t = 0.8 or 80%. During a disruptive event, if the actual water volumes received dropped to (V_1t = 0.04 m³, V_2t = 0.06 m³, V_3t = 0.09 m³, and V_4t = 0.05 m³), 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.⁴, and industrial experience, R is set to 0.75 m³/capita/day and the threshold is 75%.

For the above given example, if q₁ = q₂ = q₃ = q₄ = 900 m³/day, then pop₁ = pop₂ = pop₃ = pop₄ = 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.