Identifying Optimal Intensity Measures for Predicting Damage Potential of Mainshock–Aftershock Sequences

Large earthquakes are followed by a sequence of aftershocks. Therefore, a reasonable prediction of damage potential caused by mainshock (MS)–aftershock (AS) sequences is important in seismic risk assessment. This paper comprehensively examines the interdependence between earthquake intensity measures (IMs) and structural damage under MS–AS sequences to identify optimal IMs for predicting the MS–AS damage potential. To do this, four categories of IMs are considered to represent the characteristics of a specific MS–AS sequence, including mainshock IMs, aftershock IMs (i.e., IMMS and IMAS, respectively), and two newly proposed IMs through taking an entire MS–AS sequence as one nominal ground motion (i.e., IM1MS–AS), or determining the ratio of IMAS to IMMS (i.e., IM2MS–AS), respectively. The single-degree-of-freedom systems with varying hysteretic behaviors are subjected to 662 real MS–AS sequences to estimate structural damage in terms the Park–Ang damage index. The intensities in terms of IMMS, IMAS, and IM1MS–AS are correlated with the accumulative damage of structures (i.e., DI1MS–AS). Moreover, the ratio (i.e., DI2MS–AS) of the AS-induced damage increment to the MS-induced damage is related to IM2MS–AS. The results show that IM2MS–AS exhibits significantly better performance than IMMS, IMAS, and IM1MS–AS for predicting the MS–AS damage potential, due to its high interdependence with DI2MS–AS. Among the considered 22 classic IMs, Arias intensity, root-square velocity, and peak ground displacement are respectively the optimal acceleration-, velocity-, and displacement-related IMs to formulate IM2MS–AS. Finally, two empirical equations are proposed to predict the correlations between IM2MS–AS and DI2MS–AS in the entire structural period range.


Introduction
Past large earthquake events, such as the 1999 Chi-Chi earthquake [1], the 2008 Wenchuan earthquake [2], and the 2011 Tohoku earthquake [3] have shown that structures located in seismically active regions are commonly exposed to sequential-type excitations during a strong earthquake event. The cascading effects of mainshock (MS) and the following aftershocks (AS) could lead to a significant accumulation of structural damage, which tends to make a structure vulnerable to severe damage and even collapse. Nevertheless, current seismic design codes usually consider the design earthquake as a single event only, yet neglecting the detrimental effects caused by MS-AS sequences. Therefore, it is important to estimate the cumulative damage of structures under MS-AS sequences in an appropriate way, for accurately evaluating their dynamic performance subjected to a major earthquake event.
In recent years, numerous studies have been conducted to evaluate structural damage caused by MS-AS sequences using single-degree-of-freedom (SDOF) systems. These studies could be traced back to Mahin [4] who first examined the inelastic response of a SDOF system under a specific MS-AS sequence consisting of one MS and two ASs. It was reported that the ASs could lead to a significant increase in structural inelastic deformation and energy dissipation. After that, various structural performance indicators, e.g., peak displacement response [5][6][7], behavior factor [8,9], ductility demand [8,[10][11][12][13], input energy [14], damage index [15], and collapse capacity [16], have been adopted to examine the AS effect by subjecting SDOF systems to a suite of real [11,12] or artificial [8][9][10][11]15,16] MS-AS sequences. Recently, multiple degree of free dom systems (MDOFs) have been used to investigate the structural damage caused by MS-AS sequences; the studied structures include steel frames [17] bare and infilled reinforced concrete frames [18][19][20], wood structures [21], equipped structures [22], bridges [23] and other kind of structure [24]. The aforementioned studies compare the structural damage under MSs alone and MS-AS sequences, showing the significant effect of aftershocks on structural performance under earthquakes.
Accurate prediction of the damage potential of an earthquake has always been a well concerned issue in earthquake engineering [25]. To estimate the earthquake damage potential, it is necessary to introduce two intermediate variables, one representing the strength or severity of the earthquake and the other describing the actual damage or destructiveness due to the earthquake. Therefore, investigation of the interdependence between the earthquake intensity measures (IMs) and the structural damages plays a core role in the prediction of earthquake damage potential. An IM with a close dependence on structural damage can well reflect the damage potential of an earthquake. In recent years, different IMs have been linked with various structural damage indices to examine their correlation levels [26][27][28][29][30][31][32]. These studies have reported that the correlation between IMs and structural damage is affected by many factors, including the near-or far-fault ground motions [33,34], dimensionality of earthquake inputs [35,36], structural properties [37], and structural damage indicators adopted [38]. Due to the complex characteristics of ground motions, no IM alone can be regarded as the optimal one for assessing earthquake damage potential. Recently, several advanced IMs have been developed by combining classic IMs together to improve the capability of predicting damage potential [39][40][41].
In the studies stated above, however, the focus is mostly attached on the damage potential of single earthquake events without accounting for the effect of sequential earthquakes. In other words, the damage potential for the MS-AS sequences has not been paid adequate attention as that for the MSs alone. Recently, Elenas et al. [42] and Kavvadias et al. [43] examined the damage potential of MS-AS sequences by relating MS and AS intensities to the accumulative damage of reinforcement concrete (RC) frame structures under MS-AS sequences, respectively. Yet, only one RC frame structure was considered in both studies, leading to the limitation of the results to specific structural property. Moreover, only MS [42] or AS [43] IMs are considered to represent MS-AS sequences. This is unfavorable to reflect the real damage potential of a MS-AS sequence since the MS-AS-induced accumulative damage of structures is, to a large extent, dependent on the intensities of both MS and AS. Therefore, there is an urgent need to make a comprehensive investigation on the prediction of damage potential for MS-AS sequences.
The main purpose of this study is to identify a set of optimal IMs of MS-AS sequences which exhibit a high correlation to structural damages. To cover a wide range of general structural properties, several representative SDOF systems with varying hysteretic behaviors are taken as the study cases. The considered SDOF systems are then subjected to a suite of 662 sequential earthquake inputs selected from 13 earthquake events. The structural accumulative damage caused by the MS-AS sequences is obtained using the damage index proposed by Park et al. [44] and modified by Kunnath et al. [45]. Based on the obtained results, a comprehensive correlation analysis is conducted to search for the optimal IMs for predicting damage potential of MS-AS sequences.

Intensity Measures for MS-AS Sequences
In this study, four categories of IMs were used to represent MS-AS sequences based on 22 classic IMs. They are PGA [46], I A [47], P a [48], E a [49], a rms [50], a rs [51], I c [44], I a [52], PGV [46], P v [48], P D [53], E v [49], v rms [50], v rs [51], I v [52], I F [54], PGD [46], P d [48], E d [49], d rms [50], d rs [51], I d [52]. The 22 classic IMs have also been used in [38], which can be divided into three groups related to acceleration (see Figure 1), velocity (see Figure 2), and displacement (see Figure 3), respectively. It is noteworthy that the structure-specific IMs are not included in this study since they are commonly related to the structural fundamental period that would be elongated due to the MS-induced damage. In other words, two different structural periods are needed to define a structure-specific IM for the MS and the following AS. Thus, this study only adopts IMs associated with the basic properties of ground motions.
The considered four MS-AS IM categories include the mainshock IM, i.e., IM MS , and aftershock IM, i.e., IM AS , which are taken as the important intensity indicators for predicting the damage potential of MS-AS sequences according to Elenas et al. [42] and Kavvadias et al. [43], respectively. Besides IM MS and IM AS , two new types of IMs for MS-AS sequences were proposed herein, which are briefly represented by IM 1 MS-AS and IM 2 MS-AS . Using the considered 22 classic IMs, IM 1 MS-AS is defined by taking the entire MS-AS sequence as one nominal time history of ground motion acceleration.
In particular, the peak amplitude-related IMs, including PGA, PGV, and PGD, are used to determine IM 1 MS-AS by max (IM MS , IM AS ). As for the other IMs related to ground motion duration, they are used to define IM 1 MS-AS using the combined duration of the MS and the associated AS ground motions.
As for IM 2 MS-AS , it is defined by the ratio of IM AS to IM MS , i.e., IM 2 MS-AS = IM AS /IM MS . It is clear that IM 2 MS-AS is a dimensionless and compound parameter by combing IM AS and IM MS together. Tables 1-3 provide the considered IMs related to acceleration, velocity, and displacement, respectively, and the corresponding definitions for MS-AS ground motions. Table 1. Acceleration-related intensity measures (IMs) and the corresponding MS-AS IMs. are the duration for the MS, AS, and MS-AS sequence, respectively, between the corresponding times in terms of t 5 and t 95 corresponding to 5% and 95% of total I A , respectively. Table 2. Velocity-related IMs and the corresponding MS-AS IMs.     and IM AS by using the acceleration-, velocity-, and displacement-related IMs, respectively. It is found that most of the selected MS ground motions have larger intensities than the corresponding AS ground motions regardless of the IMs adopted, indicating significantly higher hazard caused by MSs compared to ASs. As observed, approximately over 90% (except only 89.32% of PGA MS values is greater than the corresponding PGA AS values, as shown in Figure 1a) of the considered MS ground motions show larger IMs than the corresponding AS ground motions. This result shows the larger damage potential of MSs than the corresponding ASs. Since the I A , E a and a rs have similar definition equations (see Table 1), their percentages corresponding to I MS A > I AS A (Figure 1b), E MS a > E AS a (Figure 1c), and a MS rs > a AS rs (Figure 1d), respectively, are almost the same. Similar results can be also observed between E v (Figure 2c) and v rs (Figure 2e) and that between E d (Figure 3c) and d rs (Figure 3e). Moreover, among the considered acceleration, velocity, and displacement IMs, the integral-formed IMs show slightly larger percentage with IM MS > IM AS than the other IMs.

Structural Damage under MS-AS Sequences
The damage index (DI) is used to simulate structural damage. In this study, four kinds of structural damage are considered, including the mainshock-induced damage for the initial system, i.e., DI MS , the AS-induced incremental damage for the mainshock-damaged system, i.e., DI AS|MS , the accumulative damage caused by MS-AS sequences, i.e., DI1 MS-AS , and the ratio of between DI AS|MS and DI MS , i.e., DI2 MS-AS = DI AS|MS /DI MS . In the following studies, DI1 MS-AS is related to IM MS , IM AS , and IM1 MS-AS , and DI2 MS-AS is related to IM2 MS-AS .
To determine DI1 MS-AS and DI2 MS-AS , the structural damage index proposed by Park and Ang [44] is adopted herein. The Park-Ang damage index is a linear combination of maximum deformation response and hysteretic energy dissipation. Therefore, it increases monotonically so that it is beneficial to quantify the accumulative damage caused by MS-AS sequences [15]. The Park-Ang damage index has been modified by Kunnath et al. [45], and the modified damage index, DI, is expressed by where uy and uu are the yielding and ultimate displacements of a specific SDOF system under monotonic loading, respectively; umax and EH are the maximum displacement and absorbed energy of the system under earthquake loading, respectively; and β is the energy consumption factor, which is taken as β = 0.15 herein following Zhai et al. [15].

Structural Damage under MS-AS Sequences
The damage index (DI) is used to simulate structural damage. In this study, four kinds of structural damage are considered, including the mainshock-induced damage for the initial system, i.e., DI MS , the AS-induced incremental damage for the mainshock-damaged system, i.e., DI AS|MS , the accumulative damage caused by MS-AS sequences, i.e., DI1 MS-AS , and the ratio of between DI AS|MS and DI MS , i.e., DI2 MS-AS = DI AS|MS /DI MS . In the following studies, DI1 MS-AS is related to IM MS , IM AS , and IM1 MS-AS , and DI2 MS-AS is related to IM2 MS-AS .
To determine DI1 MS-AS and DI2 MS-AS , the structural damage index proposed by Park and Ang [44] is adopted herein. The Park-Ang damage index is a linear combination of maximum deformation response and hysteretic energy dissipation. Therefore, it increases monotonically so that it is beneficial to quantify the accumulative damage caused by MS-AS sequences [15]. The Park-Ang damage index has been modified by Kunnath et al. [45], and the modified damage index, DI, is expressed by where uy and uu are the yielding and ultimate displacements of a specific SDOF system under monotonic loading, respectively; umax and EH are the maximum displacement and absorbed energy of the system under earthquake loading, respectively; and β is the energy consumption factor, which is taken as β = 0.15 herein following Zhai et al. [15].

Structural Damage under MS-AS Sequences
The damage index (DI) is used to simulate structural damage. In this study, four kinds of structural damage are considered, including the mainshock-induced damage for the initial system, i.e., DI MS , the AS-induced incremental damage for the mainshock-damaged system, i.e., DI AS|MS , the accumulative damage caused by MS-AS sequences, i.e., To determine DI 1 MS-AS and DI 2 MS-AS , the structural damage index proposed by Park and Ang [44] is adopted herein. The Park-Ang damage index is a linear combination of maximum deformation response and hysteretic energy dissipation. Therefore, it increases monotonically so that it is beneficial to quantify the accumulative damage caused by MS-AS sequences [15]. The Park-Ang damage index has been modified by Kunnath et al. [45], and the modified damage index, DI, is expressed by where u y and u u are the yielding and ultimate displacements of a specific SDOF system under monotonic loading, respectively; u max and E H are the maximum displacement and absorbed energy of the system under earthquake loading, respectively; and β is the energy consumption factor, which is taken as β = 0.15 herein following Zhai et al. [15].

SDOF Systems
A total of eight different hysteretic models are considered herein to simulate the nonlinear behavior of SDOF systems under MS-AS sequences. Similar SDOF model assumptions and considerations have been used in [41]. Figure 4 shows the hysteretic rules of the considered models, where the force and displacement are both normalized by the yield force and yield displacement, respectively. A polygonal type of hysteretic model is adopted in this study for illustration. Compared to the smooth type of hysteretic models (e.g., Bouc-Wen models [55][56][57][58], which have the advantage of using a set of differential or algebraic equations), the polygonal type of hysteretic model requires less parameters, which can also include both strength and stiffness deterioration. The SDOF systems partly or totally w/wo considering P-∆ effect, pinching, and damage-induced deterioration are used to account for the varying hysteretic behaviors of the SDOF system. Table 4 provides the considered eight SDOF models, where Delta is the parameter controlling the P-∆ effect; PinchX and PinchY are the pinching factors for strain (or deformation) and stress (or force) during reloading, respectively; and Damage1 and Damage2 denote the damage due to ductility and energy, respectively, which are used to account for the deterioration due to damage.
Among the considered models, the basic model, which is numbered as Model-1, is an elastic-plastic-with-hardening model (see Figure 4a). Based on the basic model, three important factors, including the P-∆ effect, pinching effect, and unloading degradation, are partly or totally incorporated in Models 2-8 to reflect more complex nonlinear behaviors of structures under seismic loading. In particular, Models 2-4 represent the bilinear hysteretic models by individually considering the P-∆ effect, pinching effect, and unloading degradation, respectively (see Figure 4b-d). For Models 5-7, two of the above-stated influence factors are accounted for in the hysteretic behavior, as shown in Figure 4e-g, respectively. Besides, Model-8 has the most complex hysteretic behavior since it collectively incorporates the P-∆ effect, pinching effect, and unloading degradation (see Figure 4h). It is noteworthy that the considered Models 2-8 are used for parametric study to examine the effect of hysteretic models of SDOF systems on the obtained results.

