earthquake records, which can be regarded as closest proxy for the Cascadia subduction events);

**Figure 2** illustrates a methodology for probabilistic building vulnerability assessment. It consists of five basic steps:

Salient features of the key components of the framework are explained in the following.

# Structural Model

Tesfamariam and Goda (2015) studied the effect of MSAS earthquake records on the loss assessment of a 4-story non-code conforming RC space frame structure. This study extends this investigation to archetypical structures with different story numbers reported in Liel and Deierlein (2008). The archetype structures are: 2-, 4-, 8-, and 12-story non-code conforming RC buildings; the structures were designed as a space frame, and all columns and beams were part of the lateral resisting system. The buildings were designed according to the 1967 UBC seismic provisions (ICBO, 1967). Beam and column elements have the same amount of over-strength; each element is 15% stronger than the code-minimum design level. The design is governed by strength and stiffness requirements, as the 1967 UBC had few requirements for special seismic design or ductile detailing.

Finite-element modeling of structures can be achieved using a fiber or lumped plasticity model. In the fiber model, the element cross section is discretized and corresponding non-linear material properties of the core concrete, cover concrete, and reinforcing bars are assigned. On the other hand, in the lumped plasticity model, non-linearity of the beam-column element is introduced at the two ends (hinges), which are connected by an elastic element. Advantages and disadvantages of each approach are summarized in **Table 6**. Haselton et al. (2008) indicated that the lumped plasticity model, equipped with adequate hysteretic models for plastic hinges, can simulate global collapse behavior well (note: they observed that the fiber model may be numerically unstable when the responses become highly non-linear).

**Figure 3A** shows a schematic of the 4-story building. It has a floor area of 38.1 m (125 ft) by 53.3 m (175 ft); columns are spaced at 7.6 m (25 ft), and story heights are 4.6 m (15 ft) and 4.0 m (13 ft) at the ground floor and higher floor levels, respectively. The nonductile RC models used in this paper are developed by Liel and Deierlein (2008). The models are based on a lumped plasticity approach in Open system for earthquake engineering simulation (OpenSees, McKenna et al., 2000). The lumped plasticity element models used to simulate plastic hinges in beam-column elements (**Figure 3B**) utilize a tri-linear non-linear spring model that is developed by Ibarra et al. (2005) and implemented in OpenSees by Altoontash (2004). **Figure 3B** shows the tri-linear backbone curve, coupled with the associated hysteretic rules, which is used to model the structures to post-peak response and nearcollapse response. The post-peak response enables modeling of the strain hardening behavior associated with concrete crushing, rebar buckling and fracture, and bond failure (Haselton et al.,

Table 6 | Advantages and disadvantages of fiber and lumped plasticity models (Haselton et al., 2008).

2008; Liel and Deierlein, 2008). Liel and Deierlein (2008) and Haselton et al. (2008) reported that the Ibarra et al. model was calibrated with data from 255 RC column test results. Details of the calibration process and building details are provided in Liel and Deierlein (2008) and Haselton et al. (2008); for brevity, they are not repeated here. *P*-Δ effects are modeled using a leaning column. The vibration periods for the first three modes for the 2-, 4-, 8-, and 12-story buildings are summarized in **Table 7**.

### Seismic Hazard for Victoria and Ground Motion Selection

Victoria is the provincial capital of BC and is located at the southern tip of Vancouver Island (**Figure 1**). Due to its geographical location, Victoria is affected by three types of earthquakes. The first type of the influential events is an earthquake at shallow depth in the crust; historically, the 1918 and 1946 earthquakes fall under this category. The other two types of the influential earthquakes are related to the movements of the Juan de Fuca Plate, Explorer Plate, Gorda Plate, and North American Plate in the Cascadia subduction zone. In the subducting slab, deep earthquakes occur (e.g., 2001 Nisqually earthquake), while at the plate interfaces, mega-thrust subduction earthquakes, as larger as *M*w9.0, occur (e.g., 1700 Cascadia earthquake, Hyndman and Rogers, 2010). It is important to recognize that the three types of dominant earthquakes in southwestern BC have distinct characteristics in terms of recurrence interval, earthquake magnitude, location, and depth and thus should be treated differently.

The key features of the critical earthquake scenarios for a given location can be evaluated quantitatively via probabilistic seismic hazard analysis. Atkinson and Goda (2011) conducted seismic hazard studies for southwestern BC, by incorporating recent advancements in seismology. Typical outputs from probabilistic seismic hazard analysis, which are essential for seismic performance assessment of buildings and infrastructure, are the UHS and seismic deaggregation. Currently, the UHS at 2% PE in

Tesfamariam and Goda Bivariate seismic evaluation framework

50 years (equivalent to the return period of 2500 years) is adopted as the basis for seismic design provisions for new construction in Canada. The seismic deaggregation identifies critical earthquake scenarios (for instance, in terms of magnitude, distance, and earthquake type) for a selected probability level. **Figure 4A** shows UHS for Victoria at 10, 5, and 2% PE in 50 years, where the site condition is set to site class C, which is represented by the average shear-wave velocity in the upper 30 m between 360 and 760 m/s. The three probability levels are relevant for assessing the seismic performance of structures in Canada. **Figure 5** shows the seismic deaggregation results for *T* = 1.0 and 2.0 s for 10, 5, and 2% PE in 50 years; the selected vibration periods correspond to the adopted seismic intensity measure (IM) for the 2-story building and the 4-, 8-, and 12-story buildings, respectively (**Table 7**). In **Figure 5**, relative contributions due to crustal, mega-thrust (Cascadia) interface, and deep inslab earthquakes are indicated. The seismic deaggregation results suggest that relative contributions due to the Cascadia subduction earthquakes increase with the probability level and the seismic hazard values for longer vibration periods are affected more significantly by the large subduction events. The variable characteristics of the dominant scenarios are important for seismic performance evaluations and thus should be taken into account in selecting ground motion records for non-linear dynamic analyses of structural models.

Careful record selection is of critical importance to produce unbiased estimates of seismic vulnerability. In particular, when record scaling is implemented to reach high seismic excitation levels, record selection needs to account for the spectral shape effects (Luco and Bazzurro, 2007). One practical method that is widely adopted for mitigating the record scaling bias is the CMS method (Baker, 2011). In the CMS-based record selection, the target response spectrum is modified based on dominant earthquake scenarios and relevant ground motion prediction equations at a selected performance level. Typically, the base target response spectrum for record selection is a UHS and is further modified based on the mean scenarios obtained from seismic deaggregation; several tens of ground motion records that match the modified target response spectrum (i.e., CMS) are selected as input motion. However, for the seismic environments in southwestern BC, it may be too simplistic to use a single target response spectrum for a given probability level because three dominant earthquakes with different characteristics are present (**Figure 5**). For this reason, in this study, the multiple CMS-based record selection method by Goda and Atkinson (2011) is adopted, which defines three different target spectra considering the different earthquake characteristics and ground motion prediction

Table 7 | First three fundamental periods of 2-, 4-, 8-, and 12-story buildings.

models for these earthquake types. Examples of the CMS for crustal, interface, and inslab earthquakes are shown in **Figures 4B,C**; **Figure 4B** is for the 2-story building, whereas **Figure 4C** is for the 4-, 8-, and 12-story buildings. It is noted that the CMS for the interface events have richer spectral content with respect to other two earthquake types because of larger earthquake magnitudes and longer propagation paths.

Another important aspect for record selection is to prepare a suitable ground motion dataset for the seismic environments of interest. For southwestern BC, the base ground motion dataset should contain records from large mega-thrust subduction events. Moreover, the record database should contain as-recorded MSAS sequence records. To achieve these requirements, a new composite database of real MSAS sequences is compiled by combining the database that was constructed based on the Next Generation Attenuation database (Goda and Taylor, 2012) and the new database for Japanese earthquakes from the K-NET, KiK-nt, and SK-net (Goda et al., 2015). It is noteworthy that the new Japanese database includes records from the 2011 Tohoku earthquake, which may be considered as appropriate surrogate for the Cascadia subduction events. The composite dataset consists of 606 real MSAS sequence records; 75 sequences are from the NGA database and 531 sequences are from the Japanese database (each sequence has two horizontal components). This database is the largest dataset for as-recorded MSAS sequences and is sufficient to select a suitable set of record sequences by taking into account various requirements, such as earthquake type, magnitude, distance, and site class.

#### Incremental Dynamic Analysis

Incremental dynamic analysis implements a series of non-linear dynamic analyses by scaling a set of input ground motions based on an adopted IM, and develops prediction equations of engineering demand parameters (EDP, e.g., MaxISDR and ResISDR) at different IM levels. The IM is the spectral acceleration at the fundamental period of a structure. For the different building story numbers, the maximum scaling required in IDA can vary. For the 2-story building, the spectral acceleration at 1.0 s is selected as IM (**Table 7**) and the scaling range in IDA is varied from 0.05 to 1.4 *g*. For the 4-, 8-, and 12-story buildings, the spectral acceleration at 2.0 s (i.e., IM) ranges from 0.05 to 0.7 *g*. In general, numerical instability is encountered when the inter-story drift ratio of the frames exceeds 0.10. The first occurrence of such large deformation responses is treated as "*collapse*" (Vamvatsikos and Cornell, 2002). In characterizing the inelastic demand, non-linear responses that are in "collapse" and "non-collapse" states are distinguished. The collapse results are modeled by collapse fragility curves (see Collapse Fragility Assessment), whereas the non-collapse results are represented by multivariate seismic demand models (see Coupla-Based Seismic Demand Modeling). Eventually, the overall performance of the building is assessed by integrating collapse results and non-collapse results in the Section "Seismic Performance Evaluation."

Incremental dynamic analysis is carried out for the 2-, 4-, 8-, and 12-story RC frames using the set of 50 MS records as well as a set of 50 MSAS sequences, which are selected based on the multiple CMS-based procedures. The IDA results for both MS records and MSAS sequences (i.e., EDP-IM plot) are shown in **Figures 6** and **7**; **Figure 6** is for MaxISDR, whereas **Figure 7** is for ResISDR. To present the uncertainty of the IDA results, 16th–84th percentile curves (corresponding to mean ± 1 SD), are included in the figures. The overall characteristics of the IDA curves for MaxISDR and ResISDR are different; the former increases gradually with the seismic intensity level, whereas the latter increases rapidly when the seismic intensity level reaches in the range of 0.2–0.3 *g* for the 2-story building and 0.15*–*0.20 *g* for the 4-, 8-, and 12-story buildings; similar observations are also noted in FEMA P-58-1 (2012). It is noteworthy that the uncertainty of ResISDR is much greater than that of MaxISDR,

Figure 4 | (A) Uniform hazard spectra, (B) conditional mean spectra for the anchor vibration period of 1.0 s (for the 2-story RC frame), and (C) conditional mean spectra for the anchor vibration period of 2.0 s (for the 4-, 8-, and 12-story RC frames).

Figure 5 | Seismic deaggregation results for Victoria: (A–C) *T* **=** 1.0 s for 10, 5, and 2% probability of exceedance in 50 years, and (D–F) *T* **=** 2.0 s for 10, 5, and 2% probability of exceedance in 50 years.

as noted by Ruiz-García and Miranda (2006). To appreciate the differences of the IDA curves for the buildings with different story numbers, the 50th, 16th, and 84th percentile curves for the 4-, 8-, and 12-story buildings are overlaid together in **Figure 8**, noting that the same IM is adopted for these buildings (thus the IDA results can be compared directly). The results shown in **Figure 8** indicate that for a given seismic excitation level, both MaxISDR and ResISDR decrease with the story number; therefore, for the considered non-ductile RC frames, the 4-story building is more vulnerable than the other taller buildings.

Moreover, from the EDP-IM plots, it can be observed that the impact of aftershock records is significant for the 2-story building (**Figures 7A** and **8A**), whereas such marked effects diminish with increase in story number (**Figures 7B–D**–**8B–D**). One of the main reasons for the pronounced influence of aftershock records on MaxISDR and ResISDR for the 2-story building is related to its fundamental period (≈1.0 s; **Table 7**) and the dominant spectral content of the aftershock records; generally, aftershock records have richer spectral content in the short vibration period range (Goda et al., 2015). For all cases, the impact of MSAS earthquake sequence is more significant for ResISDR as compared with MaxISDR. For instance, for the 4-story building (**Figures 7B** and **8B**), in terms of median, the consideration of MSAS sequences leads to 5–10% increase for MaxISDR and up to 100% increase for ResISDR with respect to the results for MS records.

#### Collapse Fragility Assessment

The collapse fragility can be represented by a lognormal cumulative distribution function (CDF):

$$P\_{\rm c} = \Phi\left(\frac{\ln\left(\mathbf{x} / \Theta\right)}{\beta}\right) \tag{1}$$

where *P*C is the probability that a ground motion with IM = *x* will cause the structure to collapse, Φ(•) is the standard normal CDF, θ is the median of the fragility function (the IM level with 50% probability of collapse), and β is the SD of lnIM (sometimes referred to as the dispersion parameter). **Figure 9** shows the collapse fragility results (raw data and fitted lognormal curve) for MS records and MSAS sequences. The estimated values of θ and β are also provided in the figure. The impact of aftershocks is pronounced for the 2-story building, where the median collapse capacity θ is reduced by 13% (i.e., the curve is shifted toward

Figure 6 | Incremental dynamic analysis results for MaxISDR by considering MS records and MSAS sequence records: (A) 2-story, (B) 4-story, (C) 8-story, and (D) 12-story.

left). On the other hand, the collapse fragility curves of the 4-, 8-, and 12-story buildings show no or slight differences. These results are consistent with the IDA curves shown in **Figure 6**. Furthermore, in **Figure 10**, the collapse fragility results for the 4-, 8-, and 12-story buildings are superimposed. The comparison shown in **Figure 10** indicates that the median collapse capacity θ as well as the dispersion β increases with the story number; the differences of the collapsed fragility curves are more pronounced at the greater seismic excitation levels.

# Coupla-Based Seismic Demand Modeling

MaxISDR and ResISDR are statistically dependent (Goda and Tesfamariam, 2015) and thus this should be taken into account when these EDPs are characterized. For the seismic demand modeling, first, marginal probability distributions of MaxISDR and ResISDR should be developed, and second, corresponding dependence needs to be characterized. The probabilistic modeling of MaxISDR and ResISDR is performed at individual IM levels using non-collapse MaxISDR and ResISDR data (note: the number of available data points for seismic demand modeling decreases with the IM level because more data fall into collapse states; **Figure 9**).

**Figure 11A** shows the scatter plot for the 4-story building by considering MS records at 5% PE in 50 years level. In the figure, marginal distributions of MaxISDR and ResISDR are plotted along the horizontal axis and vertical axis, respectively. Note that ResISDR has a heavy right tail. Goda and Tesfamariam (2015) considered six probability distributions, i.e., lognormal, Gumbel, Frechet, Weibull, gamma, and generalized Pareto, for marginal probability distribution modeling of MaxISDR and ResISDR. For MS records and MSAS sequences, Goda and Tesfamariam (2015) showed that the Frechet distribution (Eq. 2) and generalized Pareto distribution (Eq. 3) are suitable for MaxISDR and ResISDR, respectively. The probability density functions of the Frechet and the generalized Pareto models are given by:

$$f(\mathbf{x}) = \frac{\xi}{\sigma} \left( \frac{\mathbf{x} - \boldsymbol{\mu}}{\sigma} \right)^{-1 - \frac{\xi}{\sigma}} \exp\left[ - \left( \frac{\mathbf{x} - \boldsymbol{\mu}}{\sigma} \right)^{-\frac{\xi}{\sigma}} \right] \tag{2}$$

and,

$$f(\boldsymbol{\omega}) = \frac{1}{\sigma} \Big( 1 + \xi \frac{\boldsymbol{\omega} - \boldsymbol{\mu}}{\sigma} \Big)^{-(1/\xi + 1)} \tag{3}$$

where μ is the location parameter, and σ is the scale parameter, and ξ is the shape parameter. These marginal distributions are non-normal (in particular, ResISDR); in such cases, conventional multivariate normal (or lognormal) distribution modeling is not ideal, and a more elaborate approach is necessary.

The dependence of MaxISDR and ResISDR can be characterized by using elliptical copulas, such as normal and *t*, and Archimedean copulas, such as Gumbel, Frank, and Clayton (McNeil et al., 2005). The asymmetric Archimedean copula is a mixture of one of the Archimedean copulas and the independence copula; this copula class is useful in modeling data that exhibit uneven distribution of the data points along the upperleft-lower-right diagonal line in the transformed space. In the

context of joint probability distribution modeling of MaxISDR and ResISDR, the uneven distribution of the data is related to the physical relationship between MaxISDR and ResISDR (i.e., MaxISDR ≥ ResISDR; Goda and Tesfamariam, 2015). To model the observed dependence of MaxISDR and ResISDR (e.g., scatter plot shown in **Figure 11A**), parametric copula functions are fitted to empirical copula samples using the maximum likelihood method (McNeil et al., 2005). The copula fitting of MaxISDR and ResISDR at various IM levels suggests that overall, the Gumbel (or asymmetrical Gumbel) copula (Eq. 4) is suitable for the majority of the cases examined in this study.

$$C\_{\delta}(\mu\_1, \mu\_2) = \exp\left(-[(-\ln \mu\_1)^{\delta} + (-\ln \mu\_2)^{\delta}]^{1/\delta}\right), \delta > 1\tag{4}$$

where *u*1 and *u*2 are the uniform random variables, and δ is the model parameter.

The developed statistical seismic demand models of MaxISDR and ResISDR can be used for seismic performance evaluation of structures. For instance, considering the fitted dependence function for the 4-story building at 5% PE in 50 years, numerous copula samples are first generated; their marginal distributions are uniformly distributed with the specified dependence characteristics. Using the simulated copula samples and the fitted marginal distribution models for MaxISDR and ResISDR, pairs of MaxISDR and ResISDR samples can be obtained using the inverse transformation method. The results of 5,000,000 simulations are presented in **Figure 11B**. Indeed, similar figures can be generated for different building story numbers as well as seismic hazard levels, and can be used in the seismic performance evaluation.

#### Seismic Performance Evaluation

The collapse fragility curves and the joint probability model of non-collapse inelastic seismic demands outlined in the previous sections can now be used to carry out performance-based

the 4-story building at 5% probability of exceedance in 50 years.

evaluation of a building with the limit states provided in **Tables 3** and **4**. For example, for the NBCC normal importance buildings (*basic objective*), the acceptable limit states are LS and CP for 10 and 2% PE in 50 years, respectively. For the NBCC high importance buildings (*essential service objective*), the required limit states are more stringent and correspond to DC and LS for the same performance levels. For these cases, the corresponding values of [MaxISDR, ResISDR] are as follows (**Table 4**):

For the structural models that are considered in this study (which should meet the *basic objective*), IO and DC are not applicable to evaluate their seismic performances based on bivariate structural responses. This is because the structures, when subjected to expected ground motions at IO and DC hazard levels, are essentially linear-elastic and residual responses are very small (near zero). In other words, the seismic hazard levels corresponding to 50% PE in 30 or 50 years are mainly related to the serviceability limit state and are too low to cause significant non-linear responses. As our focus in this paper is upon the non-linear responses, LS and CP are mainly concerned and an intermediate seismic performance level between LS and CP, i.e., 5% PE in 50 years (corresponding to the return period of 1000 years), is introduced.

To illustrate the proposed seismic performance evaluation method, three performance levels, i.e., 10, 5, and 2% PE in 50 years, are considered with the limit states of [MaxISDR, ResISDR] = [2.0, 1.0%], [3.0, 1.5%], and [5.0, 2.0%], respectively. These demand levels are similar to those presented in **Table 4**. **Figure 12A** shows the scatter plots of MaxISDR and ResISDR (for non-collapse cases) at the three performance levels, noting that the collapse cases are dealt with collapse fragility curves (**Figure 9**). The corresponding limit states are indicated with red broken lines. By connecting the limit state thresholds at different performance levels (gray broken lines) and plotting the seismic demands in bivariate space (blue dots), the evolution of the seismic performance evaluation of the structure can be visualized, facilitating the better understanding of the seismic performance of the structure at multiple seismic excitation levels.

The overall performance of the building is assessed through unconditional probability of unsatisfactory seismic performance (*P*NS) (i.e., overall measure at a seismic performance level). The steps followed to compute *P*NS are outlined below, with the results shown in **Figure 12B** as an example. **Figure 12B** illustrates the calculations of the probabilities of exceedance and nonexceedance of the specified limit state thresholds for the 5% PE in 50 years performance level for non-collapse cases. First, from **Figure 12B**, four probabilities of exceedance and non-exceedance can be derived:

• the lower-left quadrant corresponds to the probability of joint non-exceedance of the MaxISDR and ResISDR limits, *P*NE,NE (=0.494),

The four probabilities are useful for assessing the causes of unsatisfactory seismic performance for non-collapse cases. A large value of *P*E,NE tends to indicate that the unsatisfactory seismic performance is due to MaxISDR, whereas a large value of *P*NE,E suggests that the structure may need to be demolished after the earthquake. It is noteworthy that *P*NE,NE, *P*E,NE, *P*NE,E, and *P*E,E are conditional probabilities upon non-collapse cases. Second, the collapse probability *P*C and the non-collapse probability *P*NC, i.e., *P*NC= (1 −*P*C), need to be evaluated for the given seismic intensity level using the corresponding collapse fragility curve (**Figure 9**). Finally, once the different probability values are obtained as outlined above, the value of *P*NS can be calculated by:

$$\begin{split} P\_{\text{NS}} &= P\_{\text{C}} + P\_{\text{NC}} \times \{ P\_{\text{NE,E}} + P\_{\text{E,NE}} + P\_{\text{E,E}} \} \\ &= P\_{\text{C}} + P\_{\text{NC}} \times \{ 1 - P\_{\text{NE,NE}} \} \end{split} \tag{5}$$

**Figure 13** shows 4 by 3 panels (i.e., four buildings and three performance levels) of the bivariate MaxISDR–ResISDR data/ performance limits for MS records; four conditional probabilities of exceedance and non-exceedance as well as collapse/ non-collapse probabilities are indicated in the figure, whereas **Figure 14** shows the same set of results for MSAS sequences. To facilitate the comparison of the calculated probabilities for different cases, values of *P*NE,NE, *P*E,NE, *P*NE,E, *P*E,E, *P*C, *P*NC, and *P*NS are summarized in **Tables 8** and **9** for MS records and MSAS sequences, respectively.

**Figures 13** and **14** show that MaxISDR and ResISDR become severer with the increase in the seismic performance level; this can be inspected from the scatter of the data points as well as the increase of the collapse probability. The 2- and 4-story buildings are more vulnerable, in comparison with the 8- and 12-story buildings. The collapse probabilities for the 2-story building are generally greater than those for the 4-story building; however, for the non-collapse cases, MaxISDR and ResISDR data are more widely distributed and consequently, conditional probabilities of unsatisfactory seismic performance (e.g., *P*E,NE, *P*NE,E, and *P*E,E) for the 4-story building are greater than those for the 2-story building. Overall, unconditional probabilities of unsatisfactory seismic performance for the 4-story building are greater than others (*P*NS in **Tables 8** and **9**). Note that the causes of unsatisfactory seismic performance vary depending on building story numbers and performance levels for the non-collapse cases. For example, for the 2-story building, unsatisfactory performance is mainly due to large residual seismic demands; in this case, the damaged building may be demolished. On the other hand, for the 4-story building (e.g., 5% PE in 50 years), the unsatisfactory performance is mainly due to excessive peak transient seismic demands. These

results indicate that different counter measures may need to be implemented for different buildings as their damage mechanisms may be different.

The comparison of the results shown in **Figures 13** and **14** as well as **Tables 8** and **9** suggests that the observations made for MS records are generally applicable to MSAS sequences. However, additional seismic demands due to major aftershocks have noticeable influence on both MaxISDR and ResISDR for the 2-story building (**Figures 6A** and **9A**). Consequently, counter measures against aftershock risks should be specific to building types (i.e., dynamic structural characteristics and susceptible failure mode).

Importantly, the calculated values of PNS listed in **Tables 8** and **9** indicate that for all four non-ductile buildings, their seismic capacities may be judged as satisfactory (because PNS is relatively low) at the LS performance level (i.e., return periods of 500–1000 years), whereas they fail to meet the CP performance level required by the current standards suggested by FEMA P-58-1 (2012). Therefore, for this class of non-ductile RC buildings, seismic retrofitting should be implemented to improve the seismic performance.

#### Discussion and Conclusion

The primary objective of the building design code was LS. In developed countries, this has been met through improved seismic design provisions. Seismic vulnerability of existing buildings remains to be a major concern because of the use of older design codes and/or poor construction practices at the time of design and construction. Most of these older buildings are currently operational and are required to be further assessed and upgraded to improve potential economic consequences due to seismic damage. An accurate assessment of potential impact of future destructive earthquakes is essential for effective disaster risk reduction. Probabilistic seismic risk analysis entails comprehensive understanding of ground shaking information, such as fault rupture process, wave propagation, and site effects, as well as vulnerability of structures, such as structural damage accumulation, seismic loss generation, and societal/economic impact (Cornell and Krawinkler, 2000). Through probabilistic calculus, it evaluates the potential damage and loss that a certain group of structures is likely to experience due to various seismic events (Tesfamariam and Goda, 2015).

The current state of the art for seismic performance assessment of buildings in North America is FEMA P-58-1 (2012). It has been developed based on generic ground motions that are applicable to the seismicity of California, which might not be compatible with the seismicity in Canada. Furthermore, the damage observed from the MSAS sequence of the 2011 *M*w6.3 Christchurch earthquake in New Zealand has highlighted the need for further study on the collapse risk of RC buildings in Canada (Elwood, 2013). The rigorous probabilistic seismic performance evaluation method can be used to aid in an informed decisionmaking by comparing performance metrics of alternative seismic risk mitigation measures quantitatively. Accurate representation

Figure 13 | Mainshock earthquake record – probabilities of exceedance and non-exceedance of the 2-, 4-, 8-, and 12-story building (rows) and 10, 5, and 2% PE in 50 years hazard levels (columns). Red circles are the IDA results.

Figure 14 | Mainshock–aftershock earthquake sequence records – probabilities of exceedance and non-exceedance of the 2-, 4-, 8-, and 12-story building (rows) and 10, 5, and 2% PE in 50 years hazard levels (columns). Red circles are the IDA results.

#### Table 8 | Collapse and non-collapse probabilities and probabilities of exceedance and non-exceedance of the different limit states for MS records.

Table 9 | Collapse and non-collapse probabilities and probabilities of exceedance and non-exceedance of the different limit states for MSAS sequence records.

of different limit states, robust ground motion selection, and multivariate inelastic seismic demands are vitally important in the assessment. In this paper, a robust seismic evaluation tool, within the performance-based earthquake engineering framework, is developed. Two EDPs, MaxISDR and ResISDR, are used to determine the severity of seismic damage and consequences. The joint probability distribution and dependency are modeled using the advanced copula technique. Following SEAOC (1995) and FEMA P-58-1 (2012), the two EDPs reaching different performance limit states are defined. Moreover, the aftershock ground motions are incorporated with the conventional seismic performance evaluation methodology, and furnished a better representation of the prevalent risk. The proposed evaluation tool can indeed be used for existing structures or design of new buildings.

The proposed framework was applied to 2-, 4-, 8-, and 12-story non-ductile RC buildings located in Victoria, BC, Canada. Considering regional seismicity in southwestern BC (i.e., shallow crustal earthquakes, off-shore mega-thrust interface earthquakes from the Cascadia subduction zone, and deep inslab earthquakes), 50 MS records and 50 MSAS sequence records were selected. Subsequently, IDA was performed and the computed MaxISDR and ResISDR data were used for developing collapse fragility curves and for developing probabilistic inelastic seismic demand models using copulas.

The general conclusions related to the aftershock effects are as follows:

The unconditional probability of unsatisfactory seismic performance *P*NS integrates the collapse and non-collapse limit states and thus can be used as an overall seismic performance measure of structures. The general conclusions related to the *P*NS results for the four non-ductile RC frames are as follows:

Finally, the proposed performance-based seismic screening criteria and methods can be used for Canadian buildings. The methodology can be extended to different building types and seismicity (e.g., Eastern Canada). The consideration of MSAS sequences as

#### References

seismic input was found to be important for the seismic risk assessment of low- to mid-rise buildings, and further investigations are warranted in the future. Furthermore, the aftershock effects should also be integrated in the design of Canadian buildings.

#### Acknowledgments

Ground motion data for Japanese earthquakes and worldwide crustal earthquakes were obtained from the K-NET/KiK-net/ SK-net databases at http://www.kyoshin.bosai.go.jp/ and http:// www.sknet.eri.u-tokyo.ac.jp/, and the PEER-NGA database at http://peer.berkeley.edu/nga/index.html, respectively. This work was supported by the Natural Science Engineering Research Council Canada (RGPIN-2014-05013) to the first author and the Engineering and Physical Sciences Research Council (EP/ M001067/1) to the second author.

*Proceedings, Ninth Pacific Conference on Earthquake Engineering, Building an Earthquake-Resilient Society*, Auckland.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Tesfamariam and Goda. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Energy-Based Seismic Risk Evaluation of Tall Reinforced Concrete Building in Vancouver, BC, Canada, under** *M***w9 Megathrust Subduction Earthquakes and Aftershocks**

#### *Solomon Tesfamariam<sup>1</sup> \* and Katsuichiro Goda<sup>2</sup>*

*<sup>1</sup> School of Engineering, The University of British Columbia, Kelowna, BC, Canada, <sup>2</sup> Department of Civil Engineering, University of Bristol, Bristol, UK*

#### *Edited by:*

*Oren Lavan, Technion – Israel Institute of Technology, Israel*

#### *Reviewed by:*

*Marie-José Nollet, École de technologie supérieure, Canada Rita Bento, Universidade de Lisboa, Portugal*

*\*Correspondence: Solomon Tesfamariam<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

> *Received: 07 March 2017 Accepted: 28 April 2017 Published: 30 May 2017*

#### *Citation:*

*Tesfamariam S and Goda K (2017) Energy-Based Seismic Risk Evaluation of Tall Reinforced Concrete Building in Vancouver, BC, Canada, under Mw9 Megathrust Subduction Earthquakes and Aftershocks. Front. Built Environ. 3:29. doi: 10.3389/fbuil.2017.00029* This article presents a seismic performance evaluation framework for reinforced concrete (RC) buildings, comprising shear walls and gravity frames. The evaluation is undertaken within a performance-based earthquake engineering framework by considering regional seismicity and site-specific ground motion selection. Different engineering demand parameters (EDPs), i.e., maximum interstory drift ratio (MaxISDR) and energy-based damage index, are considered as performance indicators. Various prediction models of EDPs are developed by considering four ground motion intensity measures (IMs), i.e., spectral acceleration at the fundamental period, Arias intensity, cumulative absolute velocity (CAV), and significant duration of ground motion. For this study, a 15-story RC building, located in Vancouver, BC, Canada, is considered as a case study. By using 50 mainshock and 50 mainshock–aftershock (MS-AS) earthquake records (2 horizontal components per record and bidirectional loading), non-linear dynamic analyses are performed. Subsequently, the calculated MaxISDRs and damage indices are correlated with suitable IMs using cloud analysis, and the most *efficient* IM-EDP prediction models are selected by comparing standard deviations (SDs) of the regression errors. The MaxISDR of the shear walls is less than 1% for the mainshock and MS-AS records. The energy-based damage index shows sensitivity to delineate impact of earthquake types and aftershocks. The CAV is showed to be the most efficient IM for the energy-based damage index.

**Keywords: seismic performance, energy-based damage index, gravity frame, shear wall, reinforced concrete building, mainshock–aftershock earthquake**

# **INTRODUCTION**

#### **Motivation**

Seismic performance of reinforced concrete (RC) shear wall systems designed with Canadian design codes has been investigated by various researchers (e.g., Tremblay et al., 2001; Adebar et al., 2010; Boivin and Paultre, 2010, 2012; Luu et al., 2014). For RC core buildings designed with the CSA standard A23.3-04, Boivin and Paultre (2010) showed that the RC core performs satisfactorily for flexural demand, while potential deficiency under significant shear demand may be a concern. Koduru and Haukaas (2010) studied the seismic performance and economic loss of a 15-story RC building constructed in 1988 and located in Vancouver, BC, Canada. Their study was comprehensive, covering from regional seismic hazard, seismic vulnerability assessment, and economic impact estimation. Nevertheless, important improvements can be made with regard to use of ground motion records that are applicable to megathrust interface records from the Cascadia subduction zone, noting that the records used by Koduru and Haukaas (2010) were calibrated based on shallow crustal records. On the other hand, adopting PEER's performancebased earthquake engineering (PBEE) framework (Cornell and Krawinkler, 2000), Yang et al. (2012) carried out a seismic loss assessment for a 42-story RC dual-system building, i.e., a centrally located core wall building with perimeter special momentresisting frames. With a design earthquake intensity level, the maximum interstory drift ratio (MaxISDR) calculated was less than 2%.

In older Canadian codes, shear wall buildings are the primary seismic force resisting systems, and detailing of gravity frames are often neglected (Adebar et al., 2010). The poor detailing associated with the gravity frames can be (Tesfamariam and Saatcioglu, 2008) inadequate lap splice length, lap splice located in a potential plastic hinge zone, poor detailing of transverse reinforcement anchorage, welded detailing, and lack of support to longitudinal bars. Gravity frames, however, located in the plastic hinge zone of the shear wall can experience excessive deformation and, if not detailed properly, can sustain severe damage (Adebar et al., 2010). This type of damage, for example, was reported in the 27 February 2010 Maule Chile earthquake (Naeim et al., 2011) and the 22 February 2011 Christchurch earthquake (Stirrat et al., 2014). Furthermore, older buildings in Canada lack consideration of large interface events in seismic design procedures. (Note: the potential risk due to the Cascadia subduction earthquakes was only recognized in late 1990s.) The problem is further compounded with the prevalence of mainshock–aftershock (MS-AS) earthquake sequences.

The structural analysis, using an appropriate structural model for the calculation of engineering demand parameter (EDP) and collapse capacity/probability, is an essential component of the seismic vulnerability evaluation. To assess the probability of attaining a specific structural response level conditioned on seismic excitation, incremental dynamic analysis (Vamvatsikos and Cornell, 2002) and cloud and stripe analyses (Jalayer et al., 2007) can be used. Different damage indices are used as a surrogate measure for EDP*s* that are categorized as follows (Williams and Sexsmith, 1995): (a) non-cumulative, (b) deformation-based cumulative, (c) energy-based cumulative, and (d) combined (non-cumulative and energy-based) indices. The most common approach to relate seismic demands to structural performance limits (i.e., capacity) is based on non-cumulative drift-based EDP, such as MaxISDR and residual interstory drift ratio. For RC shear wall buildings, however, the drift-based damage indicator may show satisfactory seismic resistance performance while underestimating overall damage in the plastic region due to cyclic loading. In this article, the cumulative energy-based damage index proposed by Mehanny and Deierlein (2000), which takes into account both peak amplitude and duration of non-linear responses of structural members in quantifying the structural damage, will be used. Alternatively, other damage indices that were proposed in the literature (e.g., Park and Ang, 1985; DiPasquale and Cakmak, 1989; Reinhorn and Valles, 1995) can be adopted.

#### **Objectives**

This article presents a seismic vulnerability evaluation of a 15 story RC shear wall structure located in Vancouver, BC, Canada. The RC shear wall building includes gravity frames, which are more vulnerable to low-amplitude repeated ground motions and is modeled in OpenSees finite element software by accounting for the non-linearity in the model. The shear walls for the tall building act as a cantilever beam, and the plastic hinge is formed at the base (assumed to be the first four stories), whereas the gravity columns within this plastic region are modeled with non-linear material elements. The seismic risk assessment that is carried out in this study accounts for:

# **PBEE FOR SEISMIC VULNERABILITY ASSESSMENT**

Disaster risk reduction against future earthquakes requires decision support tools for cost-effective risk mitigation options. For seismic risk assessment and design, a PBEE methodology, originally advocated by Cornell and Krawinkler (2000) and later extended by various researchers (e.g., Goulet et al., 2007), can be adopted. In PBEE, the mean annual rate of exceedance of earthquake impact expressed in terms of damage measures (DM) and ν(DM), is quantified, involving seismic hazard, structural, and damage analyses. Mathematically, ν(DM) can be expressed as:

$$\text{v(DM)} = \iiint G(\text{DM}|\text{EDP}) dG(\text{EDP}|\text{IM}) |d\lambda(\text{IM})|,\qquad \text{(1)}$$

where λ(IM) is the mean annual rate that a certain level of IM is exceeded, *G*(EDP|IM) is the complementary cumulative distribution function of EDP given IM, and *G*(DM|EDP) is the complementary cumulative distribution function that can be characterized through damage analysis by relating EDP to the physical extent of structural damage, represented by DM. The accuracy of the earthquake impact assessment depends on the available data and the choices of the relevant models and parameters.

## **Hazard Consideration and IMs**

In western BC, Canada, three dominant earthquake sources are present: crustal, inslab, and interface events (Hyndman and Rogers, 2010). The latter two types are originated from the Cascadia subduction zone (**Figure 1A**), where the oceanic Juan de Fuca Plate sinks beneath the continental North American Plate. When the stored strain along the fault is released, a megathrust subduction earthquake, similar to the 2010 Maule Chile and 2011 Tohoku Japan earthquakes, can happen. The structural damage potential and consequences due to these three earthquake types can be significantly different because of their ground motion characteristics, depending on buildings and infrastructure of interest. Generally, in comparison with crustal and inslab earthquakes, large interface ground motions, originated from the Cascadia subduction zone, have much longer duration (**Figure 1B**). The spectral content of the ground motion records for three earthquake types can differ significantly due to different earthquake source characteristics in terms of magnitude and distance (**Figure 1C**). For instance, the effects of long-duration ground motions on tall buildings have been highlighted for the 2011 Tohoku earthquake (Takewaki et al., 2011).

The three earthquake types have different input characteristics. Thus, besides selecting appropriate EDPs, selecting the corresponding IM is important in the overall risk assessment. As such, consideration of spectral acceleration at the fundamental period *S*a(*T*1), which is most commonly adopted in modern seismic hazard and risk studies, may not be the most suitable indicator of the energy input from the ground motion. In this article, Arias intensity (AI; Arias, 1970), CAV (Electrical Power Research Institute , EPRI), and significant duration of ground motion (D5–95%, Trifunac and Brady, 1975) are considered, and their correlations with structural damage are quantified. The definitions of AI and CAV are given as follows:

$$\text{AI} = \frac{\pi}{2g} \int\_0 \left[ a(t) \right]^2 \text{dt} \tag{2}$$

**FIGURE 1 | (A)** Regional seismicity in southwestern British Columbia, Canada. **(B)** Sample ground motion records for crustal, interface, and inslab events. **(C)** Five percent damped response spectra of the sample ground motion records [single horizontal component shown in panel **(B)**] for crustal, interface, and inslab events.

and

$$\text{CAV} = \int\_0^\infty |a(t)| \text{dt},\tag{3}$$

where *a*(*t*) is the accelerogram. D5–95% is defined as duration between times the AI of a ground motion record reaches 5 and 95% of its final value.

#### **Energy-Based Damage Index**

Energy-based damage indices are cumulative and are computed with consideration of hysteretic response (Gosain et al., 1977; Park and Ang, 1985; Kraetzig et al., 1989; Mehanny and Deierlein, 2000; see **Table 1**). Gosain et al. (1977) formulated a model to describe damage by using energy absorption normalized by yield force and displacement. Park and Ang (1985) proposed a widely used damage index *D*PA as a linear combination of deformation and absorbed energy under cyclic loading. The weighing factor used for the energy term in *D*PA is calibrated through experimental work. Kraetzig et al. (1989) developed an energybased damage index *D<sup>K</sup>* that accounts for the energy dissipated in primary half cycles (PHCs) and follower half cycles (FHCs) for positive and negative parts of the response. Mehanny and Deierlein (2000) extended the Kraetzig et al.'s damage index by introducing weights on the PHCs and FHCs and the positive and negative damage indices and associated it with extent of physical damage. In this article, following Koduru and Haukaas (2010), the Mehanny–Deierlein damage index *D*MD is used. The limit states for *D*MD are shown in **Table 2**.

To define relationships between different earthquake return periods and acceptable performance limit states, seismic performance matrices are often adopted (**Table 3**). It is highlighted that for frequent [50% probability of exceedance (PE) in 30 years], occasional (50% PE in 50 years), rare (10% PE in 50 years), and very rare (2% PE in 50 years) earthquake design levels, the corresponding design performance limit states are immediate occupancy (IO), life safety (LS), and collapse prevention (CP), respectively. For the energy-based earthquake damage evaluation, limit states indicated in **Table 3** can be used with the corresponding definitions of the damage index (**Table 2**).

#### **SEISMIC PERFORMANCE EVALUATION**

#### **Structural Model**

The reference structure considered in this study is a 15-story RC building constructed in 1988 and located in Vancouver, BC, Canada (Ventura et al., 2001; Koduru and Haukaas, 2010). The primary lateral load-resistant element of the building is shear walls (**Figure 2**). The building is fairly regular in plan, with minor

#### **TABLE 2 | Limit states for the Mehanny–Deierlein damage index**.

setbacks at the fourth and fourteenth story levels (**Figure 2A**). The mixed-use building has commercial occupancy at the first floor and residential occupancy at the remaining floors. The shear walls in the staircase and elevator shafts are concentrated at the central core and form the main lateral load-resisting system (**Figure 2B**).

The first story height varies from 2.7 to 4.7 m, and subsequent stories are 2.7 m. Mass and stiffness of the building are used as a base model, and further simplifications are considered in developing a numerical model (Koduru, 2008). For example, four levels of underground parking below grade were not considered in the model, and the foundation was considered to be fixed at base. The numerical model for this building was developed by Koduru and Haukaas (2010) in OpenSees (McKenna et al., 2000). Finite element modeling of the structure was implemented using a fiber element, where the element cross-sections are discretized, and corresponding non-linear material properties of the core concrete, cover concrete, and reinforcing bars were assigned. The structural model consists of three components, i.e., gravity support columns,

**TABLE 3 | Vision 2000 recommended seismic performance objectives for buildings (SEAOC, 1995)**.

*, basic objective—proposed National Building Code of Canada (NBCC) normal importance; , essential service objective—proposed NBCC high importance;* ♢*, safety critical objective—not proposed NBCC category; ×, unacceptable performance for new construction.*

shear walls, and header beams. The modeling assumptions made by Koduru (2008) are outlined as follows:

Modal analysis is carried out; the first three modal vibration periods and corresponding damping ratios are obtained as follows. The first mode corresponds to the sway motion in the short structural axis direction; its vibration period and damping ratio are 0.9 s and 3.0%, respectively. The second mode is related to the torsional motion, and its vibration period and damping ratio are 0.84 s and 3.0%, respectively. The third mode is the sway mode in the long structural axis direction, and the corresponding vibration period and damping ratio are 0.25 s and 5.7%, respectively. The calculated vibration periods are in agreement with the measured vibration periods of the building by Ventura et al. (2001), i.e., first mode (0.81 s, short structural axis), second mode (0.79 s, torsional), and third mode (0.69 s, long structural axis).

These dynamic characteristics are important for selecting ground motion records for use in non-linear dynamic analyses.

# **Seismic Hazard for Vancouver and Ground Motion Selection**

The development of seismic damage prediction models requires a series of non-linear dynamic analyses of a structural model subjected to a set of ground motion records, which reflect the regional seismic hazard of interest. In this article, record selection based on multiple conditional mean spectra (CMS) for different earthquake types is carried out by following the same procedures described in the studies by Tesfamariam et al. (2015) and Tesfamariam and Goda (2015). The target CMS are developed for crustal, interface, and inslab earthquakes, based on full probabilistic seismic hazard analysis (PSHA) results by Atkinson and Goda (2011). The site of interest is Vancouver, and its surface soil is classified as site class C (average top 30 m shear-wave velocity ranges from 360 to 760 m/s). **Figures 3A,B** show the uniform hazard spectrum and seismic disaggregation result, respectively, at the return period of 2,500 years. The seismic disaggregation is based on spectral acceleration at 0.9 s (i.e., same as the fundamental vibration period of the building model). To develop CMS for different earthquake types, mean record characteristics for individual earthquake types are obtained from the PSHA results. Three CMS are included in **Figure 3A**, illustrating different spectral characteristics for the three earthquake types. For the considered case (i.e., Vancouver, site class C, vibration period of 0.9 s, and return period of 2,500 years), crustal, interface, and inslab events contribute equally to the overall seismic hazard.

To select a set of suitable ground motion records that match with the constructed CMS, an extended data set of real MS-AS sequences is used, which was developed by combing the worldwide (NGA) database (Goda and Taylor, 2012) with the Japanese (K-KiK-SK) database (Goda et al., 2015). The number of available MS-AS sequences is 606; among them, there are 197 crustal earthquakes, 340 interface earthquakes, and 69 inslab earthquakes. The interface events are from the *M*w8.3 2003 Tokachi-oki earthquake and the *M*w9.0 2011 Tohoku earthquake (which have similar event characteristics as the expected Cascadia subduction earthquake).

A set of ground motion records is selected by comparing response spectra of candidate records with the target spectra (i.e., CMS). The total number of records is set to 50. (Note: each record has two horizontal components.) The contributions from the 3 earthquake types are equal; as a result, 17, 17, and 16 records are selected for the crustal, interface, and inslab earthquakes, respectively. Response spectral matching is conducted in a least squares sense by considering the geometric mean of the response spectra of two horizontal components. (Note: spectral matching is performed for mainshock records of MS-AS sequences.) The vibration period range for spectral matching is from 0.1 to 2.0 s, which is inclusive of major vibration periods of the tall building model. **Figures 3C–E** show the statistics of the response spectra of the selected ground motion records (i.e., median as well as 16th and 84th percentile curves). For comparison, **Figures 3C–E** include the target CMS as well as the CMS *±* 1 conditional standard deviation (SD) (Jayaram et al., 2011). The response spectra of the selected records and the target CMS are similar for the crustal

and interface records; for inslab records, the selected records contain richer short-period spectral content than the target spectra. Given the availability of ground motion records and the dataset size of ground motion records (i.e., 16–17 for each earthquake type), matching of the candidate response spectra with the target is judged as adequate. **Figure 3F** shows the response spectra of the unscaled mainshock ground motions that are selected based on the preceding method. It is noted that the mean spectral acceleration at 0.9 s of the 50 records is about 0.36 *g*, which corresponds to the return period of 1,300 years, while 8 of 50 records exceed the spectral acceleration at 0.9 s that corresponds to the return period of 2,500 years (i.e., 0.5 *g*).

**Figure 4A** shows the magnitude–distance distribution of the selected earthquake records; in the figure, record characteristics for mainshocks and major aftershocks (i.e., events having the second largest magnitude within individual sequences) are included. Finally, **Figures 4B,C** compare the D5–95%-AI plot and the D5–95%-CAV plot for different earthquake types (mainshocks only). **Figures 4B,C** indicate that the interface records are associated with longer duration and larger CAV values than the crustal and inslab records.

The selected ground motion records are used for seismic performance assessment of the tall building in Vancouver. The records reflect regional seismic hazard and dominant record characteristics. In particular, consideration of the 2011 Tohoku records is relevant to the seismic performance assessment in Vancouver because of the anticipated macro-level similarity between the 2011 Tohoku earthquake and possible Cascadia events. The relative contributions from the crustal, interface, and inslab events are equal (for the considered scenario), and thus, these records can also be used for evaluating the effects of ground motion records having different record characteristics (i.e., spectral content and duration) on non-linear seismic demand and earthquake damage potential. The records can be employed in cloud analysis to develop probabilistic seismic demand models. For this purpose, target spectral acceleration levels need to be defined.

# **Dynamic Analysis and Cloud Analysis**

The structural analysis is carried for the 50 MS and 50 MS-AS earthquake records discussed in the previous section. The simulations are carried out using unscaled ground motions in bidirectional horizontal excitations of the 3D model (shown in **Figure 2A**). Various structural responses are stored for postprocessing, including time history data for interstory drift ratios and floor accelerations at all 16 stories, and Mehanny–Deierlein damage indices in the plastic hinge zone. The bidirectional interstory drift ratios were combined through geometric mean, and the corresponding drift values are used in the subsequent analysis. **Figure 5** illustrates calculated time histories of structural responses subjected to the three crustal/interface/inslab ground motion records, which are the same as those shown in **Figure 1B**. **Figure 5A** shows results for the interstory drift ratio at the ninth story, while **Figure 5B** shows results for the Mehanny–Deierlein damage index for gravity column at the third story. The results of the third-story column is selected because the damage index of this structural element is in the middle of all other structural elements of the gravity frame system and is thus suitable to show

**FIGURE 4 | (A)** Magnitude–distance distribution of mainshocks and major aftershocks of the selected ground motion sequences. **(B)** Duration–Arias intensity distribution of the selected mainshock ground motion records. **(C)** Duration–cumulative absolute velocity distribution of the selected mainshock ground motion records.

**FIGURE 5 | Sample time histories of structural responses subjected to three ground motion records (crustal, interface, and inslab; same records shown in Figures 1B): (A) interstory drift ratio at the ninth story and (B) damage index for gravity column at the third story**.

the average trends. It can be inspected from **Figure 5A** that under the three input records, the structure behaves elastically as the MaxISDR is only 0.3–0.4%, and there is no significant residual drift after the earthquake sequences. On the other hand, accumulation of cumulative seismic demands to the structure can be observed in **Figure 5B**; in particular, a major aftershock following the interface mainshock record (middle panel in **Figure 5B**) shows noticeable increase of the damage index. The maximum values of *D*MD for the gravity column are still less than 0.3 for the three records; therefore, the structural damage to this element is considered to be negligible for these cases (see **Table 2**).

**Figures 6A,B** show storywise profiles of interstory drift ratios (left) and relationships between MaxISDR and spectral acceleration (right), *S*a(*T*1)-MaxISDR, for MS and MS-AS records, respectively. The results are obtained from the cloud analysis, and the responses due to different earthquake types are color coded in the figures. For both MS and MS-AS records, MaxISDR is less than 1% (minimal damage) and tends to take the largest values at the 8th to 10th stories. This is a result of higher mode effects on the response of the structure. The drift is minimal in the plastic hinge zone (the first four stories), as a combined effect of the shear walls and gravity columns. The *S*a(*T*1)-MaxISDR plots display that MaxISDR is well correlated with *S*a(*T*1). In the figure panels, fitted prediction models, having a form of log10EDP = *a* + *b*log10IM, where *a* and *b* are the regression coefficients, are also included. As a quantitative measure of *efficiency* (Luco and Cornell, 2007),

**FIGURE 6 | Peak drift demands from cloud analysis displaying storywise profiles of interstory drift ratios and maximum interstory drift ratio and spectral acceleration relationships: (A) MS records and (B) mainshock–aftershock records**.

SD of the regression residuals is indicated in the figure panels. It is noteworthy that although detailed results are not shown, similar regression analyses were performed for different IM variables, such as spectral accelerations at different vibration periods, AI, CAV, and D5–95%. It was found that for MaxISDR, *S*a(*T*1) is the most efficient; this conclusion may be due to the elastic responses of the shear wall core under the ground motion records considered.

Overall, the results shown in **Figure 6** indicate the effects of aftershocks are not significant in terms of interstory drift. It is noteworthy that the average increase due to major aftershocks is about 8%. However, this increase is mainly caused by only 6 sequences of 50; the majority of the aftershocks do not increase the overall interstory drift. This observation is in agreement with Goda and Taylor (2012) and Goda et al. (2015). Under the considered ground motions, the shear wall core structure remains to be mainly in the linear elastic range, and this indeed reiterates that the shear walls RC cores are not vulnerable (e.g., Yang et al., 2012). For this reason, the subsequent investigations focus upon the seismic damage evaluation of the plastic zone, gravity columns (first to third stories), header beams (second to fourth stories), and shear walls (first to fourth stories) based on the energy-based damage index; see **Figure 2**. This focus is justified because the members in the plastic zone area are susceptible to severe damage due to cyclic loading (Koduru and Haukaas, 2010).

#### **Efficient IM for Energy-Based Damage Index**

Effectiveness of each IM on the estimation of EDPs is assessed using the concept of *efficiency*. An efficient IM results in relatively small variability of EDP given IM (Luco and Cornell, 2007). This property can be quantified by the SD of the regression residuals of predicted EDP values for a given IM. In this section, the Mehanny–Deierlein damage indices *D*MD for the structural elements in the plastic zone are considered as EDP. More specifically, in total, 10 damage indices are focused upon; 3 are for the gravity columns (first/second/third story denoted as *D*MD-C1, *D*MD-C2, and *D*MD-C3, respectively), 3 are for the header beams (second/third/fourth story denoted as *D*MD-HB2, *D*MD-HB3, and *D*MD-HB4, respectively), and 4 are for the shear walls (first/second/third/fourth story denoted as *D*MD-SW1, *D*MD-SW2, *D*MD-SW3, and *D*MD-SW4, respectively). On the other hand, four IMs are considered to identify the efficient IM parameters: *S*a(*T*1), AI, CAV, and D5–95%.

To compute efficiency of each IM, the log-linear model, i.e., log10EDP = *a* + *b*log10IM, is fitted using a least squares method. **Table 4** summarizes the logarithmic SD of the regression residuals for all combinations of EDP and IM for both MS records and MS-AS records. To show the results graphically, scatter plots of the damage index (*D*MD-C3) versus four IMs for MS-AS records are presented in **Figures 7A–D**. The results shown in **Table 4** and **Figure 7** indicate that the SDs of the residuals are smallest for CAV, followed by AI, *S*a(*T*1), and D5–95%. The results are consistent for all EDPs and for MS/MS-AS records. It can be concluded that CAV is the most efficient IM for *D*MD. Moreover, the results shown in **Figures 7A–D** suggest that interface records tend to result in greater damage index values in comparison

**TABLE 4 | Efficiency measure (i.e., logarithmic SD of regression residuals) for different intensity measures and engineering demand parameters (EDPs)**.

with crustal and inslab records. This is because interface records are long-duration ground motions (**Figure 4C**), and thus, their cumulative damage potential is higher than other short-duration ground motions. The consideration of *D*MD as EDP facilitates the incorporation of cumulative damage modes into the seismic performance evaluation.

**Figure 7E** compares the scatter plots of CAV and *D*MD-C<sup>3</sup> for MS and MS-AS records, respectively. The results clearly show that the effects of major aftershocks for *D*MD-C<sup>3</sup> are significant, resulting in increased earthquake damage. Based on the two fitted curves, the average increase of *D*MD can be quantified as 54%, which is in sharp contrast with the increase of MaxISDR shown in **Figure 6**. It is important to emphasize the differences of the aftershock effects on *D*MD-C<sup>3</sup> and MaxISDR. For *D*MD-C3, the effects due to major aftershocks are noticeable for the majority of the cases, rather than a small fraction of the cases (which was applicable to MaxISDR). To demonstrate this clearly, a histogram of the ratios of *D*MD-C<sup>3</sup> between MS-AS records and MS records is shown in **Figure 7F**. The results highlight the widespread influence of the major aftershocks on the damage index.

#### **Combined Damage Index**

The abovementioned results clearly indicate that the CAV is the efficient IM for all *D*MD, and the effects of earthquake types (long-duration interface events versus other earthquake types) and aftershocks have major influence on the earthquake damage evaluation. To perform probabilistic seismic risk analysis of the system in the plastic zone of the tall building, a combined measure of earthquake damage needs to be defined. It is noteworthy that the damage index computed for each component of the system represents local damage.

By assigning relative importance or weight to each local damage index, a global damage index can be computed. Park et al. (1987) proposed a story damage index *D*story as:

$$D\_{\text{story}} = \frac{\sum D\_i E\_i}{\sum E\_i},\tag{4}$$

where *D<sup>i</sup>* is the local damage index at location *i* and *E<sup>i</sup>* is the corresponding energy absorbed at location *i*. The energy dissipated, however, is also incorporated in the damage computation. This index can potentially misrepresent the overall damage state (Williams and Sexsmith, 1995). A more general combination rule for computing a global damage index was proposed by Bracci et al. (1989). They proposed weighing each local damage index by importance weight *wi*, and corresponding *Dstory* is computed as follows:

$$D\_{\text{story}} = \frac{\sum \boldsymbol{w}\_{i} \boldsymbol{D}\_{i}^{(c+1)}}{\sum \boldsymbol{w}\_{i} \boldsymbol{D}\_{i}^{c}}.\tag{5}$$

Parameter *c* is used to give higher importance to the most severely damaged elements. Bracci et al. (1989) suggested *c* = 1 and an equal weight to structural elements within the same story level. The global damage index is computed for each structural element, i.e., gravity columns, header beams, shear walls, within the plastic region. For simplicity, in this article, in assigning the weight, the number of structural elements at each floor is taken into consideration. For example, at the first and fourth stories, where only two types of structural members are considered, an equal weight of 0.5 can be assigned to each. On the other hand, for the second and third stories, three structural member types are present, and the weight can be assigned as 1/3. Furthermore, the damage index computed for each floor level is extended for the global damage index over the first four stories of the plastic regions, by assigning an equal weight (=0.25) to each floor.

Alternatively, more uniform weighting schemes may be considered. For instance, *c* = 0 in Eq. 5 corresponds to the arithmetic mean. Another popular choice for a uniform combination rule is the geometric mean of all contributing elements. In the following, these combination rules, in addition to Eq. 5 with *c* = 1, as proposed by Bracci et al. (1989), will be considered as a part of epistemic uncertainty associated with the structural damage assessment.

**Figure 8A** compares the 10 local damage indices of the gravity frame system for MS records with the combined damage index based on Eq. 5 with *c* = 1. It can be observed that the coupling beam at the second story shows the lowest damage index. The shear walls and gravity columns show higher damage indices, whereas the fourth story shear walls exhibit the highest damage potential. **Figure 8B** compares three combined damage indices in the plastic hinge zone for MS records, i.e., the Bracci et al. combination rule (Eq. 5 with *c* = 1), arithmetic mean, and geometric mean. The consideration of the Bracci et al. combination rule leads to higher values of the combined damage index because more weight is given to severe damage cases. The average ratios of the arithmetic mean and the geometric mean with respect to the Bracci et al. case are 0.67 and 0.52, respectively. This comparison illustrates the importance of the combination rule for defining the global damage index based on multiple local damage indices.

Based on the combined damage index for the gravity frame system (Bracci et al.'s combination rule with Eq. 5 and *c* = 1), prediction models of *D*MD,*<sup>C</sup>* in terms of CAV are developed as:

$$
\log\_{10} D\_{\text{MD},C} = -3.877 + 1.055 \log\_{10} \text{CAV} \tag{6}
$$

for MS records, and

$$
\log\_{10} D\_{\text{MD},\text{C}} = -3.538 + 1.002 \log\_{10} \text{CAV} \tag{7}
$$

for MS-AS records. The SDs of the regression residuals βD'CAV for Eqs 6 and 7 are 0.121 and 0.121, respectively. This is considered as the base case in the subsequent analyses.

As mentioned above, the combination rule of different local damage indices into a global damage index is an influential source of uncertainty. To investigate the effects of this uncertainty on seismic risk assessment, more uniform combination rules, such as arithmetic mean and geometric mean, can also be considered. By redefining the combined damage index *D*MD,*<sup>C</sup>* as arithmetic mean of the 10 local damage indices, seismic demand perdition models for *D*MD,*<sup>C</sup>* can be obtained as:

$$\log\_{10} D\_{\text{MD,C}} = -4.052 + 1.055 \log\_{10} \text{CAV} \tag{8}$$

**FIGURE 8 | (A)** Comparison of local damage indices with the combined damage index in the plastic hinge zone for MS records. **(B)** Comparison of three combined damage indices in the plastic hinge zone for MS records.

for MS records, and

$$
\log\_{10} D\_{\text{MD},\text{C}} = -3.741 + 1.011 \log\_{10} \text{CAV} \tag{9}
$$

for MS-AS records. The SDs for Eqs 8 and 9 are 0.120 and 0.115, respectively. Moreover, by adopting the geometric mean combination rule, seismic demand prediction models for *D*MD,*<sup>C</sup>* can be obtained as:

$$
\log\_{10} D\_{\text{MD},C} = -4.159 + 1.032 \log\_{10} \text{CAV} \tag{10}
$$

for MS records, and

$$
\log\_{10} D\_{\text{MD},C} = -3.911 + 1.008 \log\_{10} \text{CAV} \tag{11}
$$

for MS-AS records. The SDs for Eqs 10 and 11 are 0.123 and 0.110, respectively.

**aftershock records based on the damage index thresholds listed in**

**Tables 2 and 3**. *D*MD,*<sup>C</sup>* is calculated based on Eq. 5 with *c* = 1.

#### **Fragility Curves for Energy-Based Damage Index**

The developed CAV-*D*MD,*<sup>C</sup>* models (e.g., Eqs 6 and 7) can be used to derive fragility curves. The distribution of seismic demand about its median is often assumed to follow a two-parameter lognormal probability distribution. After estimating dispersion βD'CAV of the demand about its median, the fragility, i.e., probability that *D*MD *> D*MD,*<sup>C</sup>* at a given CAV, can be computed as:

$$\begin{split} P(\text{DMD} > \text{C\_{\text{MD},C}} | \text{CAV} = \text{cav}) \\ = 1 - \Phi \left( \frac{\ln(\hat{C}\_{\text{MD},C}) - \ln(a \cdot \text{cav}^b)}{\mathfrak{B}\_{D|\text{CAV}}} \right), \end{split} \tag{12}$$

where *C*ˆ is the median structural capacity associated with the limit state. By taking the damage limit states shown in **Tables 2** and **3** and the combination rule for the global damage index based on Eq. 5 with *c* = 1 (i.e., Eqs 6 and 7), three fragility curves are developed for the LS threshold (i.e., *D*MD,*<sup>C</sup>* = 0.30), near collapse (NC) threshold (i.e., *D*MD,*<sup>C</sup>* = 0.60), and collapse (C) threshold (i.e., *D*MD,*<sup>C</sup>* = 0.95). The fragility curves derived using Eq. 12 for MS records and MS-AS records are compared in **Figure 9**. The comparison of the fragility curves indicates that the effects of major aftershocks can be significant.

It is noteworthy that Eq. 12 accounts for statistical uncertainty associated with seismic demand predictions only. On the other hand, there are other important uncertain elements in assessing seismic fragility, such as aleatory uncertainty of capacity *C* and epistemic modeling uncertainty, denoted by β*<sup>C</sup>* and β*M*, respectively (Ellingwood et al., 2007). To include these effects in evaluating the seismic fragility, the SD βD|CAV in Eq. 12 can be replaced by:

$$
\mathfrak{B} = \sqrt{\mathfrak{B}\_{p|\_{\text{CAV}}}^2 + \mathfrak{B}\_C^2 + \mathfrak{B}\_M^2}. \tag{13}
$$

β*<sup>M</sup>* is assumed to be 0.20, by considering that the modeling process yields an estimate of building frame response that, with 90% confidence, is within *±*30% of the actual value (Ellingwood et al.,

Tesfamariam and Goda Energy-Based Seismic Risk Evaluation

2007). β*<sup>C</sup>* is assumed to be 0.25 for IO and LS, following Celik and Ellingwood (2009). For CP, however, β*<sup>C</sup>* can be considered to be 0.17 for the four-story gravity frame and 0.13 for the shear walls. Note that the mentioned values of β*<sup>M</sup>* and β*<sup>C</sup>* are obtained from the literature and are applicable to structural models that were considered therein. Therefore, caution should be exercised before adopting these recommended values of logarithmic SDs.

# **Probabilistic Seismic Risk Analysis**

It is important to assess the seismic performance of the gravity frame system of the tall RC building by taking into account uncertainties associated with regional seismic hazards, ground motions, and seismic vulnerability. The PBEE-based risk analysis methodology, as formulated in Eq. 1, is suitable for such assessments. In carrying out such risk assessments for the tall building in Vancouver, the regional seismic hazard model by Atkinson and Goda (2011) can be used as a starting point. Changes to the ground motion models are necessary because the adopted IM for the gravity frame system is CAV, instead of *S*a(*T*1). In this study, two ground motion models for CAV are considered. A model by Campbell and Bozorgnia (2012) is applicable to crustal earthquakes because this model was developed by using strong motion data from worldwide crustal earthquakes compiled in the PEER-NGA database. For the interface and inslab earthquakes, a model by Foulser-Piggott and Goda (2015) can be used because it was developed based on the extensive strong motion database of Japanese earthquakes, including the 2011 Tohoku data. For the seismic vulnerability assessment, the CAV-*D*MD,*<sup>C</sup>* models for MS records and MS-AS records developed in the study (see Eqs 6 and 7) are adopted. The numerical evaluation of Eq. 1 is conducted using Monte Carlo simulations. More specifically, a synthetic earthquake catalog of regional seismicity having five million years is generated from the regional seismicity model of Atkinson and Goda (2011). For a given seismic event in the catalog, a value of CAV is simulated by taking into account earthquake types (i.e., crustal, interface, and inslab). Subsequently, a value of *D*MD,*<sup>C</sup>* is sampled from the developed prediction models (i.e., Eqs 6–11). This is repeated for all seismic events in the catalog. In the CAVbased seismic hazard analysis, the average shear-wave velocity in the uppermost 30 m is set to 555 m/s (i.e., site class C). Once all values of CAV and *D*MD,*<sup>C</sup>* are evaluated, annual maximum hazard and risk values can be extracted to develop a seismic hazard curve (i.e., exceedance probability curve) for CAV [i.e., λ(IM) in Eq. 1] as well as a seismic risk curve for *D*MD,*<sup>C</sup>* [i.e., ν(DM) in Eq. 1]. In addition to the hazard and risk curves, seismic disaggregation plots can be obtained through the postprocessing of the results.

**Figure 10** depicts the CAV-based seismic hazard results (i.e., hazard curve and corresponding disaggregation) at the return period of 2,500 years. The seismic hazard curve shown in **Figure 10A** indicates that the CAV values at the return periods of 500 and 2,500 years correspond to 1,230 and 2990 cm/s, respectively. The disaggregation plot at the 2,500-year return period level highlights significant contribution from the interface events (i.e., events having magnitudes greater than 8.0). The increase is in sharp contrast with the counterpart for *S*a(*T*1), shown in **Figure 3B**. This is because long-duration characteristics of the interface events result in greater values of CAV.

**Figure 11** shows the *D*MD,*C*-based seismic risk results. Two risk curves for MS records and MS-AS records are presented in **Figure 11A**, whereas the seismic risk disaggregation plots for MS records and MS-AS records at the 2,500-year return period are shown in **Figure 11B** and **Figure 11C**, respectively. The seismic demand prediction models used for obtaining the results shown in **Figures 11A–C** are Eqs 6 and 7. To relate the estimated damage potential to the limit states, four ranges of the damage limit states for *D*MD (i.e., **Tables 2** and **3**) are indicated along the upper boundary of the figure panel. The results shown in **Figures 11A–C** indicate that the influence of major aftershocks on the damage potential is significant, increasing the damage index values by approximately 40% for a given probability level. When the mainshock effects only are considered, the return periods that correspond to incipient of LS and NC damage states (i.e., *D*MD,*<sup>C</sup>* = 0.3 and 0.6, respectively) are 650 years and 2050 years. These return period levels are decreased to 350 years and 1,050 years (i.e., greater risks) when the aftershock effects are taken into account in addition to those due to the mainshocks. Comparison of the disaggregation plots for CAV and *D*MD,C suggests that they are very similar; these are because the differences of the earthquake types effectively capture the long-duration effects, as shown in **Figure 7C**.

Finally, the effects of the combination rule for defining the global damage index on seismic risk assessment are investigated. For this purpose, seismic risk assessments are carried out by considering three sets of seismic demand prediction models of *D*MD,*C*, i.e., Eqs 6 and 7 versus Eqs 8 and 9 versus Eqs 10 and 11. The results are shown in **Figures 11D,E**. The base case (i.e., Eqs 6 and 7) leads to greater damage index values, in comparison with the two other cases (i.e., arithmetic mean and geometric mean) because more weights are given to severely damaged structural elements. The differences of the seismic risk curves for the three cases are significant, highlighting the importance of capturing this uncertainty in seismic risk assessments.

# **CONCLUSION**

Seismic performance of an RC shear wall system designed with Canadian design codes has shown acceptable performance in terms of drift limits. Recent damaging earthquakes have highlighted that MS-AS earthquake records are important factors in the overall risk assessment. Furthermore, gravity columns, which are not seismically detailed and thus have exhibited severe damage, can potentially lead to localized collapse. However, driftbased limit states showed little sensitivity to MS-AS earthquake records and impact of earthquake types. In this article, an energybased damage index is considered to capture the effects of longduration earthquake ground motions. It is important to assess the seismic performance of the gravity frame system of the tall RC building by taking into account uncertainties associated with regional seismic hazards, ground motions, and seismic vulnerability. The PBEE-based risk analysis methodology is suitable for such assessments. Thus, in this article, seismic performance of the 15-story RC shear wall building located in Vancouver, BC, Canada, was investigated. For the seismicity of Vancouver, BC, Canada, scope of the study and conclusions are summarized below.

# **REFERENCES**

computed. The influence of major aftershocks on the damage potential was significant, increasing the damage index values by approximately 40% for a given probability level.

*•* The combination rule for defining a global damage index based on local damage indicators has major influence on seismic risk assessments. This kind of epistemic uncertainty should be taken into account.

# **AUTHOR CONTRIBUTIONS**

Both authors have contributed in this paper. ST performed the modeling, analysis, and computed the energy-based damage index. KG contributed in the ground motion selection, efficiency measure computation, and risk analysis.

# **ACKNOWLEDGMENTS**

Ground motion data for Japanese earthquakes and worldwide crustal earthquakes were obtained from the K-NET/KiK-net/SKnet databases at http://www.kyoshin.bosai.go.jp/ and http://www. sknet.eri.u-tokyo.ac.jp/, and the PEER-NGA database at http: //peer.berkeley.edu/nga/index.html, respectively. This work was supported by the Natural Science Engineering Research Council Canada (RGPIN-2014-05013) to the first author and the Engineering and Physical Sciences Research Council (EP/M001067/1) to the second author. Part of this research was undertaken when the first author was a Benjamin Meaker visiting professor at the University of Bristol and the financial support is acknowledged.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Tesfamariam and Goda. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Estimation of Seismic Loss for a Portfolio of Buildings under Bidirectional Horizontal Ground Motions due to a Scenario Cascadia Event**

*Taojun Liu1,2 \* and Hanping Hong<sup>3</sup>*

*<sup>1</sup>Department of Civil, Environmental and Architectural Engineering, University of Colorado Boulder, Boulder, CO, United States, <sup>2</sup>United States Geological Survey, Golden, CO, United States, <sup>3</sup>Department of Civil and Environmental Engineering, University of Western Ontario, London, ON, Canada*

Earthquake ground motions induced by a scenario event are spatially (partially) correlated and (partially) coherent. Simulated ground motion records can be used to carry out nonlinear inelastic time history analysis for a portfolio of buildings to estimate the seismic loss, which is advantageous as there is no need to develop and apply empirical ground motion prediction equations and the ductility demand rules, or to search the scenariocompatible recorded records at selected sites that may not exist. Further, if the structures being considered are sensitive to the orientation of the excitation, multiple-component ground motion records are needed. For the simulation of such ground motion records, previous studies have shown that correlation and coherency between any pair of ground motion components need to be incorporated. In this study, the seismic loss of a portfolio of hypothetical buildings in downtown Vancouver under bidirectional horizontal ground motions due to a scenario Cascadia event is estimated by using simulated bidirectional ground motion records that include realistic correlation and coherency characteristics. The hysteretic behaviors of the buildings are described by bidirectional Bouc–Wen model. The results show that the use of unidirectional ground motions and single-degreeof-freedom system structural model may underestimate the aggregated seismic loss.

#### *Edited by:*

*Katsuichiro Goda, University of Bristol, United Kingdom*

#### *Reviewed by:*

*Shinichi Matsushima, Kyoto University, Japan Alkis Daskaloudis, Mott MacDonald, United Kingdom*

> *\*Correspondence: Taojun Liu<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

*Received: 07 June 2017 Accepted: 02 October 2017 Published: 24 October 2017*

#### *Citation:*

*Liu T and Hong H (2017) Estimation of Seismic Loss for a Portfolio of Buildings under Bidirectional Horizontal Ground Motions due to a Scenario Cascadia Event. Front. Built Environ. 3:61. doi: 10.3389/fbuil.2017.00061* **Keywords: seismic risk, ground motion simulation, bidirectional excitation, 2-degree-of-freedom hysteretic model, Cascadia earthquake**

# **INTRODUCTION**

Seismic loss estimation for a portfolio of buildings under scenario events generally requires three sets of information. The first one is the scenario event and its associated multiple component ground motions at the spatially distributed sites of the buildings. The second set is associated with non-linear inelastic dynamic characteristics of the buildings and their responses or degree of damage under seismic excitations. The third set contains the damage loss functions for different structure types and degree of damage. For simplicity, seismic loss estimation for a portfolio of buildings is often carried out by using the ground motion measures such as the peak ground acceleration (PGA) and spectral acceleration (SA) for random orientation. The structural responses and damage are then represented using predetermined fragility curves based on experience, expert opinion, or numerical analysis and experimental investigation [HAZUS-Earthquake, Federal Emergency Management Agency (FEMA) and the National Institute of Building Sciences (NIBS), 2003; Whitman et al., 1997]. The seismic risk of a portfolio of buildings is then estimated by incorporating the (uncertainty in the) ground motion measures of the scenario event, the fragility curves, and the damage cost functions. Therefore, significant computing task in this approach, in terms of structural responses, is to establish the fragility curves for generic structures of different structural types. Instead of using fragility curves, Goda and Hong (2008a,b) considered that each building can be approximated as a nonlinear inelastic single-degree-of-freedom (SDOF) system, and the spatially distributed buildings are subjected to spatially correlated ground motion measures. To assess the degree of damage, use of the ductility demand rules for bilinear systems developed based on selected ground motions (Hong and Hong, 2007) were considered. It was shown that the consideration of realistic spatial correlation is crucial in assessing the tail of the probability distribution of seismic loss of a portfolio of buildings. The approach avoids the need for predetermining the fragility curve of the degree of damage conditioned on the ground motion measure such as SA, and approximately takes into account dynamic and inelastic characteristics of each of the buildings. However, the ductility demand rules could be affected by records from different earthquake types (Hong et al., 2010).

To avoid the need to develop empirical ductility demand rules for structures having different hysterical behavior and different ground motion characteristics, Liu and Hong (2015a) considered that the structural responses and damage levels can be estimated directly through the time history analysis under spatially correlated and coherent ground motion excitations. Their study, again, showed the importance of considering realistic spatial correlation in assessing the tail of the probability distribution of seismic loss of a portfolio of buildings. The approach of using time history analysis is also advantageous as there is no need to develop and apply empirical ground motion prediction equations. Furthermore, the use of the simulated ground motion records avoids the search for the scenario compatible actual records at the considered building sites that are unlikely to be available in the existing database of ground motion records.

For the simulation of the ground motion records, both the spatial correlation (Goda and Hong, 2008a; Jayaram and Baker, 2009; Liu et al., 2012) and spatial coherency (Abrahamson et al., 1991; Zerva, 2009) need to be taken into account. The coherency between two ground motion record components can be estimated from the power spectral density functions of the records; it represents the correlation between the random phase variations. The spatial correlation is used to measure the correlation of ground motion measures such as the PGA or SA at two sites. Furthermore, seismic events cause multidirectional ground motions. Methodology for the simulation of ground motion records that considers both spatial correlation and coherency for multidirectional excitations at multiple sites was presented in Hong and Liu (2014), Liu and Hong (2015a,b) based on stochastic simulations (e.g., point source model and finite-fault model) (Motazedian and Atkinson, 2005; Atkinson et al., 2009). These methods apply partially (and directional-dependent) coherent white noises generated using spectral representation method as the input for stochastic simulation techniques and incorporates spatially correlated Fourier amplitude spectra (FAS).

It must be emphasized that although multidirectional synthetic ground motion record components at multiple sites can be simulated, the estimation of seismic loss by considering multidirectional excitations for a portfolio of buildings for a scenario event has not been investigated. This can be important as the seismic response of buildings could be sensitive to multidirectional excitations (Clough and Penzien, 2003; Zerva, 2009), and a building could be modeled approximately by using a non-linear inelastic 2-degree-of-freedom (2DOF) system in lieu of a SDOF system, where each degree of freedom is associated with one of the two orthogonal horizontal directions. The hysteretic behavior of the 2DOF system could be modeled by the Bouc–Wen model (Wen, 1976; Lee and Hong, 2010), which can be used to reproduce sophisticated inelastic behavior of structural components/systems under cyclic loadings.

The main objectives of this study are to provide an overall framework to estimate seismic loss of a group of buildings under multidirectional excitations, and to investigate the effect of bidirectional ground motions on the aggregate seismic loss of buildings for a scenario event. The tasks include (a) simulating ground motion record components in two horizontal orthogonal directions at multiple sites for a scenario event such as that from the Cascadia subduction zone considering spatial (and directional) correlation and coherency models derived from historical records; (b) applying the simulated records in estimating seismic responses, and aggregate losses of a portfolio of buildings for the scenario event. In the following, first, the framework to estimate the aggregate seismic loss for a portfolio of buildings under multidirectional ground motions is described. The overall framework is then illustrated by a numerical example focused on the estimation of the seismic loss of a portfolio of buildings located in downtown Vancouver under a scenario Cascadia event, though we expect the conclusions of this study are independent of the study area.

# **FRAMEWORK FOR ESTIMATING SEISMIC LOSS OF A PORTFOLIO OF BUILDINGS CONSIDERING BIDIRECTIONAL HORIZONTAL GROUND MOTIONS**

In this section, the proposed framework to estimate the aggregate seismic loss of a group of buildings is presented. The procedure consists of three major components: the simulation of ground motion records (or field); the approximation of structural modeling and nonlinear inelastic analysis, and the estimation and characterization of aggregate seismic loss using the damage cost functions. For the simulation of ground motion records, approach and empirical correlation and coherency models among ground motion components given in the literature (Hong and Liu, 2014; Liu and Hong, 2015b) are considered, except that the reference FAS and time modulation function of the records are defined using the stochastic finite-fault method (Atkinson et al., 2009). For efficiency, the bidirectional Bouc–Wen model is used to calculate the non-linear inelastic response of the 2DOF system and the damage factor; the aggregated seismic loss for the portfolio of buildings is estimated by adopting cost functions in the literature.

# **Simulation of Bidirectional Ground Motion Records at Multiple Sites**

For the simulation of ground motion record components in two horizontal orthogonal directions at multiple sites, it is considered that the reference FAS and the time modulating functions for a random horizontal component can be defined based on the finitefault model (Motazedian and Atkinson, 2005; Atkinson et al., 2009). As the two horizontal orthogonal ground motion components at a site are considered to be random with respect to the source-to-site orientation, both the reference FAS and the time modulating function for one direction that is perpendicular to the other orientation is considered to be the identical.

To obtain the reference FAS and time modulating function at the *j*-th site, simulation by using the finite-fault model is carried out *n*<sup>G</sup> times for a considered scenario event with moment magnitude **M**. For each simulated record component, the FAS is evaluated and the time window profile (i.e., time modulating function) is estimated using the Hilbert transform. The reference FAS denoted as *yj*(**M**, *Rj*, *f*), and the time modulating function are then calculated by averaging over *n*<sup>G</sup> simulations, where *f* is the frequency in Hz and *R<sup>j</sup>* is the distance from the *j*-th site to finitefault source (i.e., closest distance to the fault plane). Based on the obtained reference *yj*(**M**, *Rj*, *f*), and the time modulating function, the simulation of the ground motion record components in two horizontal orthogonal directions are carried out as illustrated in **Figure 1** and outlined below (Hong and Liu, 2014):

The target coherency functions needed in Step (a) by considering the *j*-th and *k*-th sites is denoted by ¯γpq,jk(Δ*, f*), where *j*, *k* = 1, *. . .*, *n*<sup>R</sup> represents the sites, *p, q* = 1 or 2 represents the first

and second horizontal ground motion component (for a common coordinate system) and Δ (kilometers) is the distance between the *j*-th and *k*-th sites. By taking into account the symmetry and the fact that ¯γpp,jj (0*, f*) = 1 by definition, there are six remaining coherency functions need to be considered: ¯γ12*,*jj(0*, f*), ¯γ11*,*jk(Δ*, f*), ¯γ12*,*jk(Δ*, f*), ¯γ21,jk(Δ*, f*), ¯γ22*,*jk(Δ*, f*), and ¯γ12,kk(0*, f*). The coherency function for two ground motion components along the same orientation (i.e., ¯γ11*,*jk(Δ*, f*) or ¯γ22*,*jk(Δ*, f*)), can be expressed as (Harichandran and VanMarcke, 1986),

$$\bar{\gamma}\_{\rm pp,jk}(\Delta, f) = \left| \bar{\gamma}\_{\rm pp,jk}(\Delta, f) \right| \exp\left( -i2\pi f \Delta\_{\rm P} / \nu\_{\rm up} \right), \text{ for } j \neq k \quad \text{(1)}$$

where Δ<sup>P</sup> is the projection of the separation distance Δ in the direction of wave propagation; *vap* (kilometer per second) represents the apparent velocity; 2π*f* ΔP/*v*ap represents phase angle of the wave passage effect (Der Kiureghian, 1996); the lagged coherency ¯γpp,jk(Δ*, <sup>f</sup>*) is given by (Harichandran and VanMarcke, 1986),

$$\begin{split} \left| \bar{\gamma}\_{\text{pp,jk}} (\Delta, f) \right| &= A \exp \left( -\frac{2000 \Delta}{\alpha\_0 \Theta(f)} (1 - A + \alpha\_0 A) \right) \\ &+ (1 - A) \exp \left( -\frac{2000 \Delta}{\Theta(f)} (1 - A + \alpha\_0 A) \right), \end{split} \tag{2}$$

in which θ(*f*) = *k* ( 1 + (*f/f*0) *B* )*<sup>−</sup>*1*/*<sup>2</sup> , and *A*, α0, *k*, *f* <sup>0</sup>, and *B* are model parameters.

According to Hong and Liu (2014) and Liu and Hong (2015b), the lagged coherency for two horizontal orthogonal components can be considered to be independent of Δ, and can be approximated by

$$\left| \bar{\chi}\_{\text{pq,jk}}(\Delta, f) \right| = c\_0 - c\_l f, \text{ for } p \neq q \tag{3}$$

where *c*<sup>0</sup> and *c*<sup>1</sup> are model parameters. A set of typical parameters for the model shown in Eqs 1–3 are listed in **Table 1**. Note that these parameters are developed based on data from Taiwan. However, previous studies did not show dependency on geographical locations (Harichandran and VanMarcke, 1986; Hong and Liu, 2014; Liu and Hong, 2015a). Furthermore, to the authors knowledge, there is no literature discussed the impact of different geological and seismic settings on (spatial) coherency.

The generation of band-limited noises for given ¯γpq,jk(Δ*, f*) can be carried out by applying the spectral representation method (Shinozuka and Jan, 1972) and using Eigen decomposition (Shinozuka et al., 1990) or square root decomposition.

Using the sampled time series of the noises, the analyses for Steps (b) and (c) are straight forward. To incorporate the spatial

**TABLE 1** | Typical model parameters for spatial correlation and coherency models based on Hong and Liu (2014).

correlation structure of the FAS for two horizontal orthogonal directions, it is considered that the FAS of the *p*-th direction at the *j*-th site equals *r*Ap*,j × yj*(**M**, *Rj*, *f*), where *r*Ap*,j* (*p* = 1, 2) denotes the correlated random (scaling) disturbance of *yj*(**M**, *Rj*, *f*). Similar to the case of the spatial coherency, the (intraevent) correlation coefficient between ln(*r*Ap*,j*) and ln(*r*Aq*,k*) for *j* and *k* = 1,*. . .*, *n*R, *p* and *q* = 1 or 2, denoted as ρmn,jk(Δ) is defined by six elements. The results of statistical analysis (Liu and Hong, 2013; Hong and Liu, 2014) suggested that ln(*r*Ap*,j*) could be modeled as a normal variate with the SD equals 0.523 and ρmm,jk(Δ) can be modeled using,

$$\mathfrak{p}\_{\text{mm,jk}}(\Delta) = \exp\left(-a\_1 \Delta^{b\_1}\right),\tag{4}$$

for the record components along the same direction and

$$\mathfrak{p}\_{\text{mn,jk}}(\Delta) = r\_0 \exp\left(-a\_2 \Delta^{b\_2}\right),\tag{5}$$

for the record components along the orthogonal direction, where *a*1, *b*2, *r*0, *a*2, and *b*<sup>2</sup> are model parameters. Typical values of *r*0, *a*2, and *b*<sup>2</sup> are listed in **Table 1**. The suggested values shown in **Table 1** are developed based on records obtained from stations with separation greater than 100 m, which are considered to be adequate for the present study, although parameters for a closely separation (i.e., Δ *<* 100 m) can be found in Liu and Hong (2015b).

Samples of *r*Ap*,j* can be simulated based on the above specified probabilistic model and the values of *r*Ap*,j × yj*(**M**, *Rj*, *f*) (i.e., reference FAS) can be calculated. Using the obtained reference FAS for the *p*-th direction at the *j*-th site, scaling of the FAS is carried out in Step (d), and the application of inverse Fourier transformation in Step (e) results in a set of record components for a considered scenario event. Multiple simulation cycles for the same scenario events can be carried out by repeating Steps (a)–(e).

#### **Scenario Earthquake**

It must be noted that although the selection of a scenario event is not a trivial task, it can be carried out based on seismic hazard deaggregation for a specified probability of exceedance (Bazzurro and Cornell, 1999; Hong and Goda, 2006). It can also be assigned based on engineering judgment and emergency preparedness planning requirements. Alternatively, it can be identified based on geological and seismological investigation of seismic source zones if results of such investigation are available. In this study, a scenario earthquake event described in Atkinson and Macias (2009) was adopted. The scenario is an interface event with a moment magnitude **M**8.5 and a rupture plane of 380 km *×* 90 km, placed symmetrically about a perpendicular line from the Juan de Fuca trench to the city of Vancouver. The top corner of the fault plane is placed at [47.1°N, 124.5°W], and 10 km deep from the sea level. The strike and dip angle are equal to 310° and 10°, respectively. A map showing the surface projection of the rupture plane is produced in Figure 9 in Atkinson and Macias (2009). The parameters used in the finite-fault model for this event and the site amplification factors are shown in **Tables 1** and **2** in Liu and Hong (2015a). These model parameters differ from those used by Atkinson and Macias (2009) because a newer version of the program for the finite-fault model that included several changes (Atkinson et al., 2009; Boore, 2009) was employed.

The local site condition in downtown Vancouver is considered to be site class C according to NEHRP (National Earthquake Hazards Reduction Program) (Cassidy and Rogers, 2004), where *V*S30 (average shear wave velocity for the top 30 m soil) ranges between 360 and 760 m/s (NRCC, 2005). Therefore, it is assumed that *V*S30 = 414 m/s is adequate for sites located in downtown Vancouver, which is consistent with the amplification parameters considered in Liu and Hong (2015a).

**Figure 2** shows the comparison between response spectra based on 100 simulation cycles of the simulation in this study with that estimated by Atkinson and Macias (2009) for (49.25°N, 123.13°W). The median spectrum for the 100 simulation cycles and the spectra corresponding to 84th and 16th percentile are also included in **Figure 2** to illustrate the dispersion due to simulation. The comparison indicates an adequate match that justifies the simulation method at a single site.

#### **Non-Linear Inelastic 2DOF Systems As Proxy to Buildings and Damage Index**

If a building is approximated by a non-linear inelastic SDOF system with Bouc–Wen hysteretic model under unidirectional excitations as was done in Liu and Hong (2015a), the governing equation is expressed in the following using the normalized displacement:

$$\begin{aligned} \ddot{\mathfrak{u}}\_{\rm x} + 2\xi\_{\rm x}\mathfrak{o}\_{\rm nx}\dot{\mathfrak{u}}\_{\rm x} + \mathfrak{a}\mathfrak{o}\mathfrak{o}\_{\rm nx}^2\mathfrak{u}\_{\rm x} + (1-\mathfrak{a})\mathfrak{o}\_{\rm nx}^2\mathfrak{u}\_{\rm xx} &= -\ddot{\mathfrak{u}}\_{\rm x}(t)/\Delta\mathfrak{x}\_{\rm x} \\ \dot{\mathfrak{u}}\_{\rm xx} = \frac{1}{1+\mathfrak{G}\_{\rm \eta}\varepsilon\_{\rm nx}} \left[\dot{\mathfrak{u}}\_{\rm x} - (1+\mathfrak{G}\_{\rm \nu}\varepsilon\_{\rm nx})\,\mathfrak{u}\_{\rm xx}I\_{\rm x}\right] \\ \varepsilon\_{\rm nx} = (1-\mathfrak{a})\int\_{0}^{T} \dot{\mathfrak{u}}\_{\rm x}\mathfrak{u}\_{\rm xx}dt, \end{aligned} \tag{6}$$

where *I<sup>x</sup>* = μ˙ *x* μ*zx n−*1 (β + γ sgn(μ˙ *<sup>x</sup>*μ*zx*)), μ and μ<sup>z</sup> are the displacement and hysteretic displacement normalized by the yield displacement capacity of the inelastic SDOF system, Δ*<sup>Y</sup>* (i.e., μ = *u*/Δ*<sup>Y</sup>* and μ*<sup>z</sup>* = *z*/Δ*Y*, in which *u* and *z* are the displacement and hysteretic displacement of the SDOF system, respectively); ω<sup>n</sup> = (*k*/*m*) 0.5 is the natural vibration frequency, in which *k* and *m* are the stiffness and mass of the system; *ü*g(*t*) is the ground acceleration time history; *ε<sup>n</sup>* is the normalized dissipated energy through hysteresis; α, β, γ, and *n* are shape parameters in which β + γ = 1, α controls the post-yield stiffness, and *n* controls the smoothness of the transition from linear elastic to non-linear inelastic responses; δ<sup>η</sup> and δ*<sup>v</sup>* are stiffness and strength degradation parameters, respectively. The defined symbols with additional subscript *x* are used to denote that they represent the quantities associated with *X*-axis; and an overdot on a variable denotes its temporal derivative.

The yield displacement Δ<sup>Y</sup> of the non-linear inelastic model could be approximately related to the seismic design requirements (NRCC, 2005), where the minimum required design base shear force *V*<sup>d</sup> is given by *V*<sup>d</sup> = *C*s*W*, *W* is the total weight of the structure and *C*<sup>s</sup> is the design base shear coefficient given in **Table 2** for different building types. It can be shown that Δ<sup>Y</sup> is

$$
\Delta \mathbf{y} = \mathcal{R} \mathbf{N} \mathbf{C} \mathbf{S} \,\mathrm{W}/k \tag{7}
$$

where *R*<sup>N</sup> is the coefficient taking into account that the actual yield strength of a designed structure is greater than *V*d. μ<sup>R</sup> and *R*<sup>N</sup> are considered to be lognormally distributed with mean values shown in **Table 2** and coefficient of variation (cov) of 0.3 and 0.15 (Ellingwood et al., 1980; Ibarra, 2003), respectively. μ<sup>R</sup> and *R*<sup>N</sup> are assumed to be independent for each building.

As mentioned in the Section "Introduction," the use of the non-linear inelastic 2DOF system to represent a building is more realistic than the use of SDOF system, especially if bidirectional horizontal ground motions are considered. In such a case, the governing equation in terms of normalized displacements can be represented by Park et al. (1986), Yeh and Wen (1990), Lee and Hong (2010):

μ¨*<sup>x</sup>* + 2ξ*x*ω*nx*μ˙ *<sup>x</sup>* <sup>+</sup> αω<sup>2</sup> *nx*μ*<sup>x</sup>* + (1 *−* α)ω 2 *nx*μ*zx* = *−u*¨*gx/*Δ*Yx* μ¨*<sup>y</sup>* + 2ξ*y*ω*ny*μ˙ *<sup>y</sup>* <sup>+</sup> αω<sup>2</sup> *nx*μ*<sup>x</sup>* + (1 *−* α)ω 2 *ny*μ*zy* = *−u*¨*gy/*Δ*Yy* μ˙ *zx* = 1 1 + δη*ε<sup>n</sup>* [μ˙ *<sup>x</sup> −* (1 + δ*vεn*) μ*zxI*] μ˙ *zy* = 1 1 + δη*ε<sup>n</sup>* [μ˙ *<sup>y</sup> −* (1 + δ*vεn*) μ*zyI*] *ε<sup>n</sup>* = (1 *−* α) ∫*t* 0 ( μ*zx*μ˙ *<sup>x</sup>* + μ*zy*μ˙ *y* )( cos *n* θ  + sin*<sup>n</sup>* θ )<sup>2</sup>*/<sup>n</sup> dt* (8)

where *I* = μ˙ *x* μ*zx n−*1 [β + γ sgn(μ˙ *<sup>x</sup>*μ*zx*)] + μ˙ *y* μ*zy n−*1 [β+ γ sgn(μ˙ *<sup>y</sup>*μ*zy*)], θ = tan*−*<sup>1</sup> (μ*<sup>y</sup> /*μ*<sup>x</sup>* ), and the symbols defined previously but with an additional subscript *y* instead of *x* represents the quantities associated with the *Y*-axis. The solution of Eq. 8 can be used to evaluate the "normalized" displacement at time *t*,

$$\mu\_{\rm D}(t) = \left( \left| \mu\_{\rm x}(t) \right|^{n} + \left| \mu\_{\rm y}(t) \right|^{n} \right)^{1/n} \tag{9}$$

At the incipient yield, max(μD(*t*)) equals to 1.0; max(μ<sup>D</sup> (t)) represents the peak ductility demand if it is greater than 1.0. By considering that the ductility capacity equals μcap, collapse occurs

*a IBT is the building index. Building index is related to the structural and occupancy types defined in HAZUS-Earthquake [Federal Emergency Management Agency (FEMA) and the National Institute of Building Sciences (NIBS), 2003] (1* = *W1-RES1, 2* = *W1-RES1, 3* = *W2-RES3, 4* = *W2-COM1, 5* = *S4M-RES3, 6* = *S4M-COM4, 7* = *S4H-RES3, 8* = *S4H-COM4, 9* = *C2L-RES3, 10* = *C2L-COM1, 11* = *C2M-RES3, 12* = *C2M-COM4, 13* = *C2H-RES3, 14* = *C2H-COM4, 15* = *URMLR-RES3, 16* = *URMLR-COM1, 17* = *URMMR-RES3, 18* = *URMMR-COM2).*

*b #B* = *the number of buildings and #S* = *number of stories.*

*c The target C<sup>S</sup> is used to represent the seismic design level for existing buildings.*

if max(μ<sup>D</sup> (t)) is greater than μcap. Based on these consideration, for simplicity and being similar to the case of nonlinear inelastic SDOF system (Goda and Hong, 2008b), a damage factor is defined by using,

$$\mathfrak{S}\_{\text{DF}} = \max(\min(\mathfrak{S}\_{\text{Shift}}, 1), 0) \tag{10}$$

]/[μcap *<sup>−</sup>* <sup>1</sup>

] . If

where δShift = [ max(( μx *<sup>n</sup>* + μy *n* )<sup>1</sup>*/<sup>n</sup>* ) *−* 1

δDF equals 0, it implies the responses is within elastic range (or at most at incipient yield). Collapse is observed if δDF = 1.0, and partial damage occurs for δDF within (0, 1).

#### **Aggregate Seismic Loss for a Portfolio of Buildings**

One of the most difficult and important task in estimating seismic loss is to establish the damage cost function in terms of the damage level (e.g., in terms of damage factor defined in the previous section). This can be carried out if sufficient damage survey data from historical earthquakes are available. However, as the data are always scarce, the damage cost function is often established based on structural component testing results and expert opinion or judgment. A significant set of cost functions is available in HAZUS [Federal Emergency Management Agency (FEMA) and the National Institute of Building Sciences (NIBS), 2003]. These functions for several structural types and in terms of Canadian dollars are given in Goda and Hong (2008b). More specifically, it is considered that seismic losses associated with a building are categorized into three types: building-related loss *L*BL (δ), contents-related loss *L*CO(δ), and business-interruption related loss *L*BI(δ), where δ = δDF. These damage-loss functions can be expressed as,

$$L\_{\rm BL}(\mathfrak{G}) = \mathfrak{G}^{\mathfrak{\beta}\_{\rm BL}} L\_{\rm BL}(1), L\_{\rm CO}(\mathfrak{G}) = \mathfrak{G}^{\mathfrak{\beta}\_{\rm CO}} L\_{\rm CO}(1), \text{ and}$$

$$L\_{\rm BL}(\mathfrak{G}) = \mathfrak{G}^{\mathfrak{\beta}\_{\rm BL}} L\_{\rm BL}(1) \tag{11}$$

where the values of losses for the complete damage *L*BL(1), *L*CO(1), and *L*BI(1), as well as the model parameters β BL , β CO , and β BI are shown in **Table 2** for each building type. By using the damage-loss functions, the aggregate seismic loss *L* for *n*<sup>R</sup> buildings subjected to the scenario earthquake is calculated using:

$$L = \sum\_{j=1}^{n\_{\text{R}}} \left( L\_{\text{BL}}(\mathfrak{S}\_{j}) + L\_{\text{CO}}(\mathfrak{S}\_{j}) + L\_{\text{BI}}(\mathfrak{S}\_{j}) \right) \tag{12}$$

where δ*<sup>j</sup>* denotes the damage factor δDF for the *j*-th building. The maximum possible aggregate loss, *L*max, equals that calculated from Eq. 12 for δ*<sup>j</sup>* = 1.

#### **NUMERICAL EXAMPLE APPLICATION IN SEISMIC LOSS ESTIMATION**

In this section, an example application of the framework shown in the previous section is presented for a scenario seismic event and a portfolio of buildings. For the numerical analysis, the scenario seismic event of moment magnitude **M**8.5 elaborated previously is considered. The selection of hypothetical portfolio of 100 buildings as well as the analysis results are discussed in the following. Comparison of the aggregate seismic loss of the portfolio of buildings obtained under bidirectional excitations are compared with that obtained by considering unidirectional excitation.

### **Considered Portfolio of 100 Buildings**

A portfolio of 100 hypothetical buildings located in downtown Vancouver is considered for the numerical example. The sites of the buildings are randomly selected over a square area of 2.5 km by 2.5 km centered at (49.2°N, 123.2°W), which contains 4,000 property lots, each with an area of 25 m *×* 50 m. The selected locations are illustrated in **Figure 3**. Similar to Liu and Hong (2015a), the set of 100 buildings consists of 18 building types shown in **Table 2**. The buildings are of different structural types and occupancies (40 residential buildings and 60 commercial buildings). They are sampled based on the statistical information describing the existing building stocks in downtown Vancouver (Munich Reinsurance Company of Canada, 1992; Onur, 2001). As explained previously, the ductility capacity μcap is considered to be lognormally distributed with cov of 0.3. The mean of the ductility capacity for different building types are shown in **Table 2**. However, unlike the case in Liu and Hong (2015a), in this study, the capacities of a building in two orthogonal horizontal orientations (i.e., along *X*-axis and *Y*-axis) rather than in a single orientation are considered. The yield displacements ΔYx and ΔYx of each building are defined according to Eq. 7, which can be written as Δ<sup>Y</sup> = *R*N*C*S*gT*<sup>2</sup> <sup>n</sup>*/*(2π) <sup>2</sup> with *T*<sup>n</sup> denote the natural vibration period. For each structure, *T*<sup>n</sup> in two orthogonal directions are considered to be independent identically uniformly distributed with mean shown in **Table 2** and lower and upper bounds equal to minus and plus 10% of the mean value. This is to represent the fact that the *T*<sup>n</sup> along two horizontal orthogonal direction may differ. *R*<sup>N</sup> is considered to be independently lognormally distributed

with mean shown in **Table 2**, and cov equal to 0.15 as discussed previously. Only a single set of structural characteristics of the buildings are sampled and considered in the following numerical analysis.

## **Illustration of Simulated Ground Motion Components, Calculated Structural Responses, and Damage Cost**

For the 100 building sites marked on **Figure 3**, samples of simulated time histories for two building sites obtained by applying the procedure outlined in the previous section are presented in **Figure 4**. Note that since the adequacy of the adopted simulation procedure for ground motion records to match the target spatial correlation and coherency and FAS are already discussed extensively elsewhere (Hong and Liu, 2014; Liu and Hong, 2015a), they are not repeated in here. By applying these time histories to the buildings modeled as nonlinear inelastic 2DOF systems, the responses from the time history analysis is shown in **Figure 5** in terms of the normalized displacements μx, μy, and μD. The peak values for the normalized displacement and their corresponding time are also shown in the figure. The figure illustrates that the occurrence of peak demand along different component could

**FIGURE 5** | Nonlinear inelastic response of the buildings modeled as 2-degree-of-freedom systems, in terms of normalized displacement response μx, μy, and μD, under the bidirectional ground motions shown in **Figure 4**. The peak normalized displacement and its corresponding time is also shown.

differ. The total displacement demand μ<sup>D</sup> is greater than a single component.

Based on the non-linear inelastic response, the damage factor of the 2DOF system can be evaluated using Eq. 10. For the two examples given in **Figure 5**, Building 1 is a 2-story wood frame residential building (*T*nx = 0.366 s; *T*ny = 0.432 s); Building 2 is a 15-story reinforced concrete commercial building (*T*nx = 1.664 s; *T*ny = 1.552 s). The calculated damage factor, δDF = 0.0475 for Building 1 and δDF = 0.554 for Building 2.

# **Estimation of the Aggregate Loss**

The numerical calculation carried out in the previous section is repeated for all 100 considered building. Using the calculated damage cost for each of the 100 buildings, *L*BL(δ), *L*CO(δ), and *L*BI(δ), and summing them up according to Eq. 10, the value of the aggregate loss of the portfolio of buildings *L* is obtained. This obtained value represents a sample of the aggregate loss for the considered scenario seismic event because of the uncertainty in the ground motions even for the same scenario event. By repeating the above analysis 100 times for the same set of buildings, samples of *L* are obtained and shown in **Figure 6**.

The figure shows that the aggregated seismic loss for the 2DOF system is generally following a straight line on the Gumbel probability paper, with median equal to 0.26.

# **Effect of Approximating Building As SDOF Systems versus 2DOF Systems**

Now, reconsider the ground motions and the structures shown in **Figure 3** but only considering the excitations along the *X*-axis and the buildings modeled as SDOF with the structural properties along the *X*-axis as well. The calculated responses, the damage factors, and the damage cost for the two buildings are shown in **Figure 7**. Comparison of the results with those shown in **Figure 5** indicate that the profile of the response time history of a SDOF system is generally similar with that of a 2DOF system along the same direction. However, the maximum response could differ significantly because of the consideration of interactions between the two orthogonal directions.

By repeating the analysis for the building modeled as SDOF systems and subjected to the same set of excitations along the *X*axis used in 2DOF system case, samples of *L* are obtained and also presented in **Figure 6**. The results follow Gumbel distribution well with a median value of 0.18. Similar analysis is also carried out for *Y*-axis; the results are also plotted in **Figure 6**.

Comparison of the results for SDOF systems under unidirectional excitations versus 2DOF systems subjected to bidirectional excitations indicate that the simplification of using unidirectional excitations and SDOF models can underestimate the aggregated seismic loss in a scenario earthquake event. Such underestimation is somewhat more significant at upper tail where the probability of exceedance is small. This observation emphasizes that the importance of using properly simulated bidirectional ground

under the unidirectional ground motions shown in **Figure 4**. The peak normalized displacement and its corresponding time is also shown.

motion excitations and realistic 2DOF structural models in the seismic risk assessment of structures that are sensitive to the orientation of the excitations.

# **CONCLUSION**

In this study, we first simulate bidirectional spatially (partially) correlated and (partially) coherent ground motion records for a scenario Cascadia earthquake using a simulation procedure based on stochastic finite-fault model. The seismic loss of a portfolio of hypothetical buildings in downtown Vancouver under the simulated bidirectional horizontal ground motions is then estimated. Each building is modeled as a 2DOF system with different dynamic characteristics in two orthogonal horizontal directions. The hysteretic behaviors of the 2DOF systems are described by bidirectional Bouc–Wen model. The results indicate that if unidirectional ground motions and SDOF structural models are considered, the aggregated seismic loss could be

# **REFERENCES**

underestimated, emphasizing the importance of using realistic bidirectional ground motions that includes spatial coherency and correlation structure and modeling buildings as 2DOF systems with different characteristics in two horizontal directions. If the computation power is not limited, the framework presented in this study can be expanded to investigate the effect of tri-directional ground motions and more complicated structural models.

# **AUTHOR CONTRIBUTIONS**

TL carried out the analysis and wrote the manuscript. HH oversaw the project and revised the manuscript.

# **FUNDING**

The financial support received from the Natural Sciences and Engineering Research Council of Canada and the University of Western Ontario is gratefully acknowledged.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Liu and Hong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Tsunami Hazard Analysis of Future Megathrust Sumatra Earthquakes in Padang, Indonesia Using Stochastic Tsunami Simulation

#### *Ario Muhammad1,2\*, Katsuichiro Goda1 and Nicholas Alexander1*

*1Department of Civil Engineering, University of Bristol, Bristol, UK, 2Department of Civil Engineering, University of Narotama, Surabaya, Indonesia*

#### *Edited by:*

*Nikos D. Lagaros, National Technical University of Athens, Greece*

#### *Reviewed by:*

*Filippos Vallianatos, Technological Educational Institute of Crete, Greece David De Leon, Universidad Autónoma del Estado de México (UAEM), Mexico*

> *\*Correspondence: Ario Muhammad<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Earthquake Engineering, a section of the journal Frontiers in Built Environment*

*Received: 21 September 2016 Accepted: 07 December 2016 Published: 23 December 2016*

#### *Citation:*

*Muhammad A, Goda K and Alexander N (2016) Tsunami Hazard Analysis of Future Megathrust Sumatra Earthquakes in Padang, Indonesia Using Stochastic Tsunami Simulation. Front. Built Environ. 2:33. doi: 10.3389/fbuil.2016.00033*

This study assesses the tsunami hazard potential in Padang, Indonesia probabilistically using a novel stochastic tsunami simulation method. The stochastic tsunami simulation is conducted by generating multiple earthquake source models for a given earthquake scenario, which are used as input to run Monte Carlo tsunami simulation. Multiple earthquake source models for three magnitude scenarios, i.e., *M*w 8.5, *M*w 8.75, and *M*w 9.0, are generated using new scaling relationships of earthquake source parameters developed from an extensive set of 226 finite-fault models. In the stochastic tsunami simulation, the effect of incorporating and neglecting the prediction errors of earthquake source parameters is investigated. In total, 600 source models are generated to assess the uncertainty of tsunami wave characteristics and maximum tsunami wave height profiles along coastal line of Padang. The results highlight the influence of the uncertainty of the scaling relationships on tsunami simulation results and provide a greater range of tsunamigenic scenarios produced from the stochastic tsunami simulation. Additionally, the results show that for the future major earthquakes in the Sunda megathrust, the maximum tsunami wave height in Padang areas can reach 20 m and, therefore, significant damage and loss may be anticipated in this region.

Keywords: stochastic tsunami simulation, earthquake source modeling, uncertainty and sensitivity of tsunami hazard, Sunda megathrust, West Sumatra

#### INTRODUCTION

Located among three major plates, namely the Indian-Australian, the Pacific, and the Eurasian, Indonesia archipelago is one of the most seismically active regions in the world. In the last 20 years (1994–2014), 528 earthquakes occurred in Indonesia, i.e., about 26 earthquakes per year (United States Geological Survey (USGS), 2015). Sumatra Island is the most seismically active region since it is located at the interface between the Indian-Australian and Eurasian Plates. Two major seismic sources are the 1,900-km long Sumatran fault located along the center of Sumatra Island and the Sunda megathrust zone traversing more than 2,000 km along the western coast of Sumatra (Sieh and Natawidjaja, 2000). In the past two decades, several large earthquakes occurred along the Sunda megathrust including the Aceh-Andaman earthquake in December 2004 (*M*w 9.15), the Nias earthquake in March 2005 (*M*w 8.6), two earthquakes of Bengkulu in September 2007 (*M*w 8.4 and 7.9), the Mentawai tsunamigenic earthquake in October 2010 (*M*w 7.7), and the Indian Ocean earthquake in April 2012 (*M*w 8.6). Two of the most devastating earthquake events among those were the 2004 Aceh-Andaman earthquake triggering large tsunamis along the coastal line of Sumatra, Thailand, Sri Lanka, and India with the casualties of more than 250,000 people and the 2005 Nias earthquake which killed 2,000 people (Hsu et al., 2006; Banerjee et al., 2007).

The 2004 Aceh-Andaman earthquake was caused by the 1,600-km long rupture of the Sunda megathrust starting from North of Simeulue Island to North of Andaman Islands (Meltzner et al., 2006). This failure led to another 400 km rupture of the megathrust fault in the southern part of Simeulue Island (see **Figure 1**) causing the *M*w 8.6 2005 Nias earthquake (Briggs et al., 2006). A paleotsunami study based on a 1,000-year long record of tsunami deposits in North-West of Sumatra suggests that the occurrence interval of tsunamigenic earthquakes (*M*w 9.15) from the Sumatra-Andaman region is about 600 years (Monecke et al., 2008). Although the devastating tsunami event might not occur in the next few centuries in the Sumatra-Andaman segment, the ruptures of the megathrust fault have increased the failure probability of the Mentawai segment of the Sunda megathrust areas (see **Figure 1**) which is located in South of the fault rupture areas of the 2004 and 2005 events (Nalbant et al., 2005; Chlieh et al., 2008; Sieh et al., 2008; Collings et al., 2012, 2013). The past seismicity in the Mentawai segment indicates that there were two major tsunamigenic events occurred in 1797 and 1833 (*M*w ~8.8) that

affected the coastal areas of Padang and Bengkulu. Geodetic and paleogeodetic studies indicate that the slip deficit accumulated in the Mentawai segment has already exceeded the slip occurred during the 1797 and 1833 earthquakes (Collings et al., 2013). A large earthquake in 2007 (*M*w 8.4) that ruptured the Sunda megathrust near the 1833 rupture area was significantly smaller than the accumulated slip since the twin events of the 1797 and 1833 earthquakes (McCloskey et al., 2005; Nalbant et al., 2005). Hence, the possibility of earthquake and tsunami hazards in West of Sumatra from the Mentawai segment which has a recurrence interval of 200 years according to paleoseismological studies remains large (Sieh et al., 2008). In addition, the slip deficit is sufficient to generate a *M*w 8.8–9.0 earthquake (Zachariasen et al., 1999; Sieh et al., 2008).