Seismic damage accumulation in multiple mainshock–aftershock sequences

Earthquakes are generally clustered, both in time and space. Conventionally, each cluster is made of foreshocks, the mainshock, and aftershocks. Seismic damage can possibly accumulate because of the effects of multiple earthquakes in one cluster and/or because the structure is unrepaired between different clusters. Typically, the performance‐based earthquake engineering (PBEE) framework neglects seismic damage accumulation. This is because (i) probabilistic seismic hazard analysis (PSHA) only refers to mainshocks and (ii) classical fragility curves represent the failure probability in one event, of given intensity, only. However, for life cycle assessment, it can be necessary to account for the build‐up of seismic losses because of damage in multiple events. It has been already demonstrated that a Markovian model (i.e., a Markov chain), accounting for damage accumulation in multiple mainshocks, can be calibrated by maintaining PSHA from the classical PBEE framework and replacing structural fragility with a set of state‐dependent fragility curves. In fact, the Markov chain also works when damage accumulates in multiple aftershocks from a single mainshock of known magnitude and location, if aftershock PSHA replaces classical PSHA. Herein, this model is extended further, developing a Markovian model that accounts, at the same time, for damage accumulation: (i) within any mainshock–aftershock seismic sequence and (ii) among multiple sequences. The model is illustrated through applications to a series of six‐story reinforced concrete moment‐resisting frame buildings designed for three sites with different seismic hazard levels in Italy. The time‐variant reliability assessment results are compared with the classical PBEE approach and the accumulation model that only considers mainshocks, so as to address the relevance of aftershocks for life cycle assessment.


| INTRODUCTION
The performance-based earthquake engineering (PBEE) framework allows to assess the rate of earthquakes able to cause exceedance of a given seismic loss threshold to an engineering system of interest (e.g., a structure). 1 A more modest goal of PBEE is the failure rate, that is, the rate of earthquakes causing some undesired performance for the structure or λ f , which can be computed as follows: In the equation, P[ f |IM = im] is the probability that the structure fails the performance objective given a ground motion (GM) intensity measure (IM) value (im), to be obtained by structural response analysis, and dλ im,E is the differential of the hazard curve, that is, the function providing the rate of exceedance of im at the site of the construction, that is, λ im,E , from probabilistic seismic hazard analysis (PSHA). 2 The E subscript indicates that, in classical PSHA, the rate only accounts for the mainshock.
The rate from Equation 1 allows computing the probability that at least one earthquake, able to let the structure fail, occurs in (t,t+Δt): Equation 2 relies on two working assumptions that are also at the base of the modern seismic design approach: (i) structural failure is due to a single earthquake, that is, no progressive attainment of limit state due to damage accumulation is considered, and (ii) although earthquakes typically occur as clustered, only the mainshock, that is, the prominent event (to follow) of the cluster, is able to lead the structure to failure. The issue (i) is related to the fact that fragility curves provide the failure probability given the IM of a single earthquake, while the issue (ii) is related to the fact that PSHA, in computing λ im,E , only accounts for mainshocks. These hypotheses enable a mathematically convenient time-invariant approach to the problem according to which the occurrence of mainshocks follows a homogeneous Poisson process (HPP) characterized by a time-invariant rate. Fragilities are typically developed for the undamaged structure, based on the assumption that structural damage below the failure threshold is negligible or instantaneously repaired after the occurrence of a mainshock not causing failure. Nevertheless, earthquakes typically appear as clustered in time and space. Within each cluster, earthquakes are conventionally categorized in foreshocks, the mainshock, and aftershocks. * The mainshock is, classically, defined as the prominent magnitude event, whereas the foreshocks and the aftershocks are those preceding and following it, respectively. Therefore, in the case of long-lasting clusters, 4 the negligibility of damage accumulation within the cluster and the negligibility of the foreshocks/aftershocks' contribution to seismic hazard are questioned.
Mainly originated from the needs of life cycle assessment, damage accumulation has been the subject of a significant deal of research by those interested in stochastic modeling of degrading structures. In general, two types of deterioration mechanism are accounted for: (i) gradual deterioration of material characteristics (i.e., aging) and (ii) instantaneous damage due to sudden shocks. Some studies adopt analytically convenient reliability models, assuming independent damage increments (i.e., the damage increment does not depend on the damage history), to model one 5,6 or both issues. 7,8 However, state-of-the-art structural mechanics often requires to model structural degradation in a state-dependent fashion, that is, the future damage predictions depend on the state of the structure at the time of the evaluation. 9,10 In fact, if accumulating damage can be quantified by a number of discrete states, Markov chains can be an effective approach to solve the issue, as literature shows. 11 In particular, when earthquakes are of concern, Markovian modeling of damage accumulation replaces the failure rate from Equation 1 with a unit-time transition matrix, which collects the rates of seismic events causing eventual transitions from any damage state to any other. 12 To calibrate the stochastic model, the fragility curve must be replaced by a set of vulnerability models, that is, state-dependent fragilities, which provide the probability of exceeding a certain damage state, conditional on the GM intensity of the earthquake and the current damage state at the time of earthquake occurrence. If state-dependent fragility curves are combined with classical PSHA, the resulting transition matrix refers to damage accumulation in multiple mainshocks, as sketched in Figure 1a.
Regarding the fact that PSHA only considers mainshocks, it is discussed in a number of studies that the occurrence of seismic clusters can be treated in a similar fashion as the conventional PSHA, because the occurrence rate of the whole cluster is equal to that of the mainshocks. Then, the clusters follow the same Poisson occurrence process of the mainshocks. [13][14][15] It should be noted, finally, that if state-dependent fragility curves are combined with aftershock PSHA (APSHA) 16 instead, the Markov chain is able to model damage accumulation within a single aftershock sequence, following a mainshock of given magnitude and location, as shown already in literature. 12 Based on these grounds, the objective of the work presented hereafter is to formulate a Markovian model for seismic structural damage accumulation for both mainshocks and aftershocks. It is assumed that the structure, in a time interval of interest, can face multiple mainshock-aftershock sequences that are random in number and occur according to an HPP. It is also assumed that the path from as-new (AN) conditions to failure is represented by a series of damage states (DSs). In fact, the obtained Markov chain accounts, at the same time, for: (i) damage accumulation within any mainshock-aftershock seismic sequence, and (ii) damage accumulation in multiple sequences, as sketched in Figure 1b.
To illustrate the work, the remainder of the paper is structured in such a way that the Markov chain-based model for damage accumulation due to mainshocks only is given first, also recalling the basis of the classical PSHA. Subsequently, the damage accumulation within the aftershock sequence from a single mainshock, of given magnitude and location, is addressed. Then, the damage state transition matrix, for mainshock-aftershock sequences, is derived. The seismic reliability model is discussed through applications to a series of six-story reinforced concrete (RC) momentresisting frame code-conforming buildings, 17,18 to compare with the classical PBEE approach, and to the case when damage accumulation neglects aftershocks. Moreover, because the considered sites are exposed to different levels of seismic hazard, the applications also allow to assess the importance of the damage accumulation issue with respect to the seismicity of the construction location.

