Soil Dynamics and Earthquake Engineering

Conventional earthquake risk assessments use fragility and vulnerability models that are based on seismic demands from individual (mainshock) ground motions, and implicitly assume that a structure is intact before an earthquake hits. This study develops a suite of more realistic state-dependent seismic fragility and vulnerability models for a wide range of building taxonomies, leveraging state-of-the-art methods to account for dynamic damage accumulation in structures due to multiple earthquake events (i.e


Introduction
Recent worldwide earthquake events have highlighted the need to account for the effects of sequential earthquake-induced ground shaking in both seismic hazard analysis and risk calculations. For instance, studies suggest that the progressive damage of buildings during the entire 2010-2012 Canterbury New Zealand sequence (consisting of a 7.2 moment magnitude ( ) mainshock and three significant ≥ 5.9 aftershocks that occurred five, nine, and 15 months later) contributed more to the overall economic losses than those of the mainshock event alone [e.g., [1][2][3][4]. The 2016-2017 Central Italy sequence began with a 6.2 earthquake in August 2016, which was followed by two earthquakes ( 5.9 and 6.5) in October 2016 (considered triggered events; [5]). After the August 2016 event, 24% of buildings in Accumoli town had damages ranging from negligible to heavy, and 4% collapsed. Most partially damaged buildings in the town suffered additional damage in the two events of October 2016, increasing the total percentage of collapsed buildings to 65% [6].
The vast majority of available numerical fragility and vulnerability models are not suitable for capturing the damage and loss that can accumulate during sequential earthquake events, as they only consider * Corresponding author.
The importance of incorporating aftershock effects in damage and loss modeling has been recognized and increasingly acted upon in more recent years [e.g., [13][14][15]. Luco et al. [16] investigated the residual capacity of a mainshock-damaged three-story steel moment-resisting frame building adopting back-to-back incremental dynamic analysis (IDA; [17]), in which the same record is applied to represent both a mainshock and an aftershock ground motion. Specifically, a given record is first scaled to induce a target damage level (i.e., reach a prescribed level of maximum interstory drift ratio), then the record is applied again and scaled incrementally to calculate the post-mainshock residual capacity of the structure (i.e., the smallest spectral acceleration that would induce collapse in an aftershock). Using a similar method, Raghunandan et al. [18] computed aftershock fragility models https://doi.org/10.1016/j.soildyn.2023.107821 Received 6 October 2022; Received in revised form 17 January 2023; Accepted 2 February 2023 conditional on the mainshock damage. Non-linear dynamic analyses procedures (such as IDA or Multi-Stripe Analysis, MSA) involve significant computational effort [19], site-specific hazard-consistent record selection [e.g., 20,21], and excessive scaling of the seed records to various intensity levels (leading to ground-motion characteristics that do not represent the corresponding intensity level; e.g., [22]).
The cloud approach to seismic structural response quantification [23] is generally considered more appropriate for seismic risk assessments of building portfolios [e.g. , 12]. The original cloud analysis approach (i.e., considering only a single ground motion; [23]) uses univariate linear regression (in log scale) of an engineering demand parameter ( ) on a ground-motion intensity measure ( ). One advantage of the cloud approach is that it does not necessarily require a site-specific, hazard-consistent record selection [e.g., 24], and is based on no-to-moderate scaling of ground motions to describe a (possibly uniform) wide range of values [e.g., 25,26]. Some recent work has extended cloud analysis to the prediction of cumulative damage due to aftershocks, but only considering a limited number of building types [e.g., 12,27,28] rather than the wide range of building taxonomies that would be necessary for building portfolio analyses. These studies often only focus on the collapse limit state and generally do not derive vulnerability models that can be used in loss estimation exercises [e.g., 27,28]. In addition, they rely on groundmotion selection procedures based on random sampling approaches such as Latin Hypercube Sampling [e.g., 12,29], which do not necessarily cover a wide uniform range of s that may be important for fully capturing the cumulative damage effects of multiple earthquakes.
Furthermore, most studies on the residual capacity of mainshockdamaged buildings and related fragility model development have used peak structural response quantities like maximum displacement or interstory drift ratio as the [e.g., [30][31][32][33]. However, Gentile and Galasso [28] demonstrated that this approach may not appropriately capture damage accumulation. This is because peak response quantities may not monotonically increase with the duration and number of the applied ground motions, depending on the accuracy of non-linear hysteretic models (i.e., their ability to adequately capture cyclic and in-cycle stiffness and strength deterioration) and the intensity of the ground motions. Gentile and Galasso [28] introduced a new energybased framework for more accurately capturing accumulated damage, but this has only been demonstrated on a limited number of archetype reinforced concrete structures and deteriorating ordinary bridges [34], to date.
The aim of this study is to develop a large suite of state-dependent fragility and vulnerability models for a wide range of building taxonomies that might feature in a building portfolio, using a harmonized end-to-end methodology. The methodology leverages the newly developed energy-based framework for state-dependent fragility models by Gentile and Galasso [28], includes a collapse probability model, and incorporates a new standardized procedure to select ground motions with a wide-ranging uniform distribution of s. Using 1200 artificial ground-motion pairs, state-dependent fragility and vulnerability models are developed for 561 building classes (i.e., structural types) from the global database provided by the Global Earthquake Model (GEM; [10]) and are shared on a GitHub repository (see Data Availability). This study also discusses an approach for using the developed statedependent fragility and vulnerability models to quantify both damage and loss accumulation of building portfolios during ground-motion sequences. This approach is illustrated using ground motions from the 2010-2012 Canterbury sequence and a portfolio of building classes representative of the Christchurch central business district.
This study is organized as follows. Section 2 outlines the methodology for developing state-dependent fragility and vulnerability models. Section 3 applies this methodology to develop a global database of fragility and vulnerability models for 561 building classes. Finally, Section 4 introduces an approach to implement the developed statedependent fragility and vulnerability models in a portfolio loss assessment, before some concluding remarks (Section 5).  Fig. 1 shows a flowchart of the underlying end-to-end methodology for state-dependent fragility and vulnerability model development, given a set of building classes (and corresponding numerical models) that might be featured in a building portfolio of interest. The numerical models can be both single-(SDoF) or multi-degree-of-freedom (MDoF) models.

Methodology for state-dependent model development
A number of scaled ground-motion pairs (first ground motion, G1, and second ground motion, G2) are defined based on an input candidate ground-motion database and the characteristics of the building classes. This study considers G1 and G2 as two generic ground motions (i.e., they could represent foreshock-mainshock, mainshock-mainshock and mainshock-aftershock combinations), as opposed to most similar studies in the literature that only account for mainshock-aftershock sequences [e.g., 35]. Non-linear time-history analyses are then performed, subjecting the structural models of the building classes to each defined ground-motion pair. The probability of exceeding a set of input during G2 is estimated through cloud analysis [23] and the recently proposed hysteretic energy-based probability seismic demand model (PSDM; [28]).
are defined in terms of displacement-based thresholds, then converted to energy-based ones using the developed PSDM. The probability of collapse is also considered within the statedependent fragility model development. State-dependent vulnerability models are then developed using damage-to-loss models (DLMs) that associate each considered to a corresponding probability distribution of loss ratio [e.g., 36].
The use of ground-motion pairs to develop state-dependent fragility and vulnerability models is consistent with the Markovian assumption [11], which states that the additional damage experienced in a given earthquake only depends on the state of the structure right before the seismic event occurs, rather than the sequence of events that preceded it [37]. This assumption -which is not new in the context of modeling damage accumulation during earthquake sequences [e.g., 11,38] -greatly simplifies the development of both statedependent fragility and vulnerability models, as well as their implementation in portfolio risk assessments (Section 4).
The following sections present the details of each step of the proposed methodology. The related computer codes are developed in Python (see Data Availability).

Ground-motion selection
Pairs of unscaled records (G1 and G2 ) are selected from the input candidate ground-motion database and corresponding scaling factors ( 1 and 2 ) are defined. For each set of scaled G1 = 1 G1 and G2 = 2 G2 records, a discrete two-dimensional distribution with bins along the 1 and 2 axes (up to a maximum value, ) is computed (see Fig. 2). The simulated annealing method (i.e., dual annealing; [39,40]) is used to ensure these distributions uniformly cover the required range of 1 and 2 as best as possible. This method minimizes the following objective function: where ℎ is the number of ground motions falling in the th twodimensional bin. Eq. (1) represents the (log) mean squared error of ℎ with respect to its average. The simulated annealing method is run with the following constraints: (1) Values of 1 and 2 must strictly lie within [0.5,2.0], so that G1 and G2 are only moderately scaled, limiting potential bias in the fragility parameters [e.g., 41,42]; (2) Each unique G1 -G2 pair can only be selected once; (3) At least % of records selected for a given ground motion (G1 or G2 ) must be unique. This constraint helps towards achieving statistical independence of the data points used for fitting the PSDM. The higher the value of , the greater the computational time needed to minimize Eq. (1). Based on a number of trials with the candidate ground-motion database for state-dependent model development described in Section 3.1.2, = 70 was found to yield an adequate trade-off between excessive computational time and an adequate ground-motion selection for most explored cases, so this value is adopted herein; (4) Unless the considered structural models are three-dimensional, a G1 -G2 pair cannot include horizontal components from the same three-component ground-motion record [25]. For oneand two-dimensional structural models, this constraint also helps towards achieving statistical independence of the data points used for fitting the PSDM.
The proposed ground-motion selection procedure requires input values for and . = 20 is recommended for application of the methodology (and used in this work), to provide a reasonable discretization of levels.
should be selected such that each considered may be reached for G1 [e.g., 25], but its value is often further constrained by the limitations of the records in the input candidate ground-motion database (see Section 3.1.2). Input structural models are grouped together according to their fundamental periods and ductility (see 3.2.1), where each group is associated with one and a corresponding value. The ground-motion selection procedure is then carried out for each relevant { , } combination. Fig. 2 compares the proposed ground-motion selection procedure (lower panel) with the randomized approach based on Latin Hypercube Sampling used in Aljawhari et al. [12] (upper panel), for = 2g. The candidate ground-motion database for state-dependent model development described in Section 3.1.2 is used for both examples. A significant number of ground motions are concentrated at lower values in the upper plot. This is because the ground-motion database in Section 3.1.2 includes a large proportion of records with relatively low s; thus, a randomized G1-G2 pair selection is less likely to select large 1 − 2 pairs. The lower panel in Fig. 2 displays reasonably uniform distributions of 1 − 2 pairs between 0 and 2 g. The type used in this work is the average spectral acceleration at a range of periods of interest (AvgSA(0.2 -1.5 ), henceforth referred to as AvgSA( )), conventionally computed as the geometric mean of 10 discrete evenly-spaced spectral ordinates between 0.2 and 1.5 [e.g., 43], where is the fundamental period of the considered structure. AvgSA is relatively more sufficient [44] than other IMs, which means that site-specific, hazard-consistent ground-motion selection [e.g., 24] is relatively less crucial for cloud analysis [e.g., 45,46].
In summary, the proposed ground-motion selection procedure compiles artificial ground-motion pairs by moderately scaling groundmotion records in a uniform exploration of the 1 − 2 space. The resulting state-dependent fragility and vulnerability models can be used under multiple ground-motion sequences that could occur in the form of multiple mainshocks, mainshocks triggering additional earthquakes on nearby fault segments, mainshock-aftershock and aftershock-aftershock sequences. Real ground-motion sequences could be considered instead [e.g., 35,[47][48][49]. However, since recordings of sequences with high intensities are generally scarce, using real data would limit our ability to calibrate the models with observations of significant damage accumulation and jeopardize the statistical robustness of the models. Furthermore, the proposed procedure does not rely on site-specific hazard-consistent record selection methodologies [e.g., 21,35], making it relevant for large-scale (e.g., global) applications [10].

Probabilistic seismic demand model
Non-linear time-history analyses are performed for each input numerical model and each selected G1-G2 pair, incorporating a number of seconds of free vibration in between the ground motions. This freevibration time allows the structure to come to the rest position after G1, and it can be selected based on the fundamental period of the considered structure and engineering judgement. A free vibration time of 40 s is used in this work, which is reasonable for all building classes considered as part of the state-dependent model development process in Section 3 and widely accepted in the literature [e.g., 12,28].
The vector-valued PSDM developed by Gentile and Galasso [28] describes the structural demands obtained from the non-linear timehistory analyses using the cumulative hysteretic energy dissipated in the ground-motion sequence (i.e., G1 and G2) as the . This PSDM is the first physically-consistent model that accounts for the accumulation of damage, even though it is solely based on numerical analyses and requires further validation through experimental/field data. It is expressed as: in which is the cumulative hysteretic energy associated with G1 and G2, , is the hysteretic energy associated with ground motion GN (where is either 1 or 2), 1 is the maximum displacement (or drift) associated with G1, and 2 is as defined previously. All other variables ( , , 0 , , and ) are parameters estimated through regression. Eq. (2) with (a) 2 = 0 is the power-law relationship = , 1 = 1 , which can be calibrated with G1 data only; and (b) 1 = 0 is the power-law relationship = , 2 = 0 2 , which is commonly adopted for cloud-based single-ground-motion (e.g., mainshock) problems. The power-law relationship = , 2 = 0 2 can be calibrated using G1 data only because G1 can be interpreted as a G2 that follows a G1 producing zero drift [28]. The factor linearly reduces the intercept of the , 2 = 0 2 relationship. The steps used to fit Eq. (2) are similar to those described in Gentile and Galasso [28]: 1. Fit the relationship , 1 = 1 through linear least-squares regression in the log-log space, using G1 data only (i.e., 1 , , 1 ); 2. Fit the relationship , 2 = 0 2 through linear least-squares regression in the log-log space using G1 data only (i.e., 1 , , 1 ); 3. Determine the parameter from though non-linear least-squares regression in the log-log space, using 1 , 2 , and , 2 . Comparing two procedures for selecting ground-motion pairs (G1-G2) with values between 0 and 2 g, from the candidate ground-motion database described in Section 3.1.2. Upper panel: the histogram of G1-G2 obtained from the selection procedure specified in Aljawhari et al. [12]. Lower panel: distributions of G1-G2 obtained from the selection procedure proposed in this work. The average number of ground motions falling in the th 2D bin is equal to ∑ =1 ℎ ∕ , with bins along the 1 and 2 axes.
Step #3 represents a departure from the corresponding step in Gentile and Galasso [28], which uses non-linear least squares regression in linear space to calibrate . The modification to Step #3 is made because it is found that the log-space residuals of for the data of Section 3 are closer to a normal distribution than those in linear space, based on the Euclidean metric distance proposed by Chun et al. [50]. For practical implementation of the above steps, Eq. (2) is calibrated excluding data that result in (a) , 1 = 0 or , 2 = 0, as these represent cases where the structure does not incur any damage (i.e., linear elastic conditions); and (b) collapse, according to the details of Section 2.3.

Collapse probability model
In this study, collapse corresponds to a global dynamic instability of the numerical analysis (i.e., non-convergence) or the exceedance of a conventional maximum displacement [e.g., 10] during G1 or G2. The probability of collapse ( | 1 , 2 ) (where 1 corresponds to a damage state experienced during G1) is modeled with a formulation of the generalized logistic function [51], in which 1 and 2 are predictors (consistent with Eq. (2)): , , , and are model parameters, and * 1 is the lowest value of 1 that results in ( | 1 , 2 ) = 1 (i.e., corresponding to approaching 1 in Eq. (3)). Eq. (3) represents a univariate logistic regression model dependent on 1 in the 1 − ( | 1 , 0) plane, and a univariate logistic regression model dependent on 2 in the 2 − ( | 0, 2 ) plane (which is consistent with the logistic function commonly adopted for cloud-based single-ground-motion problems). Based on these observations and consistent with the PSDM proposed by Gentile and Galasso [28], the following steps are used to fit Eq. 2 ) )). The parameters and are estimated using G1 data only (i.e., 1 ), in line with the logic followed in Section 2.2.

State-dependent fragility models
State-dependent fragility models for different combinations of damage states due to G1 (i.e., 1 ) and G1+G2 (i.e., 1+ 2 ) are denoted as ( , representing the probability of exceeding the th damage state ( , with = 1 , … , and equal to the number of considered damage states) due to G1+G2, conditioned on 1 , 2 and no collapse (i.e., ). The fragility model can be expressed as: in which is the standard normal cumulative distribution function, is the mean and ln | 1 , 2 , is the dispersion of the fragility model. Note that 1+ 2 ≥ 1 (and the right hand-side of Eq. (4) reduces to 1 for = 1 ). If 1 = 0 (no damage occurs due to G1), Eq. (4) becomes equivalent to a single-ground-motion fragility model that only considers the possible occurrence of damage due to G2. ln | 1 , 2 , are computed with the equations provided by Gentile and Galasso [28]. Building-to-building variability ( − − ) is also included in the computation of ln | 1 , 2 , as suggested by Martins and Silva [10]. The value of is set equal to 0.3 in this work, per Martins and Silva [10].
The total probability of damage exceedance due to G1+G2 can be expressed as: in which ( | 1 , 2 ) is the probability of collapse calculated with Eq. (3).

State-dependent vulnerability models
This step of the methodology requires a DLM in the form of expected values, ( | ) and (optionally) coefficients of variation, ( | ), of the loss conditioned on each considered damage state. The DLM is usually region-and building-class-specific, and should be carefully selected based on the considered application [e.g., 53]. Given the time windows in which seismic sequences commonly take place (e.g., a few months to a few years; [54]), repair of damaged buildings between the occurrence of G1 and G2 is neglected in the development of the statedependent vulnerability models, per Papadopoulos and Bazzurro [15]. These vulnerability models can be written as [e.g., 55]: in which ( is the probability of occurrence of damage state ( = 1 , … , ) given 2 and the previous damage state 1 (derived from Eq. (4)), ( | 1 , 2 ) is the expected loss ratio given 1 and 2 , and is the probability density function of given . Eq. (6) could be extended to describe the variability around the expected loss ratio [see 55].

State-dependent model development
The methodology is now used to develop state-dependent fragility and vulnerability models for a large database of building classes. This database was developed by GEM based on a worldwide survey and comprises 561 building classes. Single-ground-motion fragility and vulnerability models for each building class have been developed by Martins and Silva [10], and then used for the assessment of economic losses in the global seismic risk model [56]. Each building class is identified following a simple taxonomy (see Table 1; [10]): (1) construction material; (2) lateral load resisting system; (3) ductility level; and (4) height (number of stories). For further details, the reader is referred to Martins and Silva [10].

Numerical models and damage state definition
Input numerical equivalent non-linear SDoF models [e.g., 57] are developed for each building class, using OpenSees ( [58]; see Data Availability). Each SDoF model is associated with a capacity curve provided by Martins and Silva [10] (see Fig. 3), expressed in terms of spectral acceleration (SA) versus spectral displacement (SD). The shape of the capacity curves depends on the expected reduction in the base shear capacity due to damage accumulation. For building classes where no significant base shear capacity reduction is expected (e.g., steel frames, timber and composite structures), a trilinear elastoplastic capacity curve is used. For other building classes, the capacity follows a quadrilinear curve (e.g., confined masonry or reinforced concrete with infills). The non-linear hysteretic behavior of the SDoF models is described with the Pinching4 material in OpenSees, using the cyclic degradation parameters from Martins and Silva [10]. For further details on the SDoF model development, the reader is referred to Martins and Silva [10].
Four building-class-specific are considered (from slight damage, 1 , to complete damage, 4 ). are defined based on the yielding ( ) and the ultimate ( ) displacements from the capacity curves (Table 2), which are derived from relevant literature [10]. Each damage-state displacement threshold is modeled with a lognormal distribution to account for aleatory uncertainty in the damage criteria. The median thresholds ( ) are provided in Table 2, and the coefficient of variation is set equal to 0.45 [10]. The conventional maximum displacement for the definition of collapse (see Section 2.3) is 1.5 times the displacement that results in 4 [10]. For further details on the definition of the , the reader is referred to Martins and Silva [10].
It is acknowledged that SDoF models might not accurately capture the response of individual buildings because they do not account for the effects of ground-motion directionality, aftershock polarity, structural Table 2 Median damage state thresholds (after Martins and Silva [10]) and description of the . and are the yielding and the ultimate displacements from the capacity curve, respectively. For further details on the definition of the damage states in terms of and , the reader is referred to Martins and Silva [10].

Damage state Description
Slight (  irregularities, higher mode behavior, localized failures, or the vertical component of earthquake-induced ground motions [e.g., 59,60]. However, SDoF models are generally deemed acceptable for describing the behavior of general building classes in the context of portfolio seismic risk assessments [e.g., 10,61,62], due to the typical lack of available detailed structure-specific data and limits on computational power [61]. Furthermore, using simpler structural modeling assumptions (i.e., SDoF instead of MDoF models) makes it easier to calibrate seismic risk models against real loss data [63].

Candidate ground-motion database
The input candidate ground-motion database consists of 7009 asrecorded three-component ground-motion records of the Pacific Earthquake Engineering Research (PEER) NGA-West2 database [64], which has been used to develop fragility models in several studies [e.g., 10,22]. The horizontal components of each three-component record are considered separately, such that the input database comprises 14018 single-component accelerograms. Fig. 4 provides a moment magnitude ( ) versus Joyner-Boore distance scatter plot corresponding to the candidate ground-motion database. ranges between 3.9 and 7.9, and the Joyner-Boore distance ranges between 0.1 and 600 km.

Damage-to-loss model
The DLM in Table 3 is used in this study for illustrative purposes only, in which the loss ratio corresponding to each follows a beta distribution. This DLM has been used by Martins and Silva [10] to develop single-ground-motion vulnerability models for all building classes considered in this study. The detailed definition of building class-specific DLMs is outside the scope of this study.

Ground-motion selection
AvgSA( ) are determined by clustering building classes according to ranges of fundamental period values (see Table 4), where each group is assigned an approximate value. A preliminary maximum AvgSA( ) value for each building class is first calculated as + 0.5, where is the yielding capacity of the considered building class, is calculated as: and is the ductility of the corresponding numerical model. , the transitional period, is taken to be 1s [65]. Eq. (7) is a conservative version of the strength reduction factor by Chopra [66], but only requires knowledge of . The preliminary values of maximum AvgSA( ) are then rounded up to the nearest 0.5 g unit, to finalize a limited set of maximum AvgSA( ) values across the considered building classes. Maximum AvgSA( ) values are also constrained between 1 g and 3.5 g for ≤ 0.8 s, and from 1 g to 3 g for > 0.8 s. Upper limit maximum values of 3.5 g ( ≤ 0.8 s) and 3 g ( > 0.8 s) are compatible with the maximum AvgSA( ) value of the input ground-motion database. 1 g is set as the lower limit of the maximum value, to avoid excessive extrapolation in the fitting of the PSDM (Section 2.2) and the collapse probability model (Section 2.3) [e.g., 25]. Fig. 5 provides the number of building classes in each final grouping, i.e., associated with each combination of assumed (Table 4) and maximum AvgSA( ).
The procedure detailed in Section 2.1 is run to select 1200 sets of G1 , G2 , 1 , and 2 for each building-class grouping. 1200 is deemed an appropriate number of ground-motion pairs to achieve statistically sufficient results in the structural response and is consistent with previous studies [e.g., 12,28]. Table 4 Period ranges used to associate building classes with an assumed fundamental period ( ) to compute AvgSA( ) (see Section 2.1).  (Table 4) and maximum AvgSA( ).

Additional steps: Example details
This section provides an example application of the rest of the model development methodology to building class CR-LFM-DUL-H2. This building class represents low-ductility buildings with a reinforced concrete frame and two stories (see Table 1) and features in several developed countries worldwide [67]. The capacity curve for CR-LFM-DUL-H2 is highlighted in Fig. 3. The numerical SDoF model for CR-LFM-DUL-H2 has a fundamental period equal to 0.56s. From Section 3.1.2, CR-LFM-DUL-H2 is assigned an assumed period of =0.6s and a maximum AvgSA( ) value equal to 1.5 g. The selected (moderately scaled) ground motions are shown in Fig. 6, where a wide range of AvgSA( ) values are observed for both G1 and G2. Fig. 7 provides the calibrated PSDM and collapse probability model for CR-LFM-DUL-H2. Fig. 7 shows that when the conventional maximum displacement of collapse for CR-LFM-DUL-H2 is exceeded (i.e., 1 ≳ 0.15), the collapse probability model is correctly equal to 1. As expected, the probability of collapse increases with increasing 1 (i.e., increasing the initial 1 ) and increasing 2 . Fig. 8 provides state-dependent fragility models developed for CR-LFM-DUL-H2 (the corresponding approximate mean, ln | 1 , 2 (0.6 s) , and dispersion, ln | 1 , 2 (0.6 s) , are presented in Table 5). As expected, ln | 1 , 2 (0.6 s) decreases with increasing 1 (e.g., the mean of the model for 4 | 3 is lower than that for 4 | 2 ). Fig. 9 provides the state-dependent vulnerability models for CR-LFM-DUL-H2, based on the fragility models shown in Fig. 8. As expected, ( | 1 , 2 (0.6 s)) increases for increasing 1 . For instance, ( | 1 , 2 (0.6 s) = 1 g)    is 9.5%, 52.5%, 106.1% greater than the single-ground-motion case (i.e., 1 = 0 ) when 1 is 1 , 2 , and 3 , respectively.

Proposed implementation procedure
This section details a procedure for using the developed statedependent fragility and vulnerability models to quantify both damage and loss accumulation during ground-motion sequences. The procedure is shown in Fig. 10 for one asset (with known location and building class); this process would be repeated for all assets of interest. It is assumed that the s of all events in the sequence are known at the asset's location.
values are generally available for past earthquakes at station sites [e.g., 68], otherwise they can be computed with simulation-based approaches [e.g., 69]. The procedure involves an iteration through realizations to sample different damage states of the asset after each ground motion. Each realization comprises of the following steps: • Set the current damage state and expected loss ratio to the respective values associated with the asset prior to the considered ground motion or sequence, i.e., = = 1 and ( ) = ( ) = ( | , 0). If the prior state of the asset is unknown, it is assumed that ( ) = 0 and = 0 ; • The total expected loss ratio increase for the th realization is computed as: The final expected loss ratio for the asset is computed as ( ) = ∑ =1 ( ) + ( ). Note that the event-specific expected increase in , ( ), can be computed as ( S. Iacoletti et al. Fig. 10. Proposed procedure for using the developed state-dependent fragility and vulnerability models to quantify both damage and loss accumulation during ground-motion sequences, for a single asset.

Example implementation
This section presents a simple case study -set in Christchurch, New Zealand at the time of the 2010-2012 Canterbury sequence -to illustrate the implementation procedure proposed in Section 4.1. The building classes used for the case-study portfolio are selected based on the findings of Uma et al. [70], which reports that the four predominant building typologies in the Christchurch central business district at the time of interest were: • Wood frames (W-LFM in Table 1 [70]. MUR-LWAL is always considered non-ductile (DNO in Table 1). For all building-typologies, building height is equally distributed between a low-rise story range (denoted as H2, 1-3 stories; [70]), and a mid-rise story range (denoted as H5, 4-7 stories; [70]). The building classes of the case-study portfolio are shown in Fig. 11; all correspond to building classes for which state-dependent fragility and vulnerability models were developed in Section 3.
The case-study sequence consists of ground motions from four earthquakes of the 2010-2012 Canterbury sequence. These four events are the mainshock ( 7.2 Darfield earthquake) and three aftershocks ( 6.2, 6.0, and 5.9, respectively). The response spectra are downloaded from the GeoNet website (see Data Availability) at the Christchurch Resthaven (REHS) station [68], which is located in the Christchurch central business district (longitude: 172.635 • ; latitude: −43.522 • ) and is characterized by site class D (deep or soft soil) as defined by the New Zealand Loadings Standard [71]. For convenience, it is assumed that all assets in the case-study portfolio are located at this site.
( ) and ( ) are calculated for each building class according to the approach outlined in Section 4.1 (henceforth called the proposed approach), and compared with values obtained using two alternative approaches: • The conventional approach used in seismic risk assessments (henceforth called alternative approach #1), which uses the single-ground-motion vulnerability model (i.e., ( ) = ( | 0, 2 )) to compute ( , ) for the mainshock event ground motion ( ) only, such that ( ) = 1 ∑ =1 ( , ); • A single-ground-motion approach that also considers aftershocks (henceforth called alternative approach #2). In this case, the underlying vulnerability calculations do not account for damage accumulation (i.e., the vulnerability module has no memory of the previous damage state). For this approach, ( ) is evaluated independently with the single-ground-motion vulnerability model (i.e., ( ) = ( | 0, 2 )) for each event and ( ) is computed accordingly. Table 6 reports ( ) and ( ) for the case-study portfolio building class CR-LFM-DUL-H2 highlighted in Section 3.2.2, along with 2 = AvgSA(0.6s) values corresponding to the four considered ground motions. Table 7 provides the ( ) values for all building classes of the case-study portfolio. From this table, the proposed approach produces ( ) values two to 20 times higher than alternative approach #1. This is at least partially explained by the fact that the first aftershock ( 6.2) produces larger values than the mainshock (see Table 6). The neglect of damage accumulation in alternative approach #2 can provide misleading ( ) values (including ( ) > 1), which may be either lower and higher than those of the proposed approach (see Tables 6 and 7).

Conclusions
This study develops state-dependent fragility and vulnerability models for a wide range of globally applicable building taxonomies, using Fig. 11. Building classes of the case-study portfolio (the taxonomy is described in Table 1). a harmonized end-to-end methodology. This methodology includes an innovative procedure (based on the simulated-annealing method) to select ground motions with an approximately uniform distribution of intensity measures that correspond to the characteristics of portfolio assets. The state-dependent fragility model development component leverages the energy-based probabilistic seismic demand model by Gentile and Galasso [28], but also incorporates a collapse probability model. The state-dependent vulnerability model development process relies on damage-to-loss models that are widely available in the literature. State-dependent fragility and vulnerability models are developed for 561 building classes of a global database provided by the Global Earthquake Model (GEM), using corresponding equivalent single-degree-of-freedom numerical models and a candidate groundmotion database comprising of around 7000 multi-component records from active shallow earthquakes. This study also presents a possible approach for using the developed state-dependent fragility and vulnerability models in portfolio loss assessments that incorporate seismic sequences. The proposed approach quantifies the time-dependent increase in the damage states of assets with the state-dependent fragility models, and estimates corresponding loss ratio increments with the state-dependent vulnerability models. The approach is illustrated with a portfolio of building classes representative of the old Christchurch central business district and ground motions of four earthquakes from the 2010-2012 Canterbury sequence. This approach is compared to two alternative approaches: (a) the conventional approach to seismic loss assessment that only considers the consequences of the mainshock ground motion; and (b) the approach of (a), but including a memoryless evaluation of damage from aftershocks. Approach (a) is the state-of-practice in portfolio seismic loss assessments. Approach (b) is simpler than the proposed approach, failing to consider typical real-life situations where assets remain unrepaired between events of a sequence. Results reveal that the proposed approach leads to final expected loss ratios two to 20 times higher than those of approach (a). The proposed approach also produces more reasonable results than those of approach (b), which can lead to expected losses that far exceed the replacement value of an asset. These findings underline the importance of a portfolio loss assessment appropriately capturing the time-dependent nature of the built environment's condition throughout a seismic sequence.
In summary, the model development methodology used in this study can support the development of extensive suites of state-dependent fragility and vulnerability models for quantifying the potential impact of damage accumulation in large-scale loss assessments of building portfolios. Although the approach proposed for implementing the developed state-dependent models in portfolio-level seismic risk assessments only captures seismic sequences in this study, it could be extended to account for the time-dependent condition of assets in more complete probabilistic seismic risk analyses. Region-specific models for building repair time could also be integrated into these analyses to realistically describe the cumulative damage of assets that remain unrepaired for years (as in the cases of the 2010-2012 Canterbury and the 2009 L'Aquila sequence; [72,73]).

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