Risk-based assessment of the post-earthquake functional disruption and restoration of distributed infrastructure systems

Journal Pre-proof Risk-Based Assessment of the Post-Earthquake Functional Disruption and Restoration of Distributed Infrastructure Systems Agam Tomar, Henry V. Burton PII:

S2212-4209(20)31504-1

DOI:

https://doi.org/10.1016/j.ijdrr.2020.102002

Reference:

IJDRR 102002

To appear in:

International Journal of Disaster Risk Reduction

Received Date: 18 August 2020 Revised Date:

28 November 2020

Accepted Date: 29 November 2020

Please cite this article as: A. Tomar, H.V. Burton, Risk-Based Assessment of the Post-Earthquake Functional Disruption and Restoration of Distributed Infrastructure Systems, International Journal of Disaster Risk Reduction, https://doi.org/10.1016/j.ijdrr.2020.102002. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Elsevier Ltd. All rights reserved.

1

Risk-Based Assessment of the Post-Earthquake Functional Disruption and

2

Restoration of Distributed Infrastructure Systems

3

Agam Tomar1 and Henry V. Burton2 1

4

Ph.D. Student, Department of Civil and Environmental Engineering, University of California Los

5 6

Angeles, [email protected] 2

Associate Professor, Department of Civil and Environmental Engineering, University of California Los

7

of

Angeles, [email protected] Abstract

9

A risk-based framework for assessing the post-earthquake functional loss and restoration of distributed

10

infrastructure systems is presented. Using an appropriate seismic source model, a stochastic earthquake

11

catalogue is generated for the target region. For each event in the catalogue, the joint distribution of

12

shaking and component-level damage maps are produced. Process-based discrete event simulation is

13

coupled with an appropriate method for linking physical damage to functional impacts to model the

14

restoration for each damage map. By considering all possible events that are impactful to the target region

15

and their associated rates of occurrence, the full probability distribution of the functional recovery

16

trajectory of the system is obtained. Dispersion disaggregation is used to isolate the individual

17

contributions of multiple random variables to the uncertainty bounds in the restoration trajectory. A case

18

study involving the city of Napa’s water system is used to demonstrate how the framework can support

19

resilience-based risk-informed decision-making for distributed infrastructure.

20

Keywords: infrastructure system functional recovery; risk analysis; stochastic event set; dispersion

21

disaggregation; community resilience planning

22

1. Introduction

23

Resilience can be described as the ability of a system to mitigate the effects of adverse events and adapt

24

and recover in a timely manner. Models that quantify the immediate functional disruption and restoration

25

of distributed infrastructure systems following a hazard event can inform resilience planning decisions.

Jo

ur

na

lP

re

-p

ro

8

One major challenge in this regard is the selection of the event(s) to use as the basis of the assessment.

27

For spatially distributed infrastructure systems, it is often impossible to find one such event. Additionally,

28

even if a single “worst-case” (in terms of consequences) scenario is determined, it is important to consider

29

its probability of occurrence. Risk-based assessment methodologies are useful in this regard because they

30

consider the full spectrum of possible consequential scenarios and their associated occurrence

31

probabilities (or rates).

32

Recent efforts to model the earthquake-induced loss and functional restoration of distributed

33

infrastructure systems have embraced simulation-based approaches (e.g. [1–6]). A common theme among

34

these studies is that only a single event was considered in their assessments. Single-scenario seismic

35

evaluations are desirable because they require less computational effort and can be useful under some

36

circumstances. However, for other applications (e.g. insurance purposes, resilience planning), it is

37

necessary to estimate the effects of many or even all possible impactful scenarios.

38

The goal of seismic risk assessment for spatially distributed infrastructure is to estimate the annual rate of

39

exceedance or return period of different levels and types (i.e. different performance metrics) of impacts.

40

Stochastic event-set models, which consider most or all significant events and their annual rates of

41

occurrence, have been used for this purpose (e.g. [7]). While there have been prior studies on stochastic

42

event-set or risk-based assessment of distributed infrastructure systems (e.g. [8–10]), they all utilized

43

immediate-impact performance metrics i.e. functional restoration of the system was not considered.

44

This study presents a methodology for risk-based assessment of post-earthquake functional disruption and

45

restoration of distributed infrastructure systems. The uniqueness of the proposed framework relative to

46

prior studies is the ability to consider a stochastic event catalogue in the context of post-earthquake

47

functional restoration modeling of spatially distributed infrastructure. In addition to this primary

48

contribution, we perform a disaggregation of the different sources of dispersion that affect the stochastic

49

functional recovery trajectories. Lastly, we propose a multi-metric optimization method for selecting

50

event subsets for infrastructure risk assessment. Through a case study involving the City of Napa’s water

51

distribution system, we demonstrate how the proposed framework provides a comprehensive

Jo

ur

na

lP

re

-p

ro

of

26

52

understanding of the impact to stakeholders and informs risk-mitigation and resilience-enhancement

53

planning.

54

2. Overview of Risk-Based Assessment Framework for Modeling Post-Earthquake Functional

55

Recovery of Distrbuted Infrastructure Figure 1 shows an overview of risk-based post-earthquake functional disruption and restoration

57

assessment framework for distributed infrastructure. First, a stochastic earthquake catalogue comprised of

58

scenarios is generated for the region of interest. The events vary based on magnitude, epicentral

59

location, rupture parameters (e.g. method of rupture, rupture type) and annual rate of occurrence. For each

60

event, a set of

61

correlation models. For each ground motion map,

62

(e.g. pipes, tanks and pump stations in a water distribution system) are sampled using the associated

63

fragility function and the shaking intensity at the relevant location. The ground motion map used for each

64

damage simulation is randomly sampled from the set of

65

of events and ground motion fields, a total of

66