| DAMAGE ACCUMULATION ONLY IN MAINSHOCKS OR ONLY IN AFTERSHOCKS
This section recalls the Markovian modeling of damage accumulation as formulated by Iervolino et al. 12 In this reliability formulation, the hazard model refers to the classical PSHA (i.e., for mainshocks only) or to APSHA (i.e., for aftershocks only). According to the model, failure can be attained not only in a single earthquake of sufficient intensity, but also in multiple seismic shocks producing progressive damage (i.e., deterioration) on the structure. The hypotheses are as follows: (i) structural damage is represented by means of DSs, (ii) time is discretized in intervals of fixed width, and (iii) it is assumed that the probability of moving from a damage state to another, given the occurrence of an earthquake, only depends on the IM of the earthquake and on the damage state of the structure at the time of the seismic shock. The following first briefly recalls the basis of PSHA and then the reliability formulation for the Markov chain-based damage accumulation process. Finally, the application in conjunction with APSHA, that is, damage accumulation during the aftershock sequence of a specific mainshock, is recapped. F I G U R E 1 Damage accumulation for a structure subjected to (a) mainshocks only and (b) mainshock-aftershock sequences, following the concept of this study 2.1 | Shock occurrence according to PSHA According to PSHA, the occurrence of mainshocks at a site with magnitude larger than a minimum threshold (m E , min ) is described by an HPP, characterized by the v E rate. As already mentioned, in the earthquake engineering context, of particular interest is the exceedance of a certain IM E threshold. The occurrence rate of mainshocks producing such exceedance, λ im,E , can be computed according to Equation 3: In the equation, P[IM E > im j m E ,r E ] is the probability that the intensity threshold is exceeded given a mainshock of magnitude M E = m E , which is separated from the site by a distance R E = r E and can be obtained from a GM prediction equation (GMPE), which provides, in fact, the distribution of the intensity measure conditional to magnitude and distance. The term f M E ,R E is the joint probability density function (PDF) of M E and R E random variables (RVs) ranging in the {m E,min , m E,max } and {r E,min , r E,max } intervals, respectively. † It is useful to note that f M E ,R E is a joint PDF of magnitude and source-to-site distance RVs, conditional to the occurrence of a generic earthquake, that is, a mainshock of unspecified characteristics. Moreover, it is fairly easy to show that the absolute value of the derivative of the hazard curve divided by the v E rate represents the PDF of the IM E (which is obviously a RV), given the occurrence of a generic mainshock (i.e., of unknown magnitude and location). Such a PDF can be indicated as f IM E E j and will be useful in the following derivations of the Markovian stochastic model: In PSHA, Equation 3 is usually computed for a wide range of im values, so as to obtain the hazard curve for the site, which provides the exceedance rate as a function of the threshold. Finally, it is worth to recall that the rate λ im,E is also the rate of an HPP, that is, the process of occurrence of mainshocks causing exceedance of im at the site.

| Transition matrix in the case of mainshocks
The model developed in Iervolino et al. 12 for mainshocks refers to a discrete-time and discrete-state homogeneous Markov process. Given that time is discretized in intervals of fixed width, the time unit (e.g., 1 year), and that the domain of the considered structural damage measure is partitioned to have a finite number (n) of damage states, DS i , i = {1,2,…n}, from as-built (i.e., as-new) conditions to failure, the Markov chain is completely characterized by the [P E ] matrix; it has size n × n and contains the rate of earthquakes causing transition from any DS i ,8i = {1,2,…n}, to any DS j ,8j = {1,2,…n}. The rates in [P E ] approximate the unit time transition probabilities. Therefore, the probability that the structure in k time units is in any of the damage states can be computed as: In the equation, {P 1 (k) P 2 (k)Á Á ÁP i (k)Á Á ÁP n (k)} is the state vector after k time units, that is, P i (k) 8i = {1,2,…n} is the probability that the structure is in DS i at time k. Analogously, {P 1 (0) P 2 (0)Á Á ÁP i (0)Á Á ÁP n (0)} is the initial state vector, collecting the probabilities that the structure is in each of the damage states at the beginning of the interval (i.e., k = 0). Factually, Equation 5 replaces Equation 2, when damage accumulation is accounted for.
The [P E ] transition matrix can be computed as per Equation 6 where [P] is the transition matrix given the occurrence of a generic event, and [I] is the n × n identity matrix. Indeed, the ν E Á [P] term provides the unit time † Herein, the subscripts E and A are used referring to mainshocks and aftershocks, respectively. When no subscript appears for IM, it means that there is no need to make a distinction. probabilities (rates in fact) that the structure moves from any damage state to any damage state and a mainshock of unspecified magnitude and location occurs, whereas the (1 − ν E ) Á [I] term collects the probabilities that the structure remains in the same state and no earthquakes occur in the unit time interval.
To be completely specified, the Markov chain needs the ν E rate and the [P] matrix to be calibrated for the case under examination. The latter collects the P i,j probabilities that the structures is found in DS j ,8j = {1,2,…n}, after the earthquake, given that it was in DS i ,8i = {1,2,…n}, before the earthquake: : ð7Þ The diagonal elements in the [P] matrix are the probabilities that the earthquake does not cause a state change and can be computed as the complement to 1 of the transition probabilities on the same row, that is, The elements in the lower triangle of the matrix are set to 0 because of the irreversible nature of damage progression. Therefore, the nth state, representing failure, is an absorbing state, 19 that is, the structure cannot escape from it: P n,n = 1.
The computation of P i,j , i < j can be performed via Equation 8, which is an application of the total probability theorem. It requires f IM E E j (the PDF of IM E in a mainshock of unspecified magnitude and location) from Equation 4 and P [DS j |DS i , IM = im ], which is the probability that the structure in DS i moves to DS j , with j > i, given the im value of the mainshock intensity: The term P[DS j |DS i , IM = im ] can be derived via Equation 9, in which P[DS ≥ DS j |DS i , IM = im ] and P [DS ≥ DS j+1 |DS i , IM = im ] are the conditional probabilities that the structure in DS i reaches or exceeds DS j and DS j +1 (given earthquake intensity), respectively. In other words, P[DS j |DS i , IM = im ] can be computed from statedependent fragility functions: According to PSHA, the RVs representing GM intensities in different mainshocks are independently and identically distributed (i.i.d.), that is, their PDF is always f IM E E j . Thus, the P i,j probabilities do not change in time, leading to a homogeneous Markov chain, if aging is neglected (Equation 5). The procedure to get the state-dependent fragility functions, needed, in turn, to compute the P i,j probabilities using Equation 8, is outlined and discussed in Section 4.3.