SDOF Systems
A total of eight different hysteretic models are considered herein to simulate the nonlinear behavior of SDOF systems under MS-AS sequences. Similar SDOF model assumptions and considerations have been used in [41]. Figure 4 shows the hysteretic rules of the considered models, where the force and displacement are both normalized by the yield force and yield displacement, respectively. A polygonal type of hysteretic model is adopted in this study for illustration. Compared to the smooth type of hysteretic models (e.g., Bouc-Wen models [55][56][57][58], which have the advantage of using a set of differential or algebraic equations), the polygonal type of hysteretic model requires less parameters, which can also include both strength and stiffness deterioration. The SDOF systems partly or totally w/wo considering P-Δ effect, pinching, and damage-induced deterioration are used to account for the varying hysteretic behaviors of the SDOF system. Table 4 provides the considered eight SDOF models, where Delta is the parameter controlling the P-Δ effect; PinchX and PinchY are the pinching factors for strain (or deformation) and stress (or force) during reloading, respectively; and Damage1 and Damage2 denote the damage due to ductility and energy, respectively, which are used to account for the deterioration due to damage.
Among the considered models, the basic model, which is numbered as Model-1, is an elasticplastic-with-hardening model (see Figure 4a). Based on the basic model, three important factors, including the P-Δ effect, pinching effect, and unloading degradation, are partly or totally incorporated in Models 2-8 to reflect more complex nonlinear behaviors of structures under seismic loading. In particular, Models 2-4 represent the bilinear hysteretic models by individually considering the P-Δ effect, pinching effect, and unloading degradation, respectively (see Figure 4bd). For Models 5-7, two of the above-stated influence factors are accounted for in the hysteretic behavior, as shown in Figure 4e-g, respectively. Besides, Model-8 has the most complex hysteretic behavior since it collectively incorporates the P-Δ effect, pinching effect, and unloading degradation (see Figure 4h). It is noteworthy that the considered Models 2-8 are used for parametric study to examine the effect of hysteretic models of SDOF systems on the obtained results.

Delta PinchX PinchY Damage1 Damage2
A total of 44 fundamental vibration periods ranging from 0.2 s to 6.0 s are considered for each SDOF system, consisting of 29 periods between 0.2 s and 3.0 s (with an increment of 0.1 s), and 15 periods between 3.2 s and 6.0 s (with an interval of 0.2 s). Four yield strength reduction factors, i.e., R = 2, 3, 4, and 5, are considered to account for the varying ductility levels of structures. Moreover, a constant damping coefficient ξ = 0.05 is utilized for all the considered SDOF systems according the past available studies, such as [15,21]. It is noteworthy that the analysis process of this paper is also suitable for the cases with other ξ values. However, the discussion on the effect of damping ratios is beyond this study and will be discussed in the future work. Consequently, a total of 1408 (8 hysteretic models × 44 vibration periods × 4 ductility levels) SDOF systems are used in subsequent calculations.

Selection of MS-AS Sequences
A set of 662 real MS-AS earthquake sequences were selected from 13 earthquake events based on the NGA-West2 database (https://peer.berkeley.edu/research/nga-west-2) (Ancheta et al [58]). The NGA-West2 database has been used to develop the latest generation of ground motion prediction equations for shallow crustal earthquakes (e.g., Campbell and Bozorgni [59], Du et al. [32] The detailed information of the selected ground motion records and the related earthquake events can be found in Zhu et al. [60]. It is noteworthy that, among the aftershock sequences triggered by a specific mainshock, only aftershock record with the largest magnitude was adopted. Moreover, the selected mainshock and associated aftershock ground motions should be recorded at the same station. Figure 5 shows the distribution of moment magnitudes Mw and rupture distances Rrup for the selected ground motions of MSs and ASs. The magnitude of MSs ranges from Mw = 5.8 to Mw = 7.62, while the magnitude of the corresponding ASs ranges from Mw = 5.01 to Mw = 7.14. As for rupture distances, the range of Rrup for MSs is between 0.32 km and 208.47 km and that for ASs is between 1.98 km and 218.16 km, respectively. As seen from Figure 5b, the Mw values of MSs are always larger than that of AS. This is compatible with the Bath's rule [61]. Figure 6 compares the elastic response spectra for the selected ground motions of MSs and ASs. It is found that the median response

Delta PinchX PinchY Damage1 Damage2
Model-1 0 A total of 44 fundamental vibration periods ranging from 0.2 s to 6.0 s are considered for each SDOF system, consisting of 29 periods between 0.2 s and 3.0 s (with an increment of 0.1 s), and 15 periods between 3.2 s and 6.0 s (with an interval of 0.2 s). Four yield strength reduction factors, i.e., R = 2, 3, 4, and 5, are considered to account for the varying ductility levels of structures. Moreover, a constant damping coefficient ξ = 0.05 is utilized for all the considered SDOF systems according the past available studies, such as [15,21]. It is noteworthy that the analysis process of this paper is also suitable for the cases with other ξ values. However, the discussion on the effect of damping ratios is beyond this study and will be discussed in the future work. Consequently, a total of 1408 (8 hysteretic models × 44 vibration periods × 4 ductility levels) SDOF systems are used in subsequent calculations.

Selection of MS-AS Sequences
A set of 662 real MS-AS earthquake sequences were selected from 13 earthquake events based on the NGA-West2 database (https://peer.berkeley.edu/research/nga-west-2) (Ancheta et al [58]). The NGA-West2 database has been used to develop the latest generation of ground motion prediction equations for shallow crustal earthquakes (e.g., Campbell and Bozorgni [59], Du et al. [32] The detailed information of the selected ground motion records and the related earthquake events can be found in Zhu et al. [60]. It is noteworthy that, among the aftershock sequences triggered by a specific mainshock, only aftershock record with the largest magnitude was adopted. Moreover, the selected mainshock and associated aftershock ground motions should be recorded at the same station. Figure 5 shows the distribution of moment magnitudes M w and rupture distances R rup for the selected ground motions of MSs and ASs. The magnitude of MSs ranges from M w = 5.8 to M w = 7.62, while the magnitude of the corresponding ASs ranges from M w = 5.01 to M w = 7.14. As for rupture distances, the range of R rup for MSs is between 0.32 km and 208.47 km and that for ASs is between 1.98 km and 218.16 km, respectively. As seen from Figure 5b, the M w values of MSs are always larger than that of AS. This is compatible with the Bath's rule [61]. Figure 6 compares the elastic response spectra for the selected ground motions of MSs and ASs. It is found that the median response spectrum of MS ground motions is significantly larger than that of AS ground motions, which is expected due to the greater shaking intensities of MSs than ASs. spectrum of MS ground motions is significantly larger than that of AS ground motions, which is expected due to the greater shaking intensities of MSs than ASs. For a given ground motion recorder, it could only detect seismic signals in a certain range of frequency due to its settings. This makes the recorded ground motions have a specific lowest usage frequency, i.e., flow. Therefore, for dynamic analysis, a ground motion can only be used for nonlinear SDOF systems with their vibration periods no larger than T = 1/flow. As for a specific MS-AS sequence, the flow values for the MS and AS ground motions are probably different. For simplicity, the representative flow value for a MS-AS sequence is taken as max (flow MS , flow AS ), where flow MS and flow AS are the flow values for the MS and AS ground motions, respectively. In this study, the highest structural period considered for the nonlinear SDOF systems is 6.0 s, and thus, the usable flow values of the selected MS-AS sequences should be no less than 0.166 Hz. Figure 7 shows the number of usable earthquake sequences at different periods with the consideration of flow requirement, where the number of usable sequences decreases with increasing periods of SDOF systems.  For a given ground motion recorder, it could only detect seismic signals in a certain range of frequency due to its settings. This makes the recorded ground motions have a specific lowest usage frequency, i.e., flow. Therefore, for dynamic analysis, a ground motion can only be used for nonlinear SDOF systems with their vibration periods no larger than T = 1/flow. As for a specific MS-AS sequence, the flow values for the MS and AS ground motions are probably different. For simplicity, the representative flow value for a MS-AS sequence is taken as max (flow MS , flow AS ), where flow MS and flow AS are the flow values for the MS and AS ground motions, respectively. In this study, the highest structural period considered for the nonlinear SDOF systems is 6.0 s, and thus, the usable flow values of the selected MS-AS sequences should be no less than 0.166 Hz. Figure 7 shows the number of usable earthquake sequences at different periods with the consideration of flow requirement, where the number of usable sequences decreases with increasing periods of SDOF systems.  For a given ground motion recorder, it could only detect seismic signals in a certain range of frequency due to its settings. This makes the recorded ground motions have a specific lowest usage frequency, i.e., f low . Therefore, for dynamic analysis, a ground motion can only be used for nonlinear SDOF systems with their vibration periods no larger than T = 1/f low . As for a specific MS-AS sequence, the f low values for the MS and AS ground motions are probably different. For simplicity, the representative f low value for a MS-AS sequence is taken as max (f low MS , f low AS ), where f low MS and f low AS are the f low values for the MS and AS ground motions, respectively. In this study, the highest structural period considered for the nonlinear SDOF systems is 6.0 s, and thus, the usable f low values of the selected MS-AS sequences should be no less than 0.166 Hz. Figure 7 shows the number of usable earthquake sequences at different periods with the consideration of f low requirement, where the number of usable sequences decreases with increasing periods of SDOF systems.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 21 spectrum of MS ground motions is significantly larger than that of AS ground motions, which is expected due to the greater shaking intensities of MSs than ASs. For a given ground motion recorder, it could only detect seismic signals in a certain range of frequency due to its settings. This makes the recorded ground motions have a specific lowest usage frequency, i.e., flow. Therefore, for dynamic analysis, a ground motion can only be used for nonlinear SDOF systems with their vibration periods no larger than T = 1/flow. As for a specific MS-AS sequence, the flow values for the MS and AS ground motions are probably different. For simplicity, the representative flow value for a MS-AS sequence is taken as max (flow MS , flow AS ), where flow MS and flow AS are the flow values for the MS and AS ground motions, respectively. In this study, the highest structural period considered for the nonlinear SDOF systems is 6.0 s, and thus, the usable flow values of the selected MS-AS sequences should be no less than 0.166 Hz. Figure 7 shows the number of usable earthquake sequences at different periods with the consideration of flow requirement, where the number of usable sequences decreases with increasing periods of SDOF systems.

Correlation Coefficient
The Pearson correlation coefficient ρ is calculated for four IM-DI pairs in the logarithm space, including IM MS vs. respectively. The basic system with R = 2 at T = 0.2 s, 1.0 s, and 5.0 s are taken herein for illustration to represent the structure with short (acceleration-sensitive), medium (velocity-sensitive), and long (displacement-sensitive) vibration periods, respectively (Riddell [38]). Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 21

Correlation Coefficient
The Pearson correlation coefficient ρ is calculated for four IM-DI pairs in the logarithm space, including IM MS vs. DI1 MS-AS , IM AS vs. DI1 MS-AS , IM1 MS-AS vs. DI1 MS-AS , and IM2 MS-AS vs. DI2 MS-AS , respectively. The basic system with R = 2 at T = 0.2 s, 1.0 s, and 5.0 s are taken herein for illustration to represent the structure with short (acceleration-sensitive), medium (velocity-sensitive), and long (displacement-sensitive) vibration periods, respectively (Riddell [38]). Figure 8 shows the scatters in terms of lnPGA MS vs. lnDI1 MS-AS , lnPGA AS vs. lnDI1 MS-AS , lnPGA1 MS-AS vs. lnDI1 MS-AS , and lnPGA2 MS-AS vs. lnDI2 MS-AS . It is clear that the correlation between the MS (or AS) intensity alone (i.e., PGA MS or PGA AS ) and the structural cumulative damage (i.e., DI1 MS-AS ) is insignificant with the ρ values less than 0.3. Moreover, there is no clear improvement in the correlation between lnPGA1 MS-AS and lnDI1 MS-AS in which the entire earthquake sequence is regarded as one nominal ground motion. Yet, PGA2 MS-AS shows a much stronger interdependence with DI2 MS-AS , with correlation coefficients exceeding 0.5. Besides PGA, similar observation can also be obtained from the other IM-related data, which are summarized in Figure 9. As observed, the correlations between lnIM2 MS-AS and lnDI2 MS-AS for different considered IMs are always much larger than those of the other pairs (e.g., between lnIM MS and lnIM1 MS-AS ). These results confirm that IM2 MS-AS performs much better than the other three types of IMs, i.e., IM M , IM AS , and IM1 MS-AS , for the prediction of the damage potential of earthquake sequences.    According to the results of Figure 10, IA, vrs and PGD are found as the optimal IMs to formulate IM2 MS-AS in the acceleration-, velocity-, and displacement-based IM group, respectively. Then the ρ values between the three pairs, including IA,2 MS-AS vs. DI2 MS-AS , vrs,2 MS-AS vs. DI2 MS-AS , and PGD2 MS-AS vs. DI2 MS-AS , are plotted together at the entire period range, as shown in Figure 11. It is noteworthy that only the results related to the basic SDOF model (i.e., Model-1) are used. It is clear that there is no single IM exhibiting the best correlation throughout the entire period range, and thus, a transition period (i.e., Tc) exists. In fact, before and after Tc, IA and vrs are the most optimal IMs (with the highest correlation to DI2 MS-AS ) to formulate IM2 MS-AS , respectively.  Correlation coefficient (ρ) Periods(s)   According to the results of Figure 10, IA, vrs and PGD are found as the optimal IMs to formulate IM2 MS-AS in the acceleration-, velocity-, and displacement-based IM group, respectively. Then the ρ values between the three pairs, including IA,2 MS-AS vs. DI2 MS-AS , vrs,2 MS-AS vs. DI2 MS-AS , and PGD2 MS-AS vs. DI2 MS-AS , are plotted together at the entire period range, as shown in Figure 11. It is noteworthy that only the results related to the basic SDOF model (i.e., Model-1) are used. It is clear that there is no single IM exhibiting the best correlation throughout the entire period range, and thus, a transition period (i.e., Tc) exists. In fact, before and after Tc, IA and vrs are the most optimal IMs (with the highest correlation to DI2 MS-AS ) to formulate IM2 MS-AS , respectively.  Correlation coefficient (ρ) Periods(s)   Figure 11.
It is noteworthy that only the results related to the basic SDOF model (i.e., Model-1) are used. It is clear that there is no single IM exhibiting the best correlation throughout the entire period range, and thus, a transition period (i.e., T c ) exists. In fact, before and after T c , I A and v rs are the most optimal IMs (with the highest correlation to DI 2 MS-AS ) to formulate IM 2 MS-AS , respectively.

Effect of Hysteretic Models
Two aspects of investigations are conducted on the effect of hysteretic models, including hysteretic model types and R levels. Figure 12 shows the calculated ρ values between the concerned four IM-DI pairs, including lnIM MS vs. lnDI1 MS-AS , lnIM AS vs. lnDI1 MS-AS , lnIM1 MS-AS vs. lnDI1 MS-AS , and lnIM2 MS-AS vs. lnDI2 MS-AS , using all the considered eight types of models with a constant R value of 2. As seen from this figure, the correlation between lnIM2 MS-AS and lnDI2 MS-AS is significantly higher than those of the other three IM-DI pairs regardless of the hysteretic model types. Among the acceleration-, velocity-, and displacement-related IMs, IA, vrs, and PGD are the most optimal IMs to formulate IM2 MS-AS , respectively, due to the associated high correlation with DI2 MS-AS . The observations stated above are consistent with those obtained in Figure 9. Moreover, varying hysteretic models does not cause significant differences in the obtained ρ values. Therefore, the effect of hysteretic model types on the obtained correlation results is limited and has no significant effect over the response for the current purposes.

Effect of Hysteretic Models
Two aspects of investigations are conducted on the effect of hysteretic models, including hysteretic model types and R levels. Figure 12 shows the calculated ρ values between the concerned four IM-DI pairs, including lnIM MS vs. lnDI1 MS-AS , lnIM AS vs. lnDI1 MS-AS , lnIM1 MS-AS vs. lnDI1 MS-AS , and lnIM2 MS-AS vs. lnDI2 MS-AS , using all the considered eight types of models with a constant R value of 2. As seen from this figure, the correlation between lnIM2 MS-AS and lnDI2 MS-AS is significantly higher than those of the other three IM-DI pairs regardless of the hysteretic model types. Among the acceleration-, velocity-, and displacement-related IMs, IA, vrs, and PGD are the most optimal IMs to formulate IM2 MS-AS , respectively, due to the associated high correlation with DI2 MS-AS . The observations stated above are consistent with those obtained in Figure 9. Moreover, varying hysteretic models does not cause significant differences in the obtained ρ values. Therefore, the effect of hysteretic model types on the obtained correlation results is limited and has no significant effect over the response for the current purposes.  Figure 11. The optimal IMs for damage potential prediction of MS-AS sequences at the entire period range.

Effect of Hysteretic Models
Two aspects of investigations are conducted on the effect of hysteretic models, including hysteretic model types and R levels. Figure 12 shows the calculated ρ values between the concerned four IM-DI pairs, including lnIM MS vs. lnDI 1 MS-AS , lnIM AS vs. lnDI 1 MS-AS , lnIM 1 MS-AS vs. lnDI 1 MS-AS , and lnIM 2 MS-AS vs. lnDI 2 MS-AS , using all the considered eight types of models with a constant R value of 2. As seen from this figure, the correlation between lnIM 2 MS-AS and lnDI 2 MS-AS is significantly higher than those of the other three IM-DI pairs regardless of the hysteretic model types.
Among the acceleration-, velocity-, and displacement-related IMs, I A , v rs , and PGD are the most optimal IMs to formulate IM 2 MS-AS , respectively, due to the associated high correlation with DI 2 MS-AS .
The observations stated above are consistent with those obtained in Figure 9. Moreover, varying hysteretic models does not cause significant differences in the obtained ρ values. Therefore, the effect of hysteretic model types on the obtained correlation results is limited and has no significant effect over the response for the current purposes.  To examine the effect of strength reduction ratio, the systems simulated by the basic model with varying R values, i.e., R = 2~5, are adopted for illustration. Figure 13  To examine the effect of strength reduction ratio, the systems simulated by the basic model with varying R values, i.e., R = 2~5, are adopted for illustration. Figure 13   This is because a system with larger R value is easier to suffer nonlinear behavior under a MS excitation, leading to a more complex effect of the AS on the structural damage.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 21 a slightly larger difference on the ρ values than that of varying hysteretic models, especially for the two IM-DI pairs in terms of vrs,2 MS-AS and DI2 MS-AS , and PGD2 MS-AS and DI2 MS-AS . This is because a system with larger R value is easier to suffer nonlinear behavior under a MS excitation, leading to a more complex effect of the AS on the structural damage.

Empirical Prediction Equations
In this section, empirical equations are developed to predict the correlation coefficients, ρ, between the most optimal IM2 MS-AS and the resulting DI2 MS-AS . As seen from Figure 11, IA,2 MS-AS and vrs,2 MS-AS are separately the most optimal IM to formulate IM2 MS-AS before and after the transition period Tc. This conclusion is verified with considering different hysteretic models and varying R levels in this section. Therefore, a step-wise relationship is feasible to represent the relationship between ρ and T. Figure 15 collects the calculated ρ values between the most optimal IM2 MS-AS and DI2 MS-AS corresponding to different hysteretic models with varying R values. From a viewpoint of median prediction, the median ρ values, i.e., ρ50%, at different T values are obtained by statistics; a stepwise equation consisting of a linear and a nonlinear section is regressed by Equation (2). It is clear that Equation (2) can appropriately fit the calculated ρ50% values with a coefficient of determination as R 2 = 0.9.

Empirical Prediction Equations
In this section, empirical equations are developed to predict the correlation coefficients, ρ, between the most optimal IM 2 MS-AS and the resulting DI 2 MS-AS . As seen from Figure 11, I A,2 MS-AS and v rs,2 MS-AS are separately the most optimal IM to formulate IM 2 MS-AS before and after the transition period T c . This conclusion is verified with considering different hysteretic models and varying R levels in this section. Therefore, a step-wise relationship is feasible to represent the relationship between ρ and T. Figure 15 collects the calculated ρ values between the most optimal IM 2 MS-AS and DI 2

MS-AS
corresponding to different hysteretic models with varying R values. From a viewpoint of median prediction, the median ρ values, i.e., ρ 50% , at different T values are obtained by statistics; a stepwise equation consisting of a linear and a nonlinear section is regressed by Equation (2). It is clear that Equation (2) can appropriately fit the calculated ρ 50% values with a coefficient of determination as R 2 = 0.9.
where the transition period T c is taken as 0. 45   The variability of ρ in terms of the standard deviation σρ is further calculated based on the obtained ρ values. Figure 16 shows the σρ values along with the period. It is clear that a stepwise type equation is suitable to fit the σρ values. Then the stepwise function shown in Equation (3) is developed, which shows an acceptable goodness-of-fitting with a value of R 2 equal to 0.84.

Discussions
In this study, a new type of IM for MS-AS sequences, i.e., IM2 MS-AS = IM AS /IM MS , is verified to be an desirable IM to predict MS-AS damage potential due to its high correlation with DI2 MS-AS = DI AS|MS /DI MS . Based on this observation, reliable log-linear relationships can be established between IM2 MS-AS and DI2 MS-AS to predict the structural damage due to MS-AS sequences. Take the cases of T = 0.4 s and T = 3.0 s as examples; both periods are selected according to Equation (2) in the ranges of T ≤ Tc and T < Tc < 6.0 s, respectively. According to the results of Sections 5 and 6, IA,2 MS-AS and vrs,2 MS-AS are proven to be the optimal IMs, respectively. The linear regressions between IA,2 MS-AS and DI2 MS-AS and between vrs,2 MS-AS and DI2 MS-AS are also provided at log-log axis in Figure 17. The variability of ρ in terms of the standard deviation σ ρ is further calculated based on the obtained ρ values. Figure 16 shows the σ ρ values along with the period. It is clear that a stepwise type equation is suitable to fit the σ ρ values. Then the stepwise function shown in Equation (3) is developed, which shows an acceptable goodness-of-fitting with a value of R 2 equal to 0.84. The variability of ρ in terms of the standard deviation σρ is further calculated based on the obtained ρ values. Figure 16 shows the σρ values along with the period. It is clear that a stepwise type equation is suitable to fit the σρ values. Then the stepwise function shown in Equation (3) is developed, which shows an acceptable goodness-of-fitting with a value of R 2 equal to 0.84.

Discussions
In this study, a new type of IM for MS-AS sequences, i.e., IM2 MS-AS = IM AS /IM MS , is verified to be an desirable IM to predict MS-AS damage potential due to its high correlation with DI2 MS-AS = DI AS|MS /DI MS . Based on this observation, reliable log-linear relationships can be established between IM2 MS-AS and DI2 MS-AS to predict the structural damage due to MS-AS sequences. Take the cases of T = 0.4 s and T = 3.0 s as examples; both periods are selected according to Equation (2) in the ranges of T ≤ Tc and T < Tc < 6.0 s, respectively. According to the results of Sections 5 and 6, IA,2 MS-AS and vrs,2 MS-AS are proven to be the optimal IMs, respectively. The linear regressions between IA,2 MS-AS and DI2 MS-AS and between vrs,2 MS-AS and DI2 MS-AS are also provided at log-log axis in Figure 17.

Discussions
In this study, a new type of IM for MS-AS sequences, i.e., IM 2 MS-AS = IM AS /IM MS , is verified to be an desirable IM to predict MS-AS damage potential due to its high correlation with   It is clear that clear linear relationships have been developed in these two figures. Based on the developed log-linear relationships, once the IMs of the mainshock and the following aftershock are determined, the damage ratio DI1 MS-AS can be predicted directly. If the initial damage of structures due to mainshock (i.e., DI MS ) has been estimated, the aftershock-induced damage increment (i.e., DI AS|MS ) and the accumulative damages (i.e., DI1 MS-AS ) due to the MS-AS sequence can be further calculated. Besides, a reliable log-linear relationship between IM2 MS-AS and DI2 MS-AS is also useful for fragility assessment under MS-AS sequences since it commonly assumes a log-linear relationship between earthquake IMs and engineering demand parameters (EDPs). To sum up, there are three folds of applications of the developed log-linear relationship between IM2 MS-AS and DI2 MS-AS , as: (1) to predict DI2 MS-AS given IM2 MS-AS ; (2) to predict DI1 MS-AS and DI AS|MS given IM2 MS-AS and DI MS ; and (3) to represent the relationship between IM and EDP for fragility assessment under MS-AS sequences.

Conclusions
The damage potential for MS-AS sequences was comprehensively conducted by examining the correlation between the intensities of MS-AS sequences and the corresponding structural damage. The following general conclusions have been drawn: (1) The correlations between IM MS and DI1 MS-AS , IM AS and DI1 MS-AS , IM1 MS-AS and DI1 MS-AS , and IM2 MS-AS and DI2 MS-AS were examined in the logarithm space for the SDOF systems with varying periods. The results show that the IM-DI pair in terms of IM2 MS-AS and DI2 MS-AS has the highest correlation among the considered four IM-DI pairs. In other words, the proposed IM, which is defined as IM AS /IM MS , has a better capability to predict damage potential caused by earthquake sequences in comparison with IM MS , IM AS , and IM1 MS-AS . (2) A total of 22 classic IMs were considered as the candidates to define the intensities of MS-AS sequences. Amongst these IMs, IA, vrs, and PGD are the most optimal ones to formulate IM2 MS-AS for the acceleration-related, velocity-related, and displacement-related IM groups, respectively, due to the high correlations with DI2 MS-AS . There is no single IM to best formulate IM2 MS-AS in the entire structural period range. IA,2 MS-AS and vrs,2 MS-AS are the best IMs before and after a transition period Tc, separately. (3) A comprehensive parametric study was conducted by considering various types of hysteretic models and different yield reduction factors. The results show that the effects of varying hysteretic models and yielding reduction factors are limited on the correlation between IM2 MS-AS and DI2 MS-AS . This result reveals that the obtained conclusion in this study is in a general prospect. (4) A step-wise equation consisting of linear and nonlinear parts was regressed as the median prediction of ρ between the best IM2 MS-AS and the caused DI2 MS-AS . Moreover, through regression analysis, another step-wise function was developed to denote the standard deviation of the ρ values along with T. The empirical equations can fit the obtained data reasonably well, providing an opportunity for the further studies to conduct a quick prediction of damage potential for a specific MS-AS earthquake sequence.

Conclusions
The damage potential for MS-AS sequences was comprehensively conducted by examining the correlation between the intensities of MS-AS sequences and the corresponding structural damage. The following general conclusions have been drawn: (2) A total of 22 classic IMs were considered as the candidates to define the intensities of MS-AS sequences. Amongst these IMs, I A , v rs , and PGD are the most optimal ones to formulate IM 2

MS-AS
for the acceleration-related, velocity-related, and displacement-related IM groups, respectively, due to the high correlations with DI 2 MS-AS . There is no single IM to best formulate IM 2 MS-AS in the entire structural period range. I A,2 MS-AS and v rs,2 MS-AS are the best IMs before and after a transition period T c , separately. (3) A comprehensive parametric study was conducted by considering various types of hysteretic models and different yield reduction factors. The results show that the effects of varying hysteretic models and yielding reduction factors are limited on the correlation between IM 2 MS-AS and DI 2 MS-AS . This result reveals that the obtained conclusion in this study is in a general prospect.
(4) A step-wise equation consisting of linear and nonlinear parts was regressed as the median prediction of ρ between the best IM 2 MS-AS and the caused DI 2 MS-AS . Moreover, through regression analysis, another step-wise function was developed to denote the standard deviation of the ρ values along with T. The empirical equations can fit the obtained data reasonably well, providing an opportunity for the further studies to conduct a quick prediction of damage potential for a specific MS-AS earthquake sequence.

Conflicts of Interest:
The authors declare no conflict of interest.