disruptive effects of damage to the affected component needs to be specified. Using a water distribution

67

system as an example, the functional impact of physical damage to pipes is the occurrence of a leak or

68

break, which leads to the loss of water, whereas damage to pumps and tanks leads to a loss of hydraulic

69

power. Component repairs are simulated on each damage map using a process-based discrete event

70

simulation (PBDES) model, which considers the available resources (inspection, repair crew and

71

materials) and the sequencing (e.g. damaged component must be inspected before it is repaired) and

72

duration of the associated activities.The metric and model for simulating functional disruption and

73

restoration at the system level is then defined. Both the metric and model depends on the system being

74

considered. For a bridge network, the metric can be based on the movement of goods and services

75

throughout the system. For a water network, which is used to demonstrate the proposed methodology, the

76

performance of the system can be assessed through time-dependent hydraulic simulation which provides

77

some measure of the ability of the network to meet the service demand. Once the recovery trajectories for

ro

of

56

-p

ground motion maps are produced using the appropriate ground motion and spatial

=

maps. In other words, considering the number ×

damage maps are produced. Next, the

Jo

ur

na

lP

re

realizations of damage to individual components

the chosen metric are generated, several static (e.g. robustness) and temporal (e.g. time to restore 100% of

79

the pre-event functionality) resilience-related metrics are quantified. By integrating the metric values

80

associated with a given damage map with the occurrence rate of the event that produced it, the full

81

probability distribution of the system performance is obtained.

lP

re

-p

ro

of

78

na

82 83

ur

Figure 1 Overview of end-to-end simulation framework 2.1 Hazard Characterization

85

The hazard characterization step produces a set of maps with ground-motion intensity realizations at each

86

location of interest and their corresponding occurrence rates. Each map captures the joint probability

87

distribution of the spatially distributed ground motion intensities. First,

88

generated from a seismic source model, which also gives the annual rate of occurrence for each event

89

with a specified magnitude, location, and faulting type. Then, an empirical ground motion model (GMM)

90

is used to compute the resulting intensity (e.g. peak ground acceleration (PGA) and peak ground velocity

91

(PGV)) for a given event and location of interest. The GMM predicts the log-mean

92 93

Jo

84

,

,

,

,…

, within-event

ground-motion intensity for the

site

and between-event = 1, 2, … ,

and the !

earthquake scenarios are

residual standard deviations of the ! = 1, 2, … ,

earthquake scenario.

94 95 96

represents the moment magnitude of the ! distance from the surface projection of the ! ,

earthquake scenario,

denotes the closest horizontal

earthquake scenario’s fault plane to the

represents the time-averaged shear wave velocity to a depth of 30$ at the

location, and

location. For each

97

earthquake scenario,

98

obtained using an appropriate model. After sampling the residuals, the total natural log of the ground

99

motion intensity ( ) is computed using Equation 1:

realizations of the spatially-correlated ground-motion intensity residual terms are

'

=

,

,

,

,… +

') '

+ ' *'

(1)

of

ln

where the index + corresponds to the ground-motion map (+ = 1,2, … ,

,ℎ./.

101

are the normalized within-event and between-event residuals for

102

location variability and *' is the event-to-event variability. Both ) ' and *' are standard normal random

103

variables. The vector of ) ' is modeled as a spatially correlated multivariate normal distribution and *' is

104

a standard univariate normally distributed random variable. The result is a set of

105

Since an equal number

106

occurrence of the earthquake scenario normalized by

107

2.2 Damage modeling

108

Calculating network performance requires an assessment of the earthquake-induced physical damage to

109

the relevant components of the network. The probabilistic link between the ground-motion intensity and

110

physical damage is provided by fragility functions in the form of 0 123 ≥ 56|

111

discrete random variable whose value represents the damage state for the 9

112

predefined threshold. The damage state is conditioned on a realization, 8, of the random variable

113

which, in this study, is the shaking intensity at the

114

2.3 Restoration Modeling

115

The DES approach [11] is used to conduct the restoration analysis. The core elements of a DES model are

116

entities, attributes, events, resources and time. Entities represent specific objects within the system, which

117

have attributes, experience events, and consume resources over time. An attribute is a feature whose

ro

100

=

× ), ) ' and *'

lP

re

-p

, respectively. ) ' is the location-to-

ground-motion maps.

is represented as ,' .

Jo

ur

na

of ground-motion maps is produced for each earthquake, the annual rate of

site corresponding to the +

'

= 8 where 123 is a

component and 56 is a ',

ground-motion map.

accumulation over time defines the state of an entity. Events affect the state of entities, which are serviced

119

by resources. The interaction between entities, events, processes, and resources, takes place in an

120

environment. The PBDES approach utilizes a set of processes, which represent a sequence of events and

121

activities through which a specific object moves. In the PBDES restoration model, the network

122

components (e.g. for a water network: pipes, tanks, pump stations, and for power network: substations,

123

transformers, transmission lines) are the entities, their damaged and/or functional states are the attributes,

124

and the inspection and repair crews are the resources. The processes include network component

125

inspection, repair and crew changes.

126

Given the damage map, the first step in the PBDES network restoration model is the initiation of

127

component inspections, which is modeled as a process. Inspection crew members (a type of resource)

128

travel to component locations within the network based on a set of predefined priorities and perform

129

inspections over some duration. The start and end of the inspection process represent two events. If an

130

inspected component is found to be damaged, it is added to a repair queue. Once the component gets to

131

the “front” of the queue, its repair takes place over some duration. Similar to inspection, the beginning

132

and end of the repair process for a component are represented as events.

133

2.4 Network Performance

134

The DES model captures the restorative change (through repairs) in the physical state of the system over

135

time. However, more important than the physical state is the evolution of system functionality during the

136

period following the disruptive event. Therefore, a model is needed to link the physical and functional

137

state of the system over time where the latter is described using an appropriate metric. As noted earlier,

138

the chosen metric and model depends on the system that is being investigated [12]. For road networks,

139

two common performance metrics are connectivity [13,14] and flow capacity [15]. Power network

140

performance can be measured using connectivity [16], serviceability ratio [9] and system flow [17].

141

Hydraulic power capacity [18], minimum surplus head [19], connectivity [20,21], and entropy-based

142

performance measures have been used to measure water network performance in prior studies.

Jo

ur

na

lP

re

-p

ro

of

118

Given the functional restoration trajectory of the system for a given event or ground motion map, various

144

resilience-related metrics can be used to summarize the performance of the system. The robustness of the

145

system can be assessed by quantifying the immediate loss of functionality :;< where : is the

146

functionality measure. Using the same restoration curve, the rapidity of the system can be measured by

147

computing the time to restore the system to some percentage of the pre-even functionality =>%@ . The

148

cumulative loss of functionality over some predefined period :A< is taken as the normalized area above

149

the restoration curve. After calculating the network performance for all simulated ground-motion maps,

150

the annual exceedance rates corresponding to each metric B Ĉ

of

143

ro

is calculated using Equation 2

'∈J

(2)

-p

B Ĉ = E ,' F[C' ≥ Ĉ ] where Ĉ a specified network performance level,

is the set of all damage maps, ,' is the weight (annual

152

occurrence rate) of damage map +, C' signifies the network performance resulting from damage map +,

153

and F[∎] is an indicator function that is equal to 1 when the network performance resulting from damage

154

map + exceeds the specified level, and 0 otherwise.

155

3. Case Study

156

Figure 2 summarizes the main steps involved in applying the risk-based assessment frameweork to the

157

Napa water network. The details of each step and the associated results are presented in Sections 3.1 to

158

3.5.

Jo

ur

na

lP

re

151

of ro -p re

159

Figure 2 Application of the risk-based assessment framework to the Napa water system

lP

160

3.1 Description of the Napa Water Network

162

Figure 3 shows the key components of the Napa water network, which supplies water to a population of

163

approximately 88,000 people via 25,000 service connections.

164

approximately 612 km of pipeline (~7400 pipe segments), 15 storage tanks with a total of 113, 652 m3

165

storage, 10 pumping stations, and 14 pressure regulating stations, represent the main components of the

166

network. The ten pumping stations serve approximately 10% of the population, and the remaining 90% is

167

served by gravity flow coming from the water treatment plants. All potable water comes from three

168

sources: Lake Hennessey (38,237,880 m3), the Milliken Reservoir (1,726,872 m3), and the Storage Water

169

Project (SWP) (27,013,212 m3). However, the Milliken Reservoir is only used during the high-demand

170

summer months [22]. Readers are referred to [2] for a detailed discussion on the Napa water network.

Three water treatment plants,

Jo

ur

na

161

of ro -p re lP na

171 172

ur

Figure 3 Schematic representation of the Napa water network 3.2 Hazard Characterization

174

The suite of possible earthquake scenarios that would affect the system is generated using the UCERF2

175

model [23–25]. The region is divided into 285 (1.45 km X 1.85 km) (1 minute) grids, which overlay the

176

water network, and the center of each grid is taken as a location of interest. The OpenSHA IM Event Set

177

Calculator [26] is used to generate a total of

178

are used to obtain the log-mean and within-event and between-event residuals for PGA and PGV (e.g.

179

[27–31]). For PGA, a weight of 0.12 is applied to the Idriss model and 0.22 to all others as recommended

180

in the 2014 U.S. National Seismic Hazard Maps and for PGV, an equal weight of 0.25 is applied to all

181

models except Idriss [32]. The correlation model for spatial ground motion intensities [33] is used to

182

generate

183

scenario. The resulting event set consists of

Jo

173

= 2,343 scenarios. For each earthquake scenario, GMMs

= 50 realizations of the spatially correlated ground-motion intensity residuals for each = 117,150 PGA and PGV maps. Figure 4 shows the

median shake maps considering all scenarios. The associated PGA and PGV values range from 0.0317 g

185

to 0.0741 g and from 2.83 cm/s to 10.66 cm/s, respectively.

(a)

(b)

re

186 187

-p

ro

of

184

Figure 4 Spatial distribution of median (a) PGA (in g) and (b) PGV (in cm/s) values considering all shake maps

190

3.3 Damage Assessment

191

A damage map is generated by sampling the damage state for each component, with probabilities

192

obtained from the fragility functions conditioned on the shaking intensity. For pipes, the fragility

193

functions developed by [2] are used to estimate the damage state conditioned on PGV and the repair rate

194

(RR). For tanks and pump stations, fragility curves from HAZUS-MH [34] are used to estimate the

195

damage state conditioned on PGA. For each scenario, a spatially correlated ground motion map is chosen

196

randomly from the set of 50, and the corresponding damage map is obtained. This process is repeated

197

1000 times for each scenario resulting in total of

198

realizations for each scenario is chosen to minimize the dispersion in network component damage across

199

different Monte Carlo estimates.

200

Of the 2,343 scenarios, only 135 caused damage to the water network. These 135 scenarios occurred on

201

the Hayward-Rodgers Creek, North San Andreas, Great Valley, Hunting Creek-Berryessa, Maacama-

202

Garberville, Mount Diablo Thrust and West Napa fault. Figure 5 shows the magnitude distribution of the

203

damaging scenarios. Most magnitudes range between 6.75 – 7.0 followed by 7.0 – 7.25. The median PGA

Jo

ur

na

lP

188 189

= 2,343,000 damage maps. The number of

(in g) and PGV (cm/s) across only the damaging maps are shown in Figure 6. The median PGA values

205

range from 0.092 g to 0.211 g and the median PGV values range from 8.37 cm/s to 34.94 cm/s. The

206

spatial distribution of mean pipe repair rate (RR) per km considering only the damaging events is shown

207

in Figure 7. A M 6.7 on the West Napa fault resulted in the highest number of damaged pipes. With an

208

annual occurrence rate of 8.288E-5, this event damaged 559 pipes (458 leaks and 140 breaks), 11 tanks

209

(10 with slight damage and 1 with moderate damage). There was no damage to pump stations.

na

lP

re

-p

ro

of

204

210

Figure 5 Magnitude distribution for damaging scenarios

Jo

ur

211

212 213

(a)

(b)

214

Figure 6 Spatial distribution of median (a) PGA (in g) and (b) PGV (in cm/s) values considering the 135

215

scenarios that caused damage

of ro -p

Figure 7 Mean pipe repair rate (per km) distribution considering only the damaging scenarios

218

3.4 Restoration Modeling and Network Performance

219

The DES model is used to simulate the functional recovery of the Napa water system. The uncertainty in

220

the inspection and repair processes are captured by assigning probability distributions to the number of

221

crews and event durations. The initial number of inspection and repair crews are sampled from a uniform

222

distribution whose values range from 10 to 120 and 20 to 40, respectively. Similarly, the initial number of

223

inspection and repair crews for tanks and pumps are sampled from a uniform distribution ranging from 10

224

to 20 and 2 to 10, respectively. Based on discussions with the Napa Water Division [22], the reduction

225

schedule for the inspection and repair crews shown in Table 1 are used. A triangular distribution is used

226

for the pipe inspection and repair durations [35], which are shown in Table 2. A normal distribution is

227

used for the tank and pump stations repair times (Table 3) [34]. Component inspections and repairs are

228

constrained by the available number of crews and working hours (9 a.m. – 6 p.m.). The inspection and

229

repair event durations and the initial distribution of crew members are taken from previous studies

230

[2,4,34].

Jo

ur

na

lP

re

216 217

Table 1: Adopted inspection and repair crew reduction schedule

Days After the Earthquake Day 0 Day 5 Day 10 Day 15 Day 20 Day 30

232 233

Number of Number of Inspection Crews Repair Crews 10 100 10 90 6 72 6 42 3 21 2 8

Min

Mode

Max

0.5 hour

1 hour

2 hours

3 hours 4 hours 4 days 6 days

4 hours 6 hours 4 days 8 days

6 hours 12 hours 6 days 10 days

lP

re

-p

Event Trunk or Distribution Damage Inspections Tank Pump Distribution Leak Distribution Break Repairs Trunk Leak Trunk Break

of

Table 2: Event duration triangular distribution parameters

ro

231

234

na

Table 3: Event duration normal distribution parameters: mean and standard deviation (σ) Event slight/minor moderate extensive complete slight/minor moderate extensive complete

ur

235

Jo

Pumping Plants

Repairs

Water Storage Tanks

236

Mean (days) 0.9 3.1 13.5 35.0 1.2 3.1 93.0 155.0

σ (days) 0.3 2.7 10.0 18.0 0.4 2.7 85.0 120.0

237

3.5 Hydraulic Simulation

238

The restoration process for the water distribution system seeks to satisfy the water demand for households

239

and other critical infrastructure. Because the system comprises a complex network of components,

240

individual repairs do not always increase the overall water supply. To ensure that the model reflects this, a

241

pressure-driven analysis is performed at each restoration time step using WNTR to quantify the

242

availability of water at each demand node. The output from this simulation combined with the

component’s recovery enables stakeholders to make informed decisions and design more effective pre-

244

mitigation strategies. Damage to a pipe is modeled as a leakage discharge, which assumes that there is a

245

hole in the pipe with a specified diameter. Similarly, tank damage is modeled by adding a leak at the base

246

with an area that is defined in terms of the tank diameter. When the pressure in the tank is reduced (as a

247

function of the leak area) below some threshold, it is taken out of operation. Damage to pump stations are

248

represented as an outage, during which the pump is not operational.

249

By coupling the PBDES model with the WNTR platform, an SSI restoration curve is generated for each

250

damaging ground motion map. A total of 2,343,000 SSI restoration curves are produced. Figure 8 shows

251

a histogram of 22O;< based only on the damaging events. The 22O;< values range from approximately 0%

252

to 88%. The 59,886 damage maps that resulted in water service disruption accounts for approximately 2.5%

253

of the stochastic event set. A histogram of the time to restore 100% SSI is shown in Figure 9. Similar to

254

Figure 8, the vertical axis values (fraction of cases) are only based on the damaging events. Note that the

255

bins of the histogram are unevenly distributed to highlight the periods of interest that are commonly used

256

by decision-makers (e.g. 3 days, 1 week, 2 weeks, 1 month, 3 months, 6 months following the

257

earthquake). The =P

258

the worst-case scenarios, the water service can remain disrupted for up to 7 months. The Applied

259

Technology Council 2016 (ATC) suggests a performance target of 100% water service restoration within

260

3 days following a “probable” earthquake [36]. It is worth noting that 2, 246, 016 damage maps (~95.7%

261

of all damage maps) meet this criterion. Figure 10 shows a histogram of the 22OA< , which is taken as the

262

area above the restoration curve. A 100% value corresponds to a complete loss of service over the

263

considered time period and a 0% value means no loss. The 22OA< values range from 0% to approximately

264

56% with a mean value (considering only the damaging events) of 1.31%.

values in Figure 9 range from 0 to 220 days, which indicates that, for some of

Jo

ur

%QQ;

na

lP

re

-p

ro

of

243

-p

ro

of

265

Figure 8 Histogram showing the distribution of 22O;< for the damaging events

lP

267