| Transition matrix in the case of aftershocks
Yeo and Cornell 16 developed APSHA for short-term risk assessment, that is, to manage the engineering risk after a significant earthquake at the daily time scale. APSHA starts from the modified Omori law, 20 in which, given that a mainshock of magnitude M E = m E occurred at t = 0, the daily rate of aftershocks ν A M E j t ð Þ is given by Equation 10: In the equation, m A,min is the minimum considered aftershock magnitude, which cannot be larger than m E , and {a, b,c,p} are parameters. According to Yeo and Cornell, considering the interval (0, ΔT A ), where ΔT A is the (assumed) duration of the sequence, the RV counting the aftershocks, Þis characterized by the following probability mass function: From ν A M E j t ð Þ, it results that the daily rate of aftershocks causing exceedance of the im threshold at the site can be computed as: where is the magnitude-distance distribution of the generic aftershock. It is indicated that this distribution is conditional on the mainshock distance; in fact, it should be, more precisely, mainshock location. The distance limits {r A,min , r A,max } are the minimum and maximum source-to-site distance of aftershocks. The resulting rate in Equation 12 also characterizes a (nonhomogeneous) Poisson process.
It can be recognized that the IM in any generic aftershock to the same mainshock is an i.i.d. RV and its PDF, , can be obtained as: where f IM A jM A ,R A is the derivative of the IM distribution of aftershocks (from a GMPE). Then, the transition matrix given the occurrence of one generic aftershock from a sequence of a mainshock characterized by {M E = m E , R E = r E } is given by: : ð14Þ The terms of the P A M E ,R E j h i matrix are calculated as in Equation 15 using the same state-dependent fragility curves discussed in the previous section, yet integrated with the distribution of the aftershock intensity, f IM A M E ,R E j , rather than f IM E E j : ‡ ‡ For the sake of illustration, it is assumed herein that the fragilities are the same in the case of mainshocks and aftershocks occurrence. This is why there is no subscript to the intensity measure. However, some considerations about sufficiency 21 of the intensity measure should be verified before making such an assumption.
Clearly, the transition matrix P A M E ,R E j h i , because of APSHA, depends on the magnitude and location of the mainshock. The damage accumulation process within the aftershocks sequence is now a nonhomogeneous Markov chain characterized by a time-variant (i.e., daily) transition matrix; this is because, according to Equation 10, the aftershocks rate is a time-variant function. In particular, given a mainshock characterized by {m E , r E }, the unit time transition probabilities are given by § Consequently, the state probabilities after k time units are calculated as The initial state vector {P 1 (0) P 2 (0)Á Á ÁP i (0)Á Á ÁP n (0)} represents the probability that the mainshock triggering the aftershock sequence puts the structure in any damage state. In the context of this model, this vector is arbitrary (i.e., it must be assigned evaluating the structural conditions after the mainshock). However, the next section will show how to account for the uncertainty in the damage due to the mainshock combining the models of Sections 2.2 and 2.3.

| DAMAGE ACCUMULATION PROCESS IN MAINSHOCK-AFTERSHOCK SEQUENCES
Herein, the final Markov model is derived; it accounts for damage in a single mainshock-aftershock sequence and in multiple sequences. In the following, the damage in a single sequence is discussed first, then multiple sequences are addressed.

| Transition matrix in one sequence from a mainshock of given characteristics
It is worthwhile to first recall that during a sequence, there can be a random number of aftershocks, identified by the . Therefore, to account for the randomness in the number of aftershocks in each sequence, the following equation can be used to get the expected value, E[Á], of this matrix: ¶ which means that it is assumed that the sequence is made by a number of aftershocks equal to the mean provided by the modified Omori law for the mainshock characterized by {m E , r E }. However, the structure can make a transition also because of the mainshock and a transition matrix is also needed. It must account specifically for the fact that the earthquake has {M E =m E , R E =r E } characteristics. Therefore, such a matrix can be built as per Equations 7 and 8, but the state-dependent fragility function has to be integrated with the distribution of IM given {m E , r E }, which can be indicated as j ; that is: § The time scale, in the case of aftershocks, is typically the day rather than the year used for mainshocks. ¶ This is an approximation of E P A ME ,RE where it is obvious that f IM E M E ,R E j is available from a GMPE, which typically provides the mean and standard deviation of a Gaussian RV, representing the logarithm of IM conditional to earthquake magnitude, source-to-site distance, and other covariates, such as site's soil conditions. Collecting these probabilities in a matrix yields: 3.2 | Transition matrix for a generic mainshock-aftershock sequence and unit time transition matrix The total probability theorem allows to compute the transition probability given the occurrence of any sequences, [P S|E ], by integrating over the range of mainshock features: where the integral of the matrix is intended as the integral of each of its terms.
In the context of damage accumulation due to mainshock-aftershock sequences, one can apply Equation 21 to calibrate a homogeneous Markov chain similar to Equation 7. Because, as noted, the mainshock-aftershock sequences occur at the same rate of the mainshocks, unit time damage transition probability matrix for seismic sequences [P S ] is given, in analogy to Equation 6, by Equation 22:

| APPLICATION: DAMAGE ACCUMULATION FOR CODE-CONFORMING STRUCTURES
In this section, the Markovian model of degrading structures due to mainshock-aftershock sequences is illustrated through the application to a set of Italian code-conforming RC buildings. In particular, the structural design refers to three Italian sites with different seismic hazard according to PSHA. In order to investigate the effects of sequential shocks, the results of the long-term seismic risk assessment are compared with (i) the approach that neglects both the effects of aftershocks and of accumulating damage (Equation 2) and (ii) the case of accumulating damage only for mainshocks, as described in Section 2.2.

| Structural models and sites
The buildings under consideration were taken from those of the Rischio Implicito Norme Tecniche per le Costruzioni (RINTC) project, 22 which assessed the seismic risk of new buildings in Italy. Among all the buildings analyzed in the project, herein, the residential six-story RC infilled frame buildings, designed for three Italian sites (Milan, Naples, and L'Aquila), are considered. The sites are representative of low, intermediate, and high seismic hazards, respectively. Local soil condition C 23 is assumed for all the sites. Each building considered herein has a 5 × 3 bays moment-resisting frame characterized by regularity in plan and elevation (Figure 2a), a floor area of 21.4 × 11.7 m 2 , and story heights of 3.05 m (except the ground floor that has a 3.4 m height). It is worth noting that although the buildings share the same architectural characteristics, they have different structures. This is because they are designed for different sites. More specifically, according to the Italian code, the buildings were designed, with respect to the damage and the life safety (LS) limit states, based on the uniform hazard (design) spectra at the construction site, referring to 50-and 475-year exceedance return periods, respectively. 24 For the purposes of this study, nonlinear dynamic analyses of the structures were performed on the equivalent single-degree-of-freedom (ESDoF) systems based on the static pushover curve (SPO) of the original three-dimensional (3D) structural models, both constructed with OpenSees. 25 For each building, two ESDoF systems were calibrated according to Fajfar 26 in conjunction with the criteria in Baltzopoulos et al. 27 approximating the original frame model with a lumped mass multiple-degree-of-freedom system and assuming the first-mode load distribution for each horizontal direction (X, Y corresponding to Figure 2a). Figure 2b shows an example of the resulting backbone curve (L'Aquila, X direction) in terms of base shear (F) and displacement (δ) of the ESDoF.
Each backbone was rendered quadrilinear via a Monte Carlo-based optimization approach, 27 and the dynamic structural properties, as well as the backbone parameters, illustrated in Figure 2b, were derived: equivalent vibration period (T * ) and equivalent mass (m * ), modal participation factor (Γ), yield strength and displacement F * y , δ * y , postyielding hardening and softening ratio (a h , a c ), capping ductility μ c = δ * c =δ * y , residual strength (F p ), and global collapse ductility μ GC = δ * GC =δ * y . For the given SPO parameters, a moderately pinching, peak-oriented hysteretic behavior without any cyclic stiffness/strength deterioration 28 was applied; an example under cyclic loading is reported in Figure 2c. All parameters, including the yielding spectral acceleration at the equivalent vibration period Sa y T * À Á = F * y =m * , are summarized in Table 1. An equivalent viscous damping ratio (ξ * ) equal to 5% was applied to all models. 29

| Definitions of damage states
The damage index considered herein is a displacement-related engineering demand parameter. This illustrative application considered five DSs based on the backbone characteristics of the considered building: (i) AN, (ii) usabilitypreventing damage (UPD) as defined in the RINTC project; (iii) LS; (iv) near collapse (NC); and (v) global collapse (GC) as also defined in the RINTC project.
The threshold displacement identifying the UPD, indicated as δ 1 , was identified through a multicriteria approach on the original 3D structural models; it corresponds to the minimum displacement under multiple damage conditions that can impede the building occupancy after a seismic event. In particular, these conditions refer to the damage in main nonstructural elements, that is, masonry infills. 22 As it regards the displacement for GC, δ 4 , it corresponds to 50% of strength deterioration with respect to the maximum base shear on the SPO curve of the original structural model in each horizontal direction. The attainment of the LS performance level (δ 2 ) corresponds to the excursion to the residual plateau of the ESDoF system backbone. The threshold of NC performance level (δ 3 ) conventionally corresponds to two third of δ 4 . The AN state is consequently specified from no damage (zero displacement) up to the attainment of UPD.
In Figure 3, these displacement thresholds are illustrated for the three buildings of L'Aquila (AQ), Naples (NA), and Milan (MI) in both X and Y directions. The corresponding numerical values of displacement thresholds are reported in Table 2.

T A B L E 1
Dynamic and static pushover parameters of the equivalent single-degree-of-freedom systems

Dir.
T

| State-dependent fragility curves
To obtain the transition probabilities P[DS j |DS i , IM = im ] in Equations 8, 15, and 18, state-dependent fragility curves need to be computed. Using the numerical model of the structure and a suite of GM records, the classical fragilities, that is, transition probabilities of the AN structure (i = 1), were computed, via incremental dynamic analysis (IDA), from the response of the intact structure under a single shock. [30][31][32] Fragility curves for the structure in a damage state (i > 1) were computed via back-to-back IDA ** (B2BIDA). 33 This method involves two-step nonlinear dynamic analysis: the intact structure is first damaged to DS i by scaling each GM record, then DS j is attained further scaling each GM record applied to the structural model, which represents the structure in DS i .
All the analyses were performed using a set of 20 large-scale, short-distance, seismic events (moment magnitude within 6.5-6.9, closest distance to fault rupture within 15-33 km, recorded on firm soil) selected by Vamvatsikos and Cornell 28 and also used in Suzuki and Iervolino. 34 Each GM includes two horizontal components, namely, x and y. For each GM, the chosen IM is the maximum spectral acceleration (Sa) between the two horizontal components, at vibration period equal to 0.5 s: Sa(0.5s) = max(Sa x (0.5s), Sa y (0.5s)). IDAs and B2BIDAs were performed utilizing an opensource graphical user interface for dynamic analysis of SDoF systems in OpenSees, DYANAS. 35 The analyses have been performed with the aim of approximately accounting for the tridimensionality of the original structures, even if the used ESDoF systems are in fact representative of the unidirectional SPO curves. In the first step of the analyses, the structural responses in the two horizontal directions were examined independently by considering the two uncoupled ESDoF systems of each structure subjected to the two components of the selected GMs. For each GM, the two horizontal components in the as-recorded orientations were randomly associated to the X and Y direction of the structural model. 22,34 The considered GM was incrementally scaled, and the value of Sa(0.5s) causing the first exceedance of the considered DS threshold in either direction was determined (i.e., the lowest causing damage in any of the two directions). This allowed to derive a single fragility function for each damage state combining the results of the dynamic analyses in the two directions. As an example, one IDA curve for the building designed in L'Aquila is rendered in Figure 4a, showing the results from a two-component GM. The incrementally scaled intensity measure value is reported on the vertical axis of the plot, whereas the horizontal axis is the maximum transient displacement, δ. The threshold values identifying the limit states in both directions of the structure are also represented. For the considered record, the first exceedance of δ 1 occurs in Y direction and the corresponding value of Sa(0.5s) is represented by the star in the figure. On the other hand, the exceedance of δ 2 occurs in X direction and the marker identifying the corresponding Sa(0.5s) is the circle. The same identification is done for the other DSs.
To obtain the fragility curves of the structure in any DS, the same GM record set is always used. The structure put in one DS i by a record was subjected to each other record to get to the target DS j . Indeed, all possible combinations of the first and second shock were considered from the 20 GM records, thus leading to 400 (20 × 20) sequential excitations in total, for each pair of DSs.
For the sake of simplicity, it was assumed that when the structure reaches a given DS in one direction, the same DS is associated to the other direction. An example of B2BIDA is reported in Figure 4b, in which the same structure of L'Aquila, after being damaged to DS 2 (i.e., UPD), is subjected to the same GM considered in Figure 4a. As shown, the response in the two directions begins from residual displacements larger than zero. Indeed, in this particular case, DS 2  Figure 4c, the structure originally in DS 3 , subjected to the same GM incrementally scaled, reaches both DS 4 and DS 5 in Y direction. Denoting as Sa i,j (with i < j ≤ 5) the RV representing the intensity measure causing the structure transition from ith to jth DS or worse, state-dependent fragility functions were derived assuming the hypothesis of lognormal distribution of Sa i,j . The median and the logarithmic standard deviation of Sa i,j , denoted asSa i,j and σ lnSa i,j , respectively, were estimated via the IM-based approach, 27,36 and the state-dependent fragility curves were computed using Equation 23: Lognormal fragility parameters for all cases are shown in Table 3. When i = 1, the initial DS is AN, then it is the case of the classical fragility functions. As expected, the worse the initial state is, the smaller median acceleration level causing the transition to a worse DS can be observed. Comparing fragility parameters of different structures, it can be

State-dependent fragility parameters
Building site seen that, in most cases, the median increases with the seismic hazard of the site. This general trend is not respected in the case of the structure in L'Aquila in LS conditions for which medians for NC and GC are lower than the corresponding values computed for the structure in Naples. The reason can be found in Figure 3 that shows the limit state thresholds for the ESDoF systems: for both directions, the structure in Naples has a residual plateau longer than the structure in L'Aquila. From the computed state-dependent fragility curves, the transition probabilities between any two DSs can be retrieved via Equation 9.

| Seismic hazard assessment
As discussed in Sections 2.1 and 2.3, the Markovian modeling of degrading structures due to mainshock-aftershock sequences requires, from the hazard side, the following terms: (i) the occurrence rate of sequences (v E ), (ii) the conditional PDF of IM E , f IM E M E ,R E j , and (iii) the conditional PDF of IM A , f IM A M E ,R E j for each possible pair of mainshock features {M E , R E }. Whereas (i) and (ii) are involved in mainshock hazard, (iii) is needed for aftershocks. The required analyses were performed on the soil condition C for each of the considered three sites, L'Aquila, Naples, and Milan, via the REASESS software. 37 Details about the models adopted for PSHA and APSHA are described in the following subsections.

| Mainshocks
The seismic source model considered herein was the one by Meletti et al. 38 which features 36 seismic source zones (numbered from 901 to 936), as represented in Figure 5a. As discussed in Stucchi et al. 39, the model of Meletti et al. lies at the base of the official hazard assessment for Italy, which is, however, implemented via a 16-branch logic tree. The branch named 921 is the only one considered for this illustrative application. It features the GMPE of Ambraseys et al. 40 that is adopted also herein. According to branch 921, the magnitude-frequency distribution for each seismic zone was modeled based on the annual activity rates for discrete magnitude bins spanning from 4.15 to 7.45 (reported in Iervolino et al. 15 ).
For each site, only the seismogenic zones dominantly contributing the hazard were considered (see Figure 5a): in the case of L'Aquila, which is within the zone numbered as 923, this single zone was considered. For Naples, the considered zones were 924, 925, 927, and 928. Finally, for Milan, it was necessary to account for eight seismic zones-902, F I G U R E 5 (a) Seismic sources mostly contributing to the considered sites; (b) annual exceedance rates of Sa(0.5s) computed via probabilistic seismic hazard analysis for the three sites [Colour figure can be viewed at wileyonlinelibrary.com] 903, 906, 907, 909, 911, 913, and 915. This selection of zones for each site allows to reproduce the hazard of the comprehensive source model with negligible approximation (see also Iervolino et al.). 41,42 Figure 5b provides the hazard curves derived by PSHA, Equation 3, for the three considered sites: L'Aquila, Naples, and Milan. These will be used to compute failure probabilities according Equation 2. mainshock events (m A,min = m E,min ). The locations of aftershocks were assumed to be uniformly distributed within a circular area centered on the mainshock location, while the size was also determined by that of the mainshock. 44 The GMPE adopted for aftershocks is the same used for mainshocks, that is, Ambraseys et al. 40 Equation 20 involves

| Aftershocks
, which depends on the duration of the aftershock sequence ΔT A . The latter was defined arbitrarily equal to 90 days since the occurrence of the mainshock although, in principle, this duration could be mainshock-magnitude-dependent. However, in Iervolino et al. 15 it was verified that increasing the duration of the aftershock sequence up to, for example 1 year, yields minor differences in hazard results.

| Transition probability matrices and time-variant seismic risk
This section provides the resulting transition probability matrices and the time-variant seismic risk for the considered three sites. The results of the sequence-based seismic risk for the highest seismicity site, L'Aquila, are presented first, in comparison with those only accounting for mainshock events. Subsequently, the trends across the three sites with different seismic hazard levels are discussed.

| L'Aquila (high hazard)
The annual transition probability matrix for sequences, [P S ], is computed via Equation 22 and is shown in Table 4. For the sake of comparison, Table 4  It is also interesting to compare the transition rate from the AN to GC accounting for cumulative damage with the annual collapse rate of the AN building according to Equation 1. The latter is equal to 1.45E-04, that is, 0.7% exceedance probability in 50 years. The same value, 1.45E-04, is the annual failure rate computed accounting for damage accumulation and mainshock-only. This is expected because, in the hypotheses of this framework, damage can hardly accumulate in a single year (which is the assumed time unit), because there is a low occurrence probability of mainshocks. On the other hand, when the mainshock-aftershocks combined effect is considered, the unit time collapse rate of the AN structure increases up to 3.63E-04; this is due to the possible damage accumulation during a single seismic sequence, which is accounted for by Equation 20. All the curves monotonically increase up to 100 years. As shown in Figure 6a, in the mainshock-only case, the GC probability of the AN structure ranges between 1.45E-04 at k = 1 and 1.81E-02 at k = 100; the corresponding probability according to the sequence-based procedures ranges between 3.63E-04 and 3.90E-02. Failure probabilities derived by Equation 2 are equal to 1.45E-04 and 1.44E-02 after 1 and 100 years, respectively. The same figure shows that discrepancies between the mainshock-only and sequence-based cases are lower when damage to a milder damage state is considered. To give an example, if UPD is of interest, the probabilities after 100 years are equal to 1.49E-01 and 1.83E-01 in the case of mainshock-only and seismic sequence, respectively. As expected, increasing the initial DS, the GC probability increases. In the mainshock-only case, the failure probabilities in 100 years are equal to 1.86E-02, 1.14E-01, and 6.86E-01, for the initial damage state corresponding to UPD (Figure 6b), LS (Figure 6c), and NC (Figure 6d), respectively. The corresponding probabilities computed in 100 years under the sequence-based damaging process are equal to 4.11E-02, 2.45E-01, and 7.94E-01, respectively. Figure 7 shows the ratio between the probabilities of reaching a given damage state accounting for the seismic sequence and the corresponding mainshock-only probabilities, identified by S and E subscripts, respectively. First, it should be noted that in most of cases, the ratios do not significantly depend on time. Ratios of GC probabilities starting from AN or UPD state are comparable and range between 2.6 and 2.1. For the same initial DS, the ratios of probabilities of reaching and remaining in LS and NC are lower than one and almost constant with time, indicating that considering seismic sequence makes the attainment of such DSs less probable, because the structure tends to proceed to worse DSs. The ratio between the probabilities that the structure moves from LS to GC (Figure 7c) is the one in which the influence of the considered time window is the most evident. Such a ratio varies from 8.14 at k = 1 to 2.15 at k = 100; it suggests that the number of seismic events has a significant effect on the LS to GC transition probability (see residual strength plateau of Figure 3). Thus, for shorter time windows, aftershocks can significantly contribute to damage progress, leading to the high value of the ratio as shown in the figure. On the other hand, enlarging the time interval, the expected number of mainshocks increases along with the probability that mainshocks cause structural collapse. Thus, the ratio in the figure decreases. This reflects the combined effect of structural behavior, in terms of the DSs definitions and state-dependent fragility curves, with seismic hazard, as quantified in the annual transition matrices of Table 4. Indeed, the ratio at k = 1 is the one between the corresponding elements of the annual transition matrices computed for sequences and mainshock only, that is, 1.88E-03 and 2.31E-04, respectively. Finally, the ratio of GC probabilities, for a structure that is already in NC, is moderately larger than 1 (ranges from 1.4 to 1.2), indicating that the sequence effect may be minor, for an already heavily damaged structure. This section discusses the results for the buildings designed in the low-and intermediate-hazard sites, Milan and Naples, respectively. The unit time transition matrices, computed via the mainshock-only and the sequence-based approaches, are reported in Table 5 for both the sites. Consistent with the findings from the RINTC project, 22 although the three structures were all designed for GM with the same return periods at all sites, the transition matrices for the buildings in Milan and Naples show smaller transition probabilities than those for L'Aquila. Figures 8 and 9 show the time-variant probabilities of damage states for the Naples and Milan buildings, respectively, which were computed through Equation 5 in the mainshock-only and sequence-based cases. As indicated by the unit time transition matrices, it can be generally seen that structural damage is expected to proceed at a lower pace. Moreover, smaller discrepancies between mainshock-only and sequence-based reliability models were observed compared with the L'Aquila site. For Naples, the GC probabilities of the AN structure after 100 years computed according to Equation 2,Equation 5 in the mainshock-only case, and Equation 5 in the sequence-based case, are equal to 9.82E-04, 1.5E-03, and 2.3E-03, respectively. The corresponding probabilities computed for the Milan case are 4.91E-05, 5.12E-05, and 6.12E-05, respectively.
The ratios between the probabilities of reaching a given damage state accounting for the seismic sequence and the corresponding probabilities for mainshock-only are shown in Figure 10. The ratio of GC probabilities for a structure in AN or UPD in Naples varies between about 2 at k = 1 to 1.6 at k = 100. In the case of Milan, the corresponding ratio is about equal to 1.2. In the case of the structure assumed in LS state at k = 0 (Figure 10c), the probability ratios for Milan are larger than the corresponding ratios in Naples. The reason is in the definition of the DS threshold ( Figure 3) from which it is apparent that the distance between maximum displacement associated to LS and GC for Milan is significantly shorter than Naples. All the other transition probability ratios suggest that in the low-hazard site, the damage contribution of the aftershock sequences is minor, whereas in the case of intermediate seismicity, the transition probabilities accounting for the mainshock-aftershock sequences may be double with respect to the mainshock-only case. This is because, roughly speaking, the larger the classical (i.e., mainshock) hazard, the larger the number of aftershocks expected in the same time interval, then the larger the probability of accumulating damage.

| CONCLUSIONS
A Markov chain-based reliability model, for damage accumulation in structures subjected to mainshock-aftershock sequences, was presented. It is a discrete-time and discrete-state Markovian process, which accounts, at the same time, for both damage accumulation (i) within a single seismic sequence and (ii) among seismic sequences.
In the Markov chain, structural deterioration because of subsequent earthquakes is represented by a finite number of damage states from the AN state to collapse. B2BIDA is employed to compute the state-dependent fragility functions, which are needed to calibrate the model. Fragility functions are coupled with mainshock and aftershock hazard terms to formulate the required unit time transition matrix.
The model was illustrated through application studies using the ESDoFs of Italian code-conforming RC frame buildings located at three Italian sites representing different seismic hazard levels in Italy. The probabilistic seismic hazard assessments were carried out adopting the same models at the base of the current Italian seismic code. For each building, five damage states were defined. All fragility functions were assumed to be lognormal and were obtained in terms of spectral acceleration (maximum between the two horizontal components of a GM record) at a period close to the fundamental vibration periods of the structures. As expected, the computed state-dependent fragility curves under the same ground motion inputs for the three different sites generally show the median IM level causing the transition between any two DSs increasing with the seismic hazard at the site. This is in accordance with the structural design reflecting the local seismicity.
The sequence-based time-variant structural reliability was compared, for each site, with that obtained, (i) neglecting damage accumulation and (ii) considering accumulation due to mainshocks only. Notable remarks are summarized in the following.
• In the site of L'Aquila (high hazard), the probability that the structure in AN condition at time 0 collapses in 100 years considering damage accumulation only from mainshocks is 26% larger than the same probability computed F I G U R E 1 0 Ratios of probabilities of damage states for the buildings designed for Naples and Milan. GC, global collapse; LS, life safety; NC, near collapse; UPD, usability-preventing damage neglecting damage accumulation. Such an increment becomes 171% when mainshock-aftershock sequences are considered. • In Naples (intermediate hazard), the increment of the collapse probability in 100 years introducing damage accumulation is about 53% for mainshocks only and becomes 134% when mainshock-aftershock sequences are considered. • The structure designed in Milan (low hazard) shows the lowest differences when aftershocks' damaging contribution is considered. The GC probability in 100 years of the AN structure has an increment of approximately 4% when damage accumulation due to mainshocks is considered and an increment of approximately 25% when mainshockaftershock sequences are considered.
The model presented herein can be of help for life cycle assessment of structures in earthquake-prone areas, in a consistent manner with respect to the PBEE framework. It allows to account for damage accumulation considering mainshocks, aftershocks, and related sources of uncertainty.