The seismic performance of a lifeline network depends on the joint damage state of all its components, which in turn is dependent on their spatial correlation. As a result, lifeline networks are typically analyzed using a scenario-based approach. In the most common method, conventional Monte Carlo simulation (MCS), a rupture, ground motion, component damage, and system performance are simulated in turn; the process is repeated many times; and the results are combined (Han and Davidson 2012). While MCS is effective and conceptually straightforward, it is computationally inefficient, if not infeasible, because it requires so many expensive simulations to capture the full probabilistic distribution of performance. Recent methods have been developed to instead identify a relatively small, computationally efficient set of ground motion hazard scenarios, each with an associated adjusted annual occurrence probability, so that when probabilistically combined they represent the hazard accurately and efficiently (e.g., Jayaram and Baker 2010, Han and Davidson 2012, Zolfaghari and Peyghaleh 2015, Manzour et al., 2016). Similarly, Brown et al. (2011Brown et al. ( , 2013)), Gearhart et al. (2011Gearhart et al. ( , 2014)), and Miller and Baker (2015) describe methods to develop a computationally efficient set of damage scenarios, each with an associated adjusted occurrence probability.
In this paper, we apply the Brown et al. (2013) method to develop a computationally efficient set of damage scenarios for the Los Angeles water transmission system. We chose this method because it has the benefit of including, within each damage scenario, all the information necessary to evaluate every possible combination of component retrofit or replacement strategies. Each damage scenario describes, for each pipe segment (or other component) in the network, both its damage state assuming it is in its current condition and its damage state if it were to be replaced with a pipe of a different resistance level. That makes it easy to use the results in a subsequent optimization to determine what pipes to retrofit/replace and how so as to meet some system-level objectives. Damage scenarios can be used for subsequent flow modeling and to improve the decisiveness for reducing consequences resulting from the infrastructure hazard exposure. The damage scenario also includes the corresponding multihazard scenario and an associated adjusted annual occurrence probability. Each multihazard scenario includes effects of multiple earthquake-related hazards, specifically ground motion, liquefaction and surface rupture. Having a set of multihazard scenarios as well allows use of that same set of multihazard scenarios for analyses of different infrastructure systems within the same geographic region so they are consistent. Unlike Brown et al. (2013), we apply the method for a water transmission system (not a highway network), consider ground motion and liquefaction (not just ground motion), and conduct the analysis for a full, computationally efficient set of multihazard scenarios developed using the Manzour et al. (2016) method.
Goal and assumptions. The goal of the damage scenario selection method is to compute a relatively small set of Q damage scenarios to be used for subsequent system-level risk analysis. Each damage scenario q defines the damage state d of each system component i for each possible resistance level r it might have. Each damage scenario also includes an associated adjusted occurrence probability Pq|j such that for multihazard scenario j, ∑ 𝑃 𝑞|𝑗 𝑞 = 1. When probabilistically combined with the adjusted occurrence probabilities, the damage distributions for each component conditional on the hazard, computed using the set of Q damage scenarios (which we will call the "reduced set" of scenarios) should match the "true" damage distributions.
The annual occurrence probability of the damage scenario q is PjPq|j, the product of the annual occurrence probability of the multihazard scenario j and the conditional probability of damage scenario q given multihazard scenario j. Each damage scenario also includes the associated multihazard scenario j that generated the damage. Finally, the method computes the errors between the "true" and reduced set marginal damage distributions for each component i and resistance level r, so the tradeoff between reduced error relative to the "true" damage probabilities and increased computational efficiency for subsequent analyses (lower Q) can be examined explicitly. Figure 1 provides a schematic summary of the type of output provided by the damage scenario development method for a set of J multihazard scenarios each of which has Qj damage scenarios associated with it such that ∑ 𝑄 𝑗 𝑗 ≤ 𝑄. Note that 'no earthquake' at bottom of Figure 1 is one damage scenario, and there may be slightly fewer than Q damage scenarios total due to rounding. Overview. The method involves three main steps: (1) generating a set of J multihazard scenarios, (2) determining Qj, the number of damage scenarios q for each multihazard scenario j; and (3) using the damage scenario optimization to generate the set of damage scenarios q and determine the adjusted occurrence probability of each, Pq|j so as to match the marginal damage function for each system component (e.g., pipe) i.
Step 1. The first step is to generate a relatively small, computationally efficient set of J multihazard scenarios, such that when combined with their associated adjusted annual occurrence probabilities, Pj, the hazard curves they produce match the "true" hazard curves as defined by probabilistic seismic hazard analysis or a conventional MCS. We used the method described in Soleimani et al. (2020) although another could be used.
Step 2. The next step is to determine how many damage scenarios Qj to generate for each multihazard scenario j. The simplest approach is to simulate the same number of damage scenarios for each multihazard scenario. Assuming a desire to end up with a total of Q damage scenarios, that means Q/J damage scenarios for each multihazard scenario. Alternatively, to be as efficient as possible, we can vary the number of damage scenarios Qj to generate for each multihazard scenario j so as to minimize the variance in the estimated system-wide damage exceedance probabilities. For multihazard scenarios with very high (or very low) intensities, the probability of exceeding a specified level of system damage does not vary as much across damage scenarios, so there is no need to sample many of them. Using a method analogous to the one used to determine the number of ground motion scenarios for every earthquake scenario derived in Han and Davidson (2012, Section 4.3.3), we use Eq. 1 to determine the number of damage scenarios, Qj, to generate for each multihazard scenario j so as to end up approximately a user-specified total of Q damage scenarios across all multihazard scenarios. To ensure each multihazard scenario includes at least one associated damage scenario, we required 𝑄 𝑗 ≥ 1 after rounding.
level r=1
In Eq. 1, Pj is the probability of multihazard scenario j, Dn is the number of damaged components in the network associated with return period n in multihazard scenario j, P(𝑑 𝑗 ≥ 𝐷 𝑛 ) is the probability of exceeding Dn in multihazard scenario j, and N is the number of return periods n at which the damage exceedance probabilities are computed (four in the case study).
To compute the Dn values, a conventional MCS is used to simulate a large number of hazard scenarios and a damage scenario for each. A probability of exceedance vs. number of damaged components in the system curve is generated, and the values of Dn read from it (Note that the hazard portion of the MCS may have been done as part of Step 1, as in the case study). To compute P(𝑑 𝑗 ≥ 𝐷 𝑛 ), for each multihazard scenario j, a large number of damage scenarios are simulated and then P(𝑑 𝑗 ≥ 𝐷 𝑛 ) is computed as the proportion for which the number of damaged components is at least 𝐷 𝑛 . If Eq. 1 is used (rather than assuming the same number of damage scenarios for each multihazard scenario), a damage model is required to compute 𝐷 𝑛 and P(𝑑 𝑗 ≥ 𝐷 𝑛 ). There are many available, including, for example, Lanzano et al. (2014), Toprak et al. (2017), andBagriacik et al. (2018). The same damage model should be used in Step 2 and Step 3. The choice depends on what form of system damage is required for the application and available data.
Step 3. For each multihazard scenario j, use the damage scenario optimization (see Damage scenario optimization section below) to generate the set of damage scenarios q and determine the adjusted occurrence probability of each, Pq|j, so as to match the component damage functions. Keep track of which multihazard scenario j goes with each damage scenario q. Results include (a) a set of approximately Q damage scenarios, each of which identifies the damage state d of each component i for each possible resistance level r it might have, its annual occurrence probability (PjPq|j), and the multihazard scenario that generated the damage (Fig. 1); (b) errors between the "true" and reduced set marginal damage distributions for each component i and resistance level r.
Step 3. It is almost the same as that in Brown et al. (2013), except that we omit the constraint that requires the reduced set damage correlation to equal the "true" damage correlation. In our application, because we (1) incorporate spatial correlation in developing multihazard scenarios, and (2) assume there is no correlation in the component damage conditional on hazard, we do not need to retain the constraint. The spatial correlation in multihazard scenarios is considered through ground motion and liquefaction (using soil type) intensities calculation. Note that it is possible that the latter assumption is not true due to similarities in the design and construction between pipe segments, but there is little data on what those correlations should be and it is likely less important than correlations due to hazard.
∑ 𝑃 𝑞|𝑗 𝑞 = 1
(5)
In the formulation, Constraint (3) defines the errors terms as the difference between 𝑚 𝑖𝑟𝑑 , the "true" probability that component i in resistance level r is in damage state d and the estimated probability component i in resistance level r is in damage state d based on the reduced set of Q damage scenarios. If the reduced set probability overestimates the "true" probability, 𝑒 𝑖𝑟𝑑 + will be positive and 𝑒 𝑖𝑟𝑑 -will be zero; if it underestimates, 𝑒 𝑖𝑟𝑑 -will be positive and 𝑒 𝑖𝑟𝑑 + will be zero. The objective function (Expression 2) minimizes the sum of those errors. Note that if desired, one could include weights in the objective function to ensure especially small errors for particular pipes, damage states, or resistance levels. Constraint (4) requires that for each damage scenario q, each component i in resistance level r must be in exactly one damage state. Constraint (5) represents the assumption that exactly one damage scenario ultimately occurs for each multihazard scenario, so the sum of the conditional damage scenario probabilities is one. By providing lower bound 𝛼 for 𝑃 𝑞|𝑗 in constraint (8), we aim to have more common damage scenarios.
The optimization is a nonlinear mixed integer program since Equation 3 includes the product of two decision variables, 𝑃 𝑞|𝑗 and 𝑤 𝑖𝑟𝑑𝑞 . We solve it using the iterative heuristic algorithm described in Gearhart et al. (2014). Briefly, it includes the following four steps: We define a neighboring solution to an existing solution as one for which all 𝑤 𝑖𝑟𝑑𝑞 are the same except two and those two are the result of moving a single pipe from one damage state to another under a single damage scenario. To compute the value of a neighboring solution requires only computation of the errors in the marginal distribution for the single pipe under this swap. The algorithm stops when there are no moves of this nature that improve the objective function. The method is implemented in Python 3.5 and the optimizations are solved by Gurobi 6.0, a commercial solver (www.gurobi.com), with all multihazard scenarios solved in parallel.
Inputs. The damage scenario development method was applied to the Los Angeles Department of Water and Power (LADWP) water transmission network, as it existed in 2007. The LADWP system in this analysis includes approximately 2,200 km of pipe serving 680,000 customers representing 4 million people in a service area of approximately 1,200 sq. km.
In applying the method to the LADWP system, we considered damage only to pipes in the transmission network for multihazard scenarios including ground motion and liquefaction. There are 57,581 pipe segments of 50m or less. One could consider pump stations or other network components if desired. We used Bagriacik et al. (2018, Table 8) damage model to compute the probability of damage to each pipe segment as a function of peak ground velocity (PGV), liquefaction resistance index (LRI) (Cubrinovski et al. 2011), pipe length, diameter, and pipe material (e.g., steel, cast iron). Since we have very large diameter pipes in our data set, we refitted the data like R4 model (considered PGV, LRI, pipe length, pipe diameter and pipe material variables) by considering diameter a binary rather than continuous variable, equal to one when pipe diameter is greater than or equal to 300 mm and zero otherwise. Consistent with the formulation of the damage model, two damage states d were considered, damaged or not. In Step 2, we assumed four values of return period n-72, 475, 975, 2475 years. The value of Dn, the number of damaged components in the network associated with return period n, was obtained by a conventional MCS over 100,000 years. To compute the exceedance probability 𝑃(𝑑 𝑗 ≥ 𝐷 𝑛 ) in Equation 1, 10,000 damage scenarios were simulated for each multihazard scenario. The analysis was run assuming J=351 and Q=750. Also, the lower bound 𝛼 for 𝑃 𝑞|𝑗 in constraint ( 8) is assumed to be 0.01.
In this application, we define each pipe resistance level r as a P(damage|hazard) vs. hazard curve, with hazard expressed in terms of return period of PGV. Specifically, we consider four resistance levels r-low, moderate, high, and current state (table 1). The low, moderate, and high resistance levels are defined by assuming a set of specified damage probabilities for each of four return periods (72, 475, 975, 2475 years). We assume, for example, that if designed to moderate resistance, a pipe segment should have a probability of damage of 0.00075, 0.001, 0.005, and 0.01 if it experiences the PGV with a 72-, 475-, 975-, or 2475-year return period, respectively. For the current state, the curve is different for each pipe, as determined by application of the damage model. To compute the 𝑚 𝑖𝑟𝑑 values for a specified multihazard scenario j, i.e., the "true" probability that component i in resistance level r is in damage state d when subjected to multihazard scenario j, we found the PGV for pipe i in multihazard scenario j, determined the return period it corresponds to, then used the curve for the appropriate resistance r to determine the probability of damage, 𝑚 𝑖𝑟𝑑 . This format for defining resistance levels is consistent with the performance-based design thinking as described in Davis (2019aDavis ( , 2019b)).
Step 1 in the case study, we used the 351 multihazard scenarios developed in Soleimani et al. (2020). The ground motion hazard was based on the 2014 National Seismic Hazard Map, in turn based on the Uniform California Earthquake Rupture Forecast v3 (UCERF3) (Petersen et al. 2011, Field et al. 2014, Powers 2015) and the liquefaction hazard was based on the LPI3 liquefaction model in Kongar et al. (2017) and liquefaction susceptibility maps from Southern California Gas (SCG 2001(SCG -2003)). The final reduced set included 710 damage scenarios derived from those 351 multihazard scenarios (less than Q=750 due to rounding of the Qj values). Among these 351 scenarios, 285 multihazard scenarios (81%) have one or two damage scenarios, 60 multihazard scenarios (17%) have three to ten damage scenarios, and 2% have more than 10 damage scenarios (up to 29). The total number of damaged pipes varies across damage scenarios. Figure 2 shows the distribution for the current resistance level. With 57,581 pipes in the system, 82% of damage scenarios include fewer than 2,500 damaged pipes; and only less than 1% have more than 12,000 damaged pipes. The adjusted annual occurrence probabilities, Pj, for the multihazard scenarios range from 6.07(10 -8 ) and 2.50(10 -3 ). The adjusted annual occurrence probabilities, Pq|j, for the damage scenarios conditional on the multihazard scenarios are 0.01 to 1.0; and the adjusted annual occurrence probabilities, Pq, for the damage scenarios range from 3.9(10 -8 ) and 1.4(10 -3 ). Figure 3 indicates the probability of exceedance for the damaged pipes over all damage scenarios, assuming current resistance level. As an example, Figure 4 shows maps associated with one example damage scenario associated with a M6.38 earthquake on the Mission Hill fault. This particular damage scenario has an adjusted annual occurrence probability of 𝑃 𝑞 = 1.28(10 -5 ), based on the product of the adjusted annual occurrence rate of the associated multihazard scenario is 𝑃 𝑗 = 2.24(10 -4 ), and the probability of the damage scenario given the multihazard scenario is 𝑃 𝑞|𝑗 = 0.057. Figures 4a and4b show maps of the rupture trace and ground motion (PGV), and liquefaction potential index (LPI), respectively. Figure 4c shows which pipes are damaged assuming the current resistance levels and Figure 4d those damaged assuming that all are replaced with pipes having a moderate resistance level. The damage scenario also includes information about which pipes would be damaged if they were replaced with pipes designed to low resistance level or high resistance level. That information is not shown due to space limitations. Comparing Figures 4c and4d shows the different damage patterns resulting from the different resistance levels. In particular, while many pipes are damaged in the northern part of the city when assuming the current resistance levels, they are not when we assume they are replaced by pipes with moderate resistance level (as defined in table 1). Note that while these two maps show damage patterns assuming all pipes have their current resistance levels or all pipes have moderate resistance level, the information for each pipe can be mixed and matched. That is, if, for example, one wanted to do an analysis to see how the system would perform if only the pipes in the northern part of the city were replaced with pipes with moderate resistance, but the rest of the pipes remained in their current resistance level, you would simply use a damage scenario created by selecting the damage state associated with r=moderate for the pipes in the northern part and selecting the damage state associated with r=current for the remaining pipes. To assess how well the set of 710 damage scenarios match the marginal damage distributions associated with each pipe for each resistance level, we can examine the error terms, 𝑒 𝑖𝑟𝑑 + and 𝑒 𝑖𝑟𝑑 -. The mean of the error terms across all damage scenarios q, pipes i, resistance levels r, and damage states d equals 1.39(10 -10 ), which is close to zero and suggesting no notable bias. To get a sense for the magnitude of the errors, 97% of the absolute values of the error terms are within 1% of the "true" value of 𝑚 𝑖𝑟𝑑 and 99.84% are within 5% of the "true value of 𝑚 𝑖𝑟𝑑 , which should be sufficiently small for most practical purposes. For a particular application, the analyst can always conduct a sensitivity analysis by varying the maximum allowable number of multihazard scenarios J and the maximum allowable number of damage scenarios Q in the reduced set, and examining how the errors change, and select the appropriate computation-error tradeoff for the application. The optimization requires less than 5 minutes to reveal the optimal solution for each damage scenario.
This paper describes development of a computationally efficient set of damage scenarios for the Los Angeles water transmission system that considers multiple earthquake hazards. Each damage scenario describes, for each pipe segment (or other component) in the network, both its damage state assuming it is in its current condition and its damage state if it were to be replaced with a pipe of a different resistance level. The damage scenario also includes the corresponding multihazard scenario and an adjusted annual occurrence probability so that when probabilistically combined the set of damage scenarios with their probabilities matches the probabilistic hazard and system functioning (e.g., percentage of demand satisfied) exceedance curves.
Once developed, one can then simulate system functioning for each damage scenario in the reduced set to obtain a fully probabilistic analysis of system functioning that is computationally efficient enough to allow sophisticated hydraulic analysis rather than a more computationally efficient but less accurate connectivity-based system functioning analysis. Since each damage scenario includes for every pipe, both its damage state assuming it is in its current condition and its damage state if it were to be replaced with a pipe of a different resistance level, the resulting set of damage scenarios can be used in a subsequent optimization to determine what pipes to retrofit/replace and how so as to meet some system-level objectives. Otherwise it would be extremely computationally intensive to do so. Possible future work includes incorporating the effects of surface fault ruptures in the damage modeling; doing sensitivity analysis to understand the tradeoff between including more multihazard scenarios, J, and more damage scenarios per multihazard scenario Q; and including damage to water system components other than pipes (e.g., pump stations).