re

266

Jo

ur

na

268

269 270

Figure 9 Histogram showing the distribution of =P

%QQ;

for the damaging events

ro

of

271

-p

272

Figure 10 Histogram showing the distribution of 22OA< for the damaging events

re

273

3.6 Disaggregation of Dispersion by Source of Uncertainty

275

Disaggregation of the network performance dispersion is an effective way to identify which random

276

variables contribute the most to the uncertainty bounds of the predicted outcomes. Within the end-to-end

277

simulation framework, there are several sources of uncertainty embedded in different stages (Figure 11).

278

In the hazard characterization, there is uncertainty in the earthquake scenario (addressed by considering

279

all possible events) and the spatial distribution of shaking intensity conditioned on a single scenario

280

(addressed by generating multiple spatially correlated ground motion maps). In the damage assessment,

281

multiple maps are generated to propagate the uncertainty embedded in the fragility functions. Within the

282

PBDES restoration model, the initial number of inspection and repair crews, temporal input parameters

283

(i.e., inspection and repair times) and the availability of repair crews, are all considered as random

284

variables.

285

Four sources of uncertainty are considered in the disaggregation: spatial distribution of shaking

286

conditioned on the scenario, component damage fragility, the available crew members and the process

287

durations. For a given scenario, the end-to-end simulation is performed by considering the uncertainty in

Jo

ur

na

lP

274

the random variable of interest while using the median values for all others. An M 6.85 event occurring

289

on the West Napa fault with an annual rate of occurrence of 3.01E-5 is considered for this purpose.

of

288

ro

290

Figure 11 Sources of uncertainty considered in the risk-based assessment framework

292

To quantify the uncertainty associated with each type of parameter, 100 end-to-end simulations are

293

performed. A total of 100 end-to-end simulations is chosen as the variance converges after approximately

294

60 simulations for all the network performance metrics. The mean plus and minus one standard deviation

295

SSI curves when all sources of uncertainty are considered is shown in Figure 12 and Table 4 shows the

296

associated mean R and coefficient of variation S values of the three network performance metrics.

297

=P

298

22O;< has the lowest S = 0.18 . This result suggests that the dispersion associated with metrics related to

299

recovery time are generally higher than the ones that only consider the initial loss of service. Additionally,

300

the metric that considers both the initial loss of service and recovery time have a higher dispersion than

301

initial loss of service metric but lower than the recovery time metric.

na

lP

re

-p

291

ur

has the highest dispersion value as measured by the coefficient of variation S = 0.82 and

Jo

%QQ;

302

Table 4 Mean R and coefficient of variation S of the network performance metrics for a M 6.85

303

earthquake scenario with all sources of uncertainty considered

Performance Metric

VVWWX Z[\\%VVW VVWYX

Statistic Value 65.6% 0.18 71.9 days 0.82 13.1% 0.69

307

na

Figure 12 Mean and one-standard deviation bounded SSI curve for the M 6.85 earthquake scenario with

ur

306

all sources of uncertainty considered

Jo

305

lP

re

-p

ro

of

304

Statistic Type μ δ μ δ μ δ

308

The 100 end-to-end simulations are repeated while considering the four sources of uncertainty (ground

309

shaking, component damage, crew size and duration parameters) individually. Figure 13 shows the

310

associated SSI curves and the mean R and coefficient of variation S values for the time required (in

311

days) to reach 25%, 50%, 75%, 100% of the pre-event SSI levels are shown in Table 5. The coefficient of

312

variation values in Table 5 and plots in Figure 13 indicate that the crew size contributes the most to the

313

overall uncertainty bounds for restoration times associated with higher pre-event SSI values (75% and

314

100%). Whereas for lower pre-event SSI levels (25% and 50%), the overall uncertainty in the restoration

315

time is dominated by the ground shaking intensity. In general, the dispersion in the restoration time

316

increases with the pre-event SSI level. For the process duration variables, Figure 13 shows that most of

317

the uncertainty is in the tail end of the restoration curve. Further evidence of this can be found by

318

comparing the coefficient of variation for the 100% SSI pre-event restoration time (0.88) to the next

319

highest value (0.06), which corresponds to the 75% SSI level.

320

Table 5 Mean R and coefficient of variation S values for the time to restore 25%, 50%, 75%, and 100%

321

of the pre-event SSI levels for the M 6.85 earthquake scenario

75 100

Jo

ur

na

lP

re

322

of

50

μ δ μ δ μ δ μ δ

Restoration Time (days) Conditioned on Source of Uncertainty Ground Shaking Component Damage Crew Size Process Duration 4.9 4.8 4.7 4.9 0.20 0.07 0.15 0.00 7.9 6.9 8.4 11.3 0.33 0.11 0.31 0.05 10.8 10.8 38.6 12.5 0.55 0.15 1.37 0.06 102.7 108.7 138.0 115.3 0.50 0.33 1.05 0.88

ro

25

Statistic

-p

Percentage of Pre-Event SSI

323 324

Figure 13 Mean and one-standard deviation bounded SSI curve for the M 6.85 earthquake scenario while

325

considering the uncertainty in the (a) ground shaking intensity, (b) component damage fragility, (c) crew

326

sizes and (d) event duration parameters

327

3.7 Risk-Based Performance Assessment

328

Figure 14, Figure 15 and Figure 16 show the annual exceedance rates for 22O;< , =P

329

respectively, which can be used to evaluate time-based performance targets i.e. probability of exceedance

330

over some duration. An example of a time-based target is that a desired performance level is maintained

331

within a duration ] (e.g. 100 years). By assuming that the failure to meet this target is an event that

332

follows a Poisson distribution, the associated probability is computed using Equation 3: bc

22OA< ,

(3)

of

0 = ^ ] = 1 _ . `a

%QQ; ,

A more specific hypothetical performance target is that the immediate post-event loss of SSI does not

334

exceed 50% within a 100-year period. Using the B Ĉ in Figure 14 that corresponds to 22O;< = 50%, the

335

probability of not meeting that target is computed to be 6%. Similarly, there may be a desire to limit

336

=P

337

this target is computed to be 26%. Using the annual exceedance relationship in Figure 16, the probability

338

that the cumulative SSI loss exceeds 20% 22OA< ^ 20% over a 100-year duration is 1%.

-p

lP

re

to less than 3 months within a 100-year period. From Figure 15, the probability of not meeting

na

%QQ;

Jo

ur

339

ro

333

340 341

Figure 14 Annual exceedance curve for 22O;<

of Figure 15 Annual exceedance curve for =P

%QQ;

-p

343

ro

342

345 346

Jo

ur

na

lP

re

344

Figure 16 Annual exceedance curve for 22OA<

347

4. Ground Motion Intensity and Damage Map Selection using Multi-Metric Optimization

348

Event-based methods are commonly used to assess the seismic risk to distributed infrastructure systems

349

because the relationship between the earthquake hazard and system performance metrics cannot be

350

expressed analytically. To reduce the computational expense associated with event-based methods, Miller

351

and Baker [8] proposed an optimization-based approach to selecting a subset of damage and

352

corresponding shake maps that provides a reasonable estimate of the full probability distribution of a

353

single network performance metric. However, as noted earlier, there is often a desire to characterize

performance of a distributed system using multiple metrics. This section discusses an extension of the

355

Miller and Baker optimization-based map selection algorithm to include multiple network performance

356

metrics in the objective function.

357

4.1 Map Selection Problem

358

The underlying goal of the map-selection problem is to choose a subset of d damage maps and

359

corresponding ground motion fields, each with an adjusted annual rate of occurrence , ′ . From an

360

optimization perspective, the difference between the network performance (as quantified by the threshold

361

exceedance curves) obtained from the subset d and the complete set of

362

minimized. More specifically, the multiple performance metric optimization problem can be expressed as:

ro

minimize

damage maps, should be

re

n

f E‖5 hi jk

`P

jk _ lk m ‖P

(4)

lP

oP

subject to

na

364

e

-p

363

of

354

ur

‖m‖P ≤ , q

(4a)

rs

\ ≤m

(4b)

where each element of m ∈ ℝJ × P , , ′ , represents the adjusted annual occurrence rate corresponding to

366

the + ′ damage map, ‖∎‖P is the L1-norm or sum of absolute values of the vector ∎, , q

367

the original occurrence rates of all

damage maps, \ ∈ ℝJ × P is a vector of zeros, jk is a vector of

368

annualized exceedance rates for the

network performance metric corresponding to

369

period values. lk ∈ ℝu ×J is a matrix of constants with each entry 9

370

value that controls the weight assigned to the

371

considered metrics.

372

The linear objective function in Equation 4 minimizes the difference between the exceedance rate curves

373

for the baseline set and the selected subset while considering all network performance metrics. The first

Jo

365

'

v,'

′

rs

is the sum of

discretized return

= F{C ′ ≥ Ĉv }, f is a scalar '

network performance metric and y is the total number of

374

constraint in Equation 4a specifies that the sum of , ′ should be less than or equal to , q

375

words, the sum of occurrence rates for the selected subset of maps should be less than or equal to that of

376

all maps in the baseline set. The second constraint in Equation 4b enforces non-negativity to every value

377

of , ′ to satisfy the physical interpretation of annual occurrence rates. The objective function and

378

constraints are linear in terms of m, therefore an optimal solution is guaranteed to exist [37]. CVXPY [38],

379

a python-embedded modeling language for convex optimization is used to solve the linear programming

380

problem. The map subset obtained after solving the optimization problem may result in more than d non-

381

zero values. The heuristic method [8,39] is used to select the d largest weights, which are renormalized

382

such that the sum of normalized weights is equal to , q

383

To evaluate the error in the network performance metric exceedance curves between the subset and

384

baseline set, the Mean Performance Measure Curve Error (MPMCE) shown in Equation 5 is used [8].

'

rs .

In other

lP

re

-p

rs .

ro

of

'

0 z{ =

P

u

∑uvoP }

b~ `b•€~ b~

}

(5)

where C•v is the estimated performance metric value using the selected subset of maps for the /

386

period.

387

4.2 Map Selection Results

388

The 2,343,000 damage maps for the Napa water network is used as the baseline and a subset of 1000

389

maps is selected using multi-objective optimization with equal weights assigned to all performance

390

metrics. Table 6 tabulates the MPMCE values obtained from optimizing on a single and multiple

391

performance metrics while testing on the remaining others. The values suggest that if the subset is

392

optimized for a single metric, for example, 22O;< , the MPMCE value is minimum for the metric used as

393

the basis of the optimization and yields high MPMCE values for the others. On the other hand, if the

394

subset is chosen while considering all network performance metrics, a universally high accuracy reflected

395

by low MPMCE values is obtained.

return

Jo

ur

na

385

396

Table 6: Mean performance measure curve error (MPMCE) for various target performance measures Optimized For All

VVWWX

Z[\\%VVW

1.31

1.31

1.32

1.31

Z[\\%VVW

3.47

3.54

3.48

3.61

VVWYX

3.47

3.47

3.48

3.47

Tested On

MPMCE (%)

VVWWX

397

VVWYX

5. Conclusion

399

A framework for evaluating the post-earthquake functional loss and restoration of a water distribution

400

system using a regional stochastic event set catalogue is presented. An end-to-end simulation is

401

implemented for each scenario, which begins with the characterization of ground shaking seismic hazard

402

and ends with a risk-based quantification of multiple resilience-related network performance metrics.

403

The methodology was applied to the city of Napa’s water distribution system. Of the 2,343 considered

404

scenarios, only 135 with magnitudes ranging from 6.25 to 8.0 caused damage to the Napa water network.

405

By coupling process-based discrete-event and pressure-driven hydraulic simulation (PBDES) models, the

406

complete functional recovery trajectory in terms of the system serviceability index (SSI) was modeled.

407

This formed the basis of the network performance assessment using metrics related to the initial loss of

408

SSI ( 22O;< ), the time to restore the pre-event SSI =P

409

Approximately 2.5% of the damage maps from the entire stochastic catalogue resulted in post-earthquake

410

water service disruption. Additionally, 95% of the damage maps achieved full SSI restoration within 3

411

days i.e. =P

412

Disaggregation of the dispersion in the network performance outcomes was used compare the relative

413

contribution of different sources of uncertainty including ground shaking intensity, component fragility,

414

crew size and event duration. When all sources of uncertainty were considered, the dispersion associated

415

with recovery time metrics (e.g. =P

416

results from the disaggregation showed that the ground shaking intensity had the highest contribution to

Jo

ur

na

lP

re

-p

ro

of

398

%QQ;

%QQ;

and the cumulative SSI loss (22OA< ).

≤ 3 days.

%QQ; )

was found to be higher than that of 22O;< and 22OA< . The

the dispersion in the restoration time corresponding to lower SSI levels (e.g. time to restore 25% and 50%

418

of pre-event level). Whereas at higher SSI levels (75% and 100% of pre-event level), the number of crew

419

members dominated the dispersion in the restoration time. In general, the dispersion in restoration time

420

increases with the SSI level.

421

Annual exceedance curves for the network performance metrics were used to evaluate time-based targets,

422

which can be directly used by decision makers for resilience planning. The Miller and Baker [8]

423

optimization-based map selection algorithm was also extended to consider multiple network performance

424

metrics in the objective function. On aggregate (when considering all metrics individually), the set of

425

maps selected based on multiple metrics performed better than those selected based on individual metrics.

426

Future work can extend the proposed framework to incorporate multiple hazards, operational adaptations,

427

the interaction with other lifeline systems and additional network performance indicators (like centrality

428

measures). With regards to interdependent lifeline systems, more research is needed to understand and

429

propagate the uncertainty associated with different types of dependencies (e.g physical, functional,

430

information).

431

6. Data Availability

432

Some or all data, models, or code used during the study were provided by a third party.

ro

-p

re

lP

na ur

Jo

433

of

417

1. Napa Water Network – Provided by NWD

434

Direct requests for these materials may be made to the provider as indicated in the

435

Acknowledgements.

436

7. Acknowledgements

437

This work was partially supported by the Nuclear Regulatory Commission (NRC) Research Grant

438

No. NRC-HQ-60-17-G-0028. Any opinions, findings, and conclusions expressed in the material

439

are those of the authors and do not necessarily reflect the views of the NRC. The authors would

440

also like to acknowledge Charles Scawthorn (SPA Risk), Katherine Klise (Sandia National

441

Laboratories), Joy Eldredge (Napa Water Division) and Douglas DeMaster (Napa Water Division)

442

for their valuable input. This work used computational and storage services associated with the Hoffman2

443

Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology

444

Group.

445

8. References

446

[1]

Masoomi H, Burton H, Tomar A, Mosleh A. Simulation-Based Assessment of Postearthquake Functionality of Buildings with Disruptions to Cross-Dependent Utility Networks. J Struct Eng

448

2020;146:4020070. [2]

Tomar A, Burton H V., Mosleh A, Lee JY. Hindcasting the Functional Loss and Restoration of the

ro

449

of

447

Napa Water System Following the 2014 Earthquake using Discrete Event Simulation. ASCE J

451

Infrastruct Syst 2020. doi:10.1061/(ASCE)IS.1943-555X.0000574. [3]

Guidotti R, Chmielewski H, Unnikrishnan V, Gardoni P, McAllister T, van de Lindt J. Modeling

re

452

-p

450

the resilience of critical infrastructure: The role of network dependencies. Sustain Resilient

454

Infrastruct 2016;1:153–68.

456

Tabucchi T, Davidson R, Brink S. Simulation of post-earthquake water supply system restoration. Civ Eng Environ Syst 2010;27:263–79. doi:10.1080/10286600902862615.

[5]

Çaǧnan Z, Davidson RA. Discrete event simulation of the post-earthquake restoration process for

Jo

457

na

[4]

ur

455

lP

453

458

electric

459

doi:10.1504/IJRAM.2007.015298.

460

[6]

461 462

systems.

Int

J

Risk

Assess

Manag

2007;7:1138–56.

Masoomi H, van de Lindt JW. Restoration and functionality assessment of a community subjected to tornado hazard. Struct Infrastruct Eng 2018;14. doi:10.1080/15732479.2017.1354030.

[7]

463 464

power

Crowley H, Bommer JJ. Modelling seismic hazard in earthquake loss models with spatially distributed exposure. Bull Earthq Eng 2006;4:249–73.

[8]

Miller M, Baker J. Ground-motion intensity and damage map selection for probabilistic

465

infrastructure network risk assessment using optimization. Earthq Eng Struct Dyn 2015;44:1139–

466

56.

467

[9]

Adachi T, Ellingwood BR. Serviceability of earthquake-damaged water systems: Effects of

468

electrical power availability and power backup systems on system vulnerability. Reliab Eng Syst

469

Saf 2008;93:78–88.

470

[10]

471

Ouyang M. Review on modeling and simulation of interdependent critical infrastructure systems. Reliab Eng Syst Saf 2014;121:43–60.

472

[11]

Fishman GS. Discrete-Event Simulation: Modeling, Programming, and Applications 2001.

473

[12]

Hosseini S, Barker K, Ramirez-Marquez JE. A review of definitions and measures of system

474

476 477

of

Rokneddin K, Ghosh J, Dueñas-Osorio L, Padgett JE. Bridge retrofit prioritisation for ageing

ro

[13]

transportation networks subject to seismic hazards. Struct Infrastruct Eng 2013;9:1050–66. [14]

Basöz N, Kiremidjian AS. A bridge prioritization method based on transportation system

-p

475

resilience. Reliab Eng Syst Saf 2016;145:47–61.

performance using GIS. Proc. 6th US-Japan Work. Earthq. Disaster Prev. Lifeline Syst. Tsukuba

479

Sci. City, Japan, 1995, p. 437–49.

considering structural deterioration of bridges. Struct Infrastruct Eng 2011;7:509–21. [16]

483 484

Earthq Eng Struct Dyn 2007;36:285–306. [17]

485 486

[18]

[19]

493

Prasad TD, Park N-S. Multiobjective genetic algorithms for design of water distribution networks. J Water Resour Plan Manag 2004;130:73–82.

[20]

491 492

Park JI, Lambert JH, Haimes YY. Hydraulic power capacity of water distribution networks in uncertain conditions of deterioration. Water Resour Res 1998;34:3605–14.

489 490

Winkler J, Duenas-Osorio L, Stein R, Subramanian D. Performance assessment of topologically diverse power systems subjected to hurricane events. Reliab Eng Syst Saf 2010;95:323–36.

487 488

Dueñas-Osorio L, Craig JI, Goodno BJ. Seismic response of critical interdependent networks.

ur

482

lP

481

Lee Y-J, Song J, Gardoni P, Lim H-W. Post-hazard flow capacity of bridge transportation network

na

[15]

Jo

480

re

478

Fragiadakis M, Christodoulou SE. Seismic reliability assessment of urban water networks. Earthq Eng Struct Dyn 2014;43:357–74.

[21]

Kessler A, Ormsbee L, Shamir U. A methodology for least-cost design of invulnerable water distribution networks. Civ Eng Syst 1990;7:20–8.

494

[22]

Napa Water Division. Urban Water Management Plan: 2015 Update. 2017.

495

[23]

Field EH, Dawson TE, Felzer KR, Frankel AD, Gupta V, Jordan TH, et al. Uniform California

496 497

earthquake rupture forecast, version 2 (UCERF 2). Bull Seismol Soc Am 2009;99:2053–107. [24]

Field EH, Arrowsmith RJ, Biasi GP, Bird P, Dawson TE, Felzer KR, et al. Uniform California

498

earthquake rupture forecast, version 3 (UCERF3)—The time-independent model. Bull Seismol

499

Soc Am 2014;104:1122–80.

501 [26]

503

for seismic hazard analysis. Seismol Res Lett 2003;74:406–19. doi:10.1785/gssrl.74.4.406. [27]

505

crustal regions. Earthq Spectra 2014;30:1025–55.

[28] Boore DM, Stewart JP, Seyhan E, Atkinson GM. NGA-West2 equations for predicting PGA, PGV,

507

and 5% damped PSA for shallow crustal earthquakes. Earthq Spectra 2014;30:1057–85. [29]

Campbell KW, Bozorgnia Y. NGA-West2 ground motion model for the average horizontal

ur

508

Abrahamson NA, Silva WJ, Kamai R. Summary of the ASK14 ground motion relation for active

lP

506

Field EH, Jordan TH, Cornell CA. OpenSHA: A developing community-modeling environment

-p

504

ro

Based Seismic Hazard Model for Southern California. Pure Appl Geophys 2011.

re

502

Graves R, Jordan TH, Callaghan S, Deelman E, Field E, Juve G, et al. CyberShake: A Physics-

of

[25]

na

500

components of PGA, PGV, and 5% damped linear acceleration response spectra. Earthq Spectra

510

2014;30:1087–115.

511

[30]

512 513

Chiou BS-J, Youngs RR. Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra. Earthq Spectra 2014;30:1117–53.

[31]

514 515

Jo

509

Idriss IM. An NGA-West2 empirical model for estimating the horizontal spectral values generated by shallow crustal earthquakes. Earthq Spectra 2014;30:1155–77.

[32]

Rezaeian S, Petersen MD, Moschetti MP, Powers P, Harmsen SC, Frankel AD. Implementation of

516

NGA-West2 ground motion models in the 2014 U.S. National Seismic Hazard Maps. Earthq

517

Spectra 2014;30:1319–33. doi:10.1193/062913EQS177M.

518 519

[33]

Jayaram N, Srinivasan K. Performance-based optimal design and rehabilitation of water distribution networks using life cycle costing. Water Resour Res 2008;44.

520

[34]

FEMA FEMA. Hazus-MH 2.1 Technical Manual: Earthquake Model 2012.

521

[35]

Tabucchi THP. Modeling Post-earthquake Restoration of the Los Angeles Water Supply System.

522 523

2007 2007. [36]

ATC. Critical Assesment of Lifeline System Performance: Understanding Societal Needs in

524

Disaster Recovery. US Department of Commerce, National Institute of Standards and Technology;

525

2016. [37]

Boyd S, Boyd SP, Vandenberghe L. Convex optimization. Cambridge university press; 2004.

527

[38]

Diamond S, Boyd S. CVXPY: A Python-embedded modeling language for convex optimization. J

[39]

Joshi S, Boyd S. Sensor selection via convex optimization. IEEE Trans Signal Process

-p

530

Mach Learn Res 2016;17:2909–13.

2008;57:451–62.

re

529

ro

528

of

526

Jo

ur

na

lP

531

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Jo

ur

na

lP

re

-p

ro

of

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: