- Email: [email protected]
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:
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: