The influence of coupled horizontal–vertical ground excitations on the collapse margins of modern RC-MRFs

With the increasing interest in vertical ground motions, the focus of this study is to investigate the effect of concurrent horizontal–vertical excitations on the seismic response and collapse fragilities of RC buildings designed according to modern seismic codes and located near active faults. It must be stressed that only mid- to high-rise buildings are of significant concern in the context of this research. The considered structures are categorized as intermediate and special RC-MRFs and have been remodeled using distributed and lumped plasticity computational approaches in nonlinear simulation platforms, so that the utilized NL models can simulate all possible modes of deterioration. For better comparison, not only was the combined vertical and horizontal motion applied, but also a single horizontal component was considered for direct evaluation of the effect of the vertical ground motions (VGMs). At the member level, axial force variation and shear failure as the most critical brittle failure mechanisms were studied, while on the global level, adjusted collapse margin ratios (ACMRs) and mean annual frequency of collapse (λCollapse) using a new vector-valued intensity measure were investigated. Findings from the study indicate that VGMs have significant effects on both local and global structural performance and cannot be neglected.


Introduction
Earthquakes in the past have indicated that enormous damage to the building structures and human casualties will result, in the case of severe seismic events. Hence, vulnerability assessment and seismic loss estimation are the primary concerns for regulatory agencies and civil engineers. To realistically asses the structural vulnerability and to incorporate the effects of uncertainties involved in the load-structure system, a probabilistic framework should be utilized. The main components of this framework are presented in Fig. 1.
From a historical point of view, the horizontal component amplitude of ground motions normally plays a dominant role compared to the vertical counterpart. However, acceleration records from the (1989) Loma Prieta earthquake and the (1994) Northridge earthquake in the USA, the (1995) Hyogoken-Kobe earthquake in Japan, (2003) Bam earthquake in Iran, and the (2011) Christchurch earthquake in New Zealand, among others showed that the magnitudes of the vertical component can be as large as, or exceed, the horizontal component. The report from Elnashai et al. (1995) also highlighted cases of brittle failure induced by direct compression, or by reduction in shear strength and ductility due to variation in axial forces arising from the vertical motion in the (1994) Northridge earthquake. In such situations, most existing code specifications assume that the ratio of vertical component of the ground motion to that of the horizontal component (V/ H) varies from 1/2 to 2/3, which must be considered unconservative and needs to be investigated.
Recent studies (Bozorgnia and Campbell 2004;Elgamal and He 2004) on horizontal and vertical ground motions have indicated that such a simple approach is not valid and appropriate for the near-fault regions anymore. The main reasons can be categorized as follows: • The attenuation rate for vertical ground motion is much higher than that of the horizontal ground motion. This rate increase in the far-field areas. Thus, structures built in the near-fault regions experience higher vertical excitations. • Vertical ground motion includes more high-frequency content than horizontal ground motion. The difference increases with the decrease in the soil stiffness.
It should also be noted that the higher values of V/H ratio do not necessarily imply more energy content on the desired structure. The reason is that the two components may not coincide in time to cause strong interaction effects.
Besides these, many of the current seismic design codes and damage estimation tools do not include the effect of vertical ground motions on the seismic response of structures and especially columns. However, the observed damage on the columns (diagonal shear cracks) during historical seismic events such as the 1994 Northridge earthquake and the 1995 Kobe earthquake was partly attributed to the effect of vertical motions (Broderick et al. 1994;Elnashai et al. 1995). Field and analytical evidence by Papazoglou and Elnashai (1996) indicated that strong vertical earthquakes can cause a significant fluctuation in the axial force in columns, resulting in a reduction in their shear capacity and compression failure of some of the columns. During the 1995 Kobe earthquake in Japan, the RC structures exhibited very high amplifications of the vertical component of more than two times. The main reasons were the low damping mechanism in the vertical direction and the absence of supplement seismic energy dissipating systems in this direction. On the other hand, because of the high stiffness in the vertical direction, a quasi-resonant response was observed in these structures. High-frequency pulses from vertical motion were recognized as the other reason for such a phenomena (AIJ 1995). Iyengar and Shinozuka (1972) investigated the effect of self-weight and vertical accelerations on the behavior of tall structures. The structures have been idealized as cantilevers and the ground motion as a random process. The main conclusion from their study was that inclusion of self-weight simultaneously with the vertical ground acceleration can increase or decrease the global peak response. These fluctuations in structural response had been considerable in most cases. On the local level, beams were identified as the most critical elements and the effect of vertical ground motion on them had been pronounced. Iyengar and Sahia (1977) investigated the effect of vertical ground motion on the response of cantilever structures using the mode superposition method; their main conclusion is that the consideration of the vertical component is essential in analyzing towers. Anderson and Bertero (1977) used numerical methods to evaluate the inelastic response of a ten-story unbraced steel frame subjected to a horizontal component of earthquake and to combinations of this component with the vertical one; they deduced the following points; the inclusion of the vertical motion on one hand does not increase the displacements, but on the other increases the girder ductility requirement by 50 % and induces plastic deformations in columns. Mostaghel (1974) and Ahmadi (1978Ahmadi ( , 1980 studied the effect of vertical motion on columns and tall buildings which have been idealized as cantilevers, using the mathematical theory of stability of Liapunov. Their main conclusion was that, in the inelastic region, if the maximum applied earthquake loading would be less than the Euler buckling load, it is guaranteed that the column will remain stable irrespective to the type of earthquake loading, and the inclusion of vertical ground excitation can be neglected. But this is unlikely to be the case for reinforced concrete columns because of the crushing of concrete in compression and the buckling of the yielded reinforcement. Munshi and Ghosh (1998) investigated the seismic performance of a 12-story RC building under a combination of horizontal and vertical ground motions. This analysis showed a slight increase in the maximum deformation when the vertical ground motion was included. The formation patterns revealed that vertical accelerations induced a slightly different hinge formation pattern and hinge rotation magnitude, and the response of the frame-wall system did not show sensitivity to the vertical acceleration in this case. Antoniou (1997) studied the effect of vertical accelerations on RC buildings by analyzing a eight-story reinforced concrete building designed for high ductility class in Euro code (EC8) with a design acceleration of 0.3 g. This analysis showed that the vertical ground motion can increase the compressive forces by 100 % or even more and lead to the development of tensile forces in columns. These fluctuations in axial forces can result in shear failure in these elements. Kim et al. (2011) studied the effect of various peak V/ H ground acceleration ratios and the time lag between the arrival of the peak horizontal and vertical accelerations on the inelastic vibration period and column response for infrastructures. It was observed that the inclusion of vertical motions notably influenced the inelastic response vibration periods and considerably increased or decreased the lateral displacement. It was also noticed that the arrival time had a minimal effect on the axial force variation and shear demand.
None of the previous studies have investigated the codeconforming RC-MRFs utilizing fragility curves and reliability methods. As the seismic vulnerability assessment of high-rise structures is a complex task, it is important to consider that both the lower and higher structural modes might be excited, because of the wide range of frequency content of the applied earthquake loads. On the other hand, the imposed displacements to these structures can be very significant, since the fundamental period of many high-rise  Int J Adv Struct Eng (2016) 8:169-192 171 structures are within the period range of 1-5 s, which corresponds to the peak displacement spectra of the standard earthquakes. To this end, in the current study, both distributed fiber-based and lumped plasticity approaches and various modes of collapse are considered in the simulation process. To show the significance of vertical ground excitations and to get the most accurate results, a new definition for V/H ratio and an optimum intensity measure are proposed. Various mass distributions are considered in the eigenvalue analysis to determine the most accurate and computationally efficient structural model. Seismic fragility curves as the main component of the current study can be derived using various approaches: observational, experimental, analytical and hybrid techniques to quantify damage and estimate monetary losses (Calvi et al. 2006). While the observational method is the most realistic and rational one, as the entire inventory is taken into consideration it is usually difficult to be utilized because no or insufficient observational-based date are available from the past events. The experimental method is not a feasible option in many cases, because of its cost and the time needed, since a wide range of structures should be tested. In the current study, the third approach based on extensive nonlinear analytical simulations is adopted. This option is the most feasible and possible methodology which can be used in many cases.
The main objective of this study is to calculate the collapse margins and mean annual frequency of collapse as the performance metrics employing displacement-based fragility curves for multiple limit states from concrete cracking to structural collapse in the near-fault areas. The collapse of structures is determined on the basis of the global failure mechanism of the structural system rather than the failure of a structural element. To achieve this goal, numerical models that capture the axial-shear-flexural behavior of the columns are created in nonlinear seismic simulation platforms, Zeus-NL and OpenSees (Elnashai et al. 2004;McKenna 2014).

Selection and characterization of input ground motions
A major stage in the process of fragility estimation is to use appropriate ground motion (GM) records. If the ground motion selection would be done in a way that the hazard consistency is ensured, then the results from the corresponding simulation and analysis would be rational. According to the recommendations of FEMA P-695 (2009), 40 earthquake ground motion records at varying hazard levels from 20 earthquake events are selected from the PEER NGA-WEST2 (Ancheta et al. 2012(Ancheta et al. , 2013Pacific Earthquake Engineering Research Center 2015) database (Table 1). This database is the most recent and very suitable for adequate fragility analysis.
The criteria for selection of the analyses records include medium to high vertical component, having large magnitude (M w C 6.0) and recorded at near-fault rupture distances (R B 25 km), with a frequency range to excite the periods of vibration of the structure in both horizontal and vertical directions. The response spectra of the selected records, the median of the acceleration response spectra and the dominant periods (T 1 ) of the reference structures which will be defined afterward are shown in Fig. 2.
Based on the results illustrated in Fig. 2, the vertical component of a ground motion tends to concentrate all its energy content in a narrow, high-frequency band, while the frequency range for the horizontal component is much wider. Hence, this phenomenon will amplify the structural responses in the short period range, which usually coincide with the vertical periods of RC elements/structures. After the GMs selection, they are amplitude scaled, using the procedure outlined in ASCE 7 (ASCE 2010) to match the 5 % damped site-specific target spectrum (corresponding to the maximum credible earthquake, MCE) given in Fig. 3.

Representative set of structures
Four RC-MRFs ranging from 7 to 20 stories are selected to represent medium-and high-rise buildings. The frames are designed and detailed according to ACI building code ACI-318 (2011) and ASCE 7 (2010) provisions. Two categories of RC-MRFs, special and intermediate, are used in the current study. The ordinary MRFs, because of their low level of ductility during an earthquake, are not considered here. The special MRF employs the strong column weak beam (SCWB) concept and specifies elaborate detailing of joints. Thus, the SMRF is expected to form the sway mechanism and possesses a high degree of ductility. On the other hand, the intermediate MRF has enough strength as well as reasonable ductility and can be used throughout most of the seismic-prone areas. 7-and 12-story buildings are designed as intermediate MRFs, while 15-and 20-story buildings are designed as special MRFs. The behavior factors (R) are considered as 5 and 8, respectively (ASCE 7 2010).
Span lengths are identical in both directions equal to 6 m, and story heights are 3.5 m. Damping is set to be 5 % in the first three modes and is considered as of Rayleigh mass and stiffness proportional type, based on the recommendations of Zareian and Medina (2010). The reference structures are shown in Fig. 4 and the beam and column dimensions and longitudinal and transverse steel ratios are listed in Table 2.

Structural modeling approaches
In general, most of the current models provided by other researchers cannot be utilized to predict the accurate behavior of the RC elements in the presence of vertical excitations. The main reason is that such models do not account for important response features, such as the interaction between shear, flexure and axial forces. One of the main failure modes of RC columns due to vertical excitations as mentioned previously is the shear failure in these members. To include this type of failure in the analytical models, some modification should be incorporated in the modeling approaches. To this end, two nonlinear simulation platforms; ZEUS-NL (Elnashai et al. 2004) and OpenSees (McKenna 2014), are used to simulate the seismic response of the reference structures. Fixed-base models are used in the analysis stage; as a result, soil-structure foundation interaction is neglected. A leaning column to account for the P-D effect from loads on the gravity system is also considered in both platforms.
In ZEUS-NL, the columns are modeled with distributed nonlinear fiber sections (flexural response) and an NL zerolength shear spring (shear response) at the end of each column. The idealization adopted in the first approach effectively models reinforcing steel, and unconfined and confined concrete (Mander et al. 1988). This approach allows monitoring the stress-strain response at each fiber over several Gauss sections through the integration of the nonlinear stress-strain response of different fibers in which the section is subdivided, as shown in Fig. 5 (shorter fiberbased elements at the end of the column and longer fiberbased elements away from the end of the column are used to capture inelastic flexural response in plastic hinge zones). This modeling approach reduces the modeling uncertainty. Material properties used from the large database of ZEUS-NL for including the damage plasticity in both concrete and reinforcements. ''stl1'' material model is chosen for reinforcement steel, which is a bilinear elasto- ''con2'' material model is chosen for the concrete, which is a uniaxial concrete model that includes the confinement provided by hoops or ties.
To simulate the shear response of the columns Elnashai (2001, 2002) developed two shear models: hysteretic shear model under constant axial force (Lee and Elnashai 2001) and hysteretic shear model under axial force variation (Lee and Elnashai 2002). In case of vertical excitations, as the axial force is not constant, the first shear model cannot capture the effect of axial force variation on the seismic response. Figure 6a shows the envelope curve of the hysteretic shear model under constant axial force, which is defined by a quadrilinear symmetric relationship comprising cracking, yielding and ultimate conditions. The response parameters on the curve are determined by the modified compression field theory (MCFT) developed by Vecchio and Collins (1986). To address the effect of varying axial force on the envelope curve, the basic formulation can be extended and the curves extracted for multiple levels of axial force, as shown in Fig. 6b. 3 @ 6m c/c 7 @ 3.5m c/c 3 @ 6m c/c 4 @ 6m c/c 4 @ 6m c/c 12 @ 3.5m c/c 15 @ 3.5m c/c 20 @ 3.5m c/c A fully concentrated plasticity (Krawinkler et al. 2004) approach is utilized in OpenSees to consider both the flexural and shear behavior of the columns in terms of vertical excitations. To this end, all the beams, columns and joints are modeled with NL springs. Beams and joints are modeled with rotational springs and the columns are modeled using rotational springs along with a zero-length shear spring located at one end of the column to consider the effect of shear failure as shown in Fig. 7. The plastic hinge behavior of the beam and column elements are calibrated through the large set of experimental data by Haselton and Deierlein (2007) and Haselton et al. (2008).
The hysteretic uniaxial material model monitors the shear/flexural behavior of the columns. In this model, shear deformations are simulated using a shear spring, while the flexural deformations are monitored using beam/column elements. As illustrated in Fig. 8, the solid line indicates the total response of the RC column. This line captures the degrading shear behavior, when the shear strength would be less than the shear corresponding to the plastic hinges development. A dashed line is presented in Fig. 8, which shows the case in which the shear strength is higher than the shear corresponding to the plastic hinges development and the model does not capture any shear degradation.  Watanabe and Ichinose (1992), Aschheim and Moehle (1992), Priestley et al. (1994) and Sezen (2002) have shown that the shear strength in RC elements decays with increased plastic deformations. Hence, the dashed line provided in Fig. 8 cannot be realistic and accurate, if the column yields in flexure close to its estimated shear strength. The main deficiency of the hysteretic uniaxial model is that it determines the point of shear failure based only on the column shear, while it should be determined based on both force and deformation.
To resolve this problem, the shear load versus deformation model for the shear spring was developed using the  existing OpenSees limit-state material model and shear limit curve developed by Elwood (2004). As shown in Fig. 8, the shear limit curve is activated and shear failure is initiated once the column shear demand exceeds the column shear capacity. In this case, the limit-state material model simulates and captures the RC column response to detect the possibility of shear failure. To this end, the shear limit curve is defined according to both the column shear and the total displacement or drift ratio (Fig. 8). In case the columns are vulnerable to shear failure after flexural yielding, then the drift capacity model proposed by Elwood and Moehle (2005) can be utilized to define the accurate limit curve. Other shear failure criterion such as plastic rotation at the two ends of the column is also proposed by LeBorgne (2012). The beams' and columns' moment-rotation behavior is simulated utilizing an NL hinge model including strength and stiffness deterioration developed by Ibarra et al. (2005) and has been used in PEER/ATC (2010) as well. As the reference structures are code conforming and have ductile behavior, the peak-oriented model is utilized. This model combines a post-peak negative stiffness branch of the backbone curve to capture in-cycle strain softening and a cyclic model to capture the strength and stiffness deterioration based on the cumulative hysteretic energy dissipated (Fig. 9). The individual modes of deterioration are described briefly below: Basic strength degradation This mode of deterioration is defined in a way to show a reduction in the yield strength. This mode also includes the strain hardening slope degradation. Based on the formulation, the values for positive and negative directions are defined independently (Fig. 9a).
Post-capping strength degradation It is defined by translating the post-capping branch toward the origin; however, the branch slope will keep constant and it will move inward to reduce the reference strength. The process will be done each time the X axis is crossed (Fig. 9b).
Unloading stiffness degradation This mode of deterioration indicates the reduction in both positive and negative unloading stiffness. The unloading stiffness (K u,i ) in each cycle depends on the unloading stiffness in the previous excursion (K u,i-1 ) (Fig. 9c).
Accelerated reloading stiffness degradation This mode of deterioration escalates the target displacement according to the loading trend in both positive and negative directions (Fig. 9d).
For ductile RC frames, beam-column joints are often modeled with rigid joint zones. However, Shin and LaFave (2004) argued that the joint regions are not fully rigid and may experience some shear deformations that can contribute to global deformation. Therefore, the analytical model must predict the inelastic behavior of ductile beam-column joints. In the current study, the pinching material model is utilized for the joint rotational springs. The modeling parameters of the pinching material for ductile exterior and interior beamcolumn joints have been calibrated by Jeon (2013) and used in the modeling procedure (Fig. 10).

Eigenvalue analysis
The detailed analytical models were subjected to eigenvalue analyses to determine the fundamental period of the structures (Chopra 1995). Zero or small errors for the horizontal modes imply that distributing the lumped-mass mesh over beams constitutes no differences, but for the case of vertical excitation the results show a lot of discrepancies. The vertical and horizontal periods for several lumped-mass models and a distributed one are calculated and compared ( Fig. 11; Table 3). It is then assumed that if a simplified lumped-mass model gives very similar natural mode shapes and vibration periods compared to the distributed one, this simplified model is reliable to simulate vertical motions and can be used in the NL-RHA.
Compared to L-Mass 1, L-Mass 2 is much more accurate. However, there are clearly differences in the vertical periods. L-Mass 3 and L-Mass 4 show both significantly similar periods in the horizontal and vertical directions. Even though L-Mass 3 has a coarse lumped-mass approach compared to the distributed mass model, the eigenvalues of this model imply very similar tendencies to the exact solution. As a result, because of the small differences in L-Mass 3 compared to L-Mass 4 results, L-Mass 3 is the most simplified lumped-mass model to cover realistic vertical motion with minimum computational effort and is implemented in the NL-RHA of the studied structures to provide fragility curves.

V/H ratio (conventional and proposed approaches)
The current and most common design practice is to take the V/H spectral ratio as 2/3 as proposed by Newmark et al. (1973) and is also used by FEMA (Bozorgnia and Campbell 2004). However, this approach is inaccurate for nearsource moderate and large earthquakes (Friedland et al. 1997;Bozorgnia et al. 1999;Bozorgnia and Campbell 2004;Button et al. 2002). To clarify this issue, a large set of earthquake records were extracted from NGA-WEST2 (Ancheta et al. 2012(Ancheta et al. , 2013Pacific Earthquake Engineering Research Center 2015) dataset to compare the peak vertical to horizontal ground acceleration (V/H ratio) at different ground motion magnitudes and distances. Figure 12 shows the detailed comparison. For a better illustration, a linear trend line is fitted to all data points.
The results indicate that the 2/3 rule is unreasonable and confirms that the V/H ratio may show significant variations, which depend on the source and site characteristics and seismic radiation pattern. The shortage of the conventional V/H calculation methods is that they do not differentiate the horizontal and vertical frequency content and cannot include the influence of dynamic properties of the considered structure. To this end, a new approach is proposed by the authors, which calculates the V/H based on the dominant vertical and horizontal spectral acceleration of the studied structures. NGA-WEST2 (Ancheta et al. 2012(Ancheta et al. , 2013; Pacific Earthquake Engineering Research Center 2015) dataset composed of more than 20,000 records is used to provide reliable calculations. Based on the eigenvalue analyses, the horizontal and vertical dominant periods of the case studies were evaluated and used to calculate the vertical to horizontal ratio for each structure exclusively. The benefit of the new approach compared to the conventional methods is that for each structure with a specific structural system, the V/H will be calculated based on its dynamic properties and would be unique. The results for the new approach are presented in Fig. 13. Red lines indicate the central 50 % (median), while lower and upper boundary lines show the 10 and 90 % quantile of the data. Two vertical lines extend from the central box, indicating that the data remaining outside the central box extends maximally to 1.5 times the height of the central box, but not past the range of the data. The reaming data points plotted by red markers are the outliers. The given results show that there is a correlation between the structure's height and V/H ratio, and the high-rise RC-MRFs experience higher levels of seismic vertical excitations. However, it does not prove that these structures may be more

Performance criteria
In this study, interstory drift ratio (IDR) is considered as the engineering demand parameter (EDP). This is particularly a suitable choice for RC-MRFs, since it relates the global response of the structure to joint rotations where most of the inelastic behavior in the moment-resisting frames is concentrated. At the first stage, limit states are considered corresponding to different performance levels as specified in FEMA 356 (2000). The IDR values of 1, 2 and 4 % are used for IO, LS and CP performance levels, respectively. In the next stage, these performance levels are calibrated and verified for each reference structure, based on its structural capacity through nonlinear pushover analysis. The IDR limit for each individual structure from pushover analysis combining of the first four natural modal shapes, weighted by modal mass participation factors, is used for the lateral load distribution pattern. IDR values corresponding to the first concrete cracking/steel yielding, concrete strain corresponding to maximum confined concrete stress and maximum confined concrete strain are considered as IO, LS and CP performance levels, respectively. Table 4 lists the median capacity values against each performance level for all buildings. The results show good agreement of the NL-pushover capacities with FEMA-356 performance levels.
Besides the limit states mentioned above, the collapse limit state can be defined on the basis of a different approach. The onset of 'collapse' for a ground motion record is identified as the point where maximum IDR response increases 'drastically' when the spectral acceleration of the record is increased by a 'small' amount (Vamvatsikos and Cornell 2002). To this end, in the next stages, incremental dynamic analysis (IDA) is performed on the reference MRFs and the collapse is defined as the point of dynamic instability when IDR increases without bounds for a small increase in the ground motion intensity for each structure individually.

Proposing an optimum intensity measure for fragility analysis
The ground motions are characterized by intensity measures (IMs), and their choice plays a crucial role in the seismic fragility estimation. An optimal IM is the one that has good efficiency, sufficiency, practicality, hazard computability and predictability among other characteristics (Mackie and Stojadinovic 2001;Luco and Cornell 2007;Giovenale et al. 2004;Padgett et al. 2008). Efficiency means the ability to accurately predict the response of a  structure subjected to earthquakes (i.e., small dispersion of structural response subjected to earthquake ground motions for a given IM). A sufficient IM is defined as one that renders structural responses subjected to earthquake ground motions for a given IM conditionally, independent of other ground motion properties (i.e., no other ground motion information is needed to characterize the structural response). Previous studies have shown that PGA is not an accurate and ideal IM for evaluating the geotechnical phenomenon, as it cannot consider the ground motion duration (Kramer et al. 2008). It is found that the fragility curves based on vector-valued IM are better able to represent the damage potential of earthquake (Baker 2015). Thus, an optimized vector-valued intensity measure, which includes the geometric mean of spectral accelerations over a range of period, is considered in the current study. The parameters in the vector-valued intensity measure should be chosen to convey the most possible information between the ground motion hazard and the structural response stages of analysis. This requires identifying parameters that most affect the structure under consideration. The adopted IM is a suitable choice, especially for the high-rise RC-MRFs where the effect of higher modes is significant. Given that S a (T 1 ) has been verified as an effective predictor of structural response for a wide class of structures, it will be used as the first element of the vector, while the effect of higher modes (HM) is considered as the second parameter (Eqs. 1 and 2).  where the constant k 1 is chosen to vary between T low /T 1 and 1, and k u between 1 and T upp /T 1 , where T low and T upp are, respectively, the lower and the upper period of the elastic spectrum. It is worth mentioning that k 1 and k u values can be optimized for each structure separately based on its dynamic characteristics using a trial and error procedure. An efficient IM reduces the amount of variation in the estimated demand for a given IM value (Giovenale et al. 2004;Padgett et al. 2008). Employing an efficient IM yields less dispersion about the estimated median in the results of the NL-RHA. The comparison of results in terms of standard deviation for conventional scalar and vector IMs are illustrated in Fig. 14. In the comparison, S a (T 1 ) is considered as the most conventional scalar IM, while [S a (T 1 ), S a (T 2 )/S a (T 1 )] is considered as a conventional vector IM.
As can be seen from Fig. 14, S a (T 1 ) can be an efficient predictor for structures with short period, but not for the structures with medium to large periods. On the other hand, the vector IM can predict the results in an efficient way for both short and large period structures; however, the results have more standard deviations compared to the proposed IM. In all cases, IM (New) has been more efficient and sufficient, as it is always associated with small values of dispersion. The maximum dispersion values are 20, 31 and 50 % for IM (New), IM (Vector) and S a (T 1 ), respectively.

Probabilistic demand models
To formulate and correlate the earthquake intensity measure (IM) to the building-specific demand measure (DM), a probabilistic seismic demand model (PSDM) should be developed. Selecting an appropriate pair of IM and DM is a difficult task, as the PSDM results are practically sufficient and efficient (Gardoni et al. 2002(Gardoni et al. , 2003. In the current study, response history analyses (RHAs) are used to quantify the seismic demand of the reference structures. The maximum interstory drift ratio (Max IDR) is considered as a suitable DM, since it is closely related to damage states of the RC-MRFs. The final results are illustrated in terms of IDR-IM in Fig. 15, under 40 sets of ground motions chosen in the previous sections.

NL-incremental dynamic analyses and seismic fragility estimate
In the fragility estimation process, a suitable analysis procedure should be implemented. Based on the recommendations provided by Vamvatsikos and Cornell (2002) and FEMA P-695 (2009), incremental dynamic analysis (IDA) is recognized as one of the best and most common procedure, where a suite of earthquake records is scaled repeatedly to find the intensity measure in which the structural collapse will occur. In this study to isolate the effect of each component, horizontal and vertical ground motions were applied separately and simultaneously.
Using the IDA approach, information about variability in ground motions can be directly incorporated into the collapse performance assessment. However, this process only captures the record-to-record variability and does not account for how well the nonlinear simulation model represents the collapse performance of the reference structures; hence, model uncertainties should also be accounted in the collapse simulation, which will be discussed later.

Effect of vertical excitations on the structural responses (local level)
Seismic performance evaluation of RC structures subjected to (medium to strong) multi-component ground motions has some complexities. One of these difficulties is the increment in axial force variation in RC columns, which can be superimposed in the overturning forces. Since there is a direct relation between flexural, shear and axial forces, the fluctuation in axial force increases the possibility of shear failure. This is mainly due to the significant variation in strength and stiffness in the columns caused by vertical excitation. NL-RHA was performed using the selected natural horizontal ground motions applied with and without the vertical component. Seismic loading may be applied upward as well as downward, thus subjecting structural members to unaccounted-for action, if the naturally existing vertical ground excitations are neglected in the design procedures. Among the structural members, columns are the most vulnerable and may be adversely affected, if high compressive forces are developed or if axial forces change to tension. Sample results of the effect of vertical ground motions on the first-story column of the 7S3B reference model are shown in Fig. 16. Table 5 provides the ratio of the vertical seismic force (maximum axial force in the column) to axial gravity load on the column with and without vertical ground motion for the reference structures for three V/H ratios in the firststory interior columns. From the obtained results, it is observed that the effect of vertical ground motion on the ratio of axial force to gravity load increases significantly for all the models in the range of 4.19 and 108.59 % with increase of the V/H ratio.
Another complexity which occurs in the presence of high vertical forces is the increment in ductility demand and reduction in ductile capacity; which may result in extensive damages. The shear capacity may be significantly affected by the variation in column axial forces, which may cause loss or reduction of the axial load contribution to shear strength. As shown in Fig. 17, the shear capacity of an interior first-story column in 7S3B frame was calculated using ACI-318 (2011) and ASCE-41 (2013) formulas and compared to shear demand with and without vertical ground motion. The priority of ASCE-41 (2013) formula over the ACI-318 (2011) formula is that it is based on the displacement ductility demand, which is more realistic. It is very evident that in the case of VGM, the shear demand exceeds the shear capacity. Hence, the possibility of shear failure would be increased before the structure reaches the global collapse state. It is worth mentioning that similar results are observed for all the studied MRFs.

Seismic fragility estimation (global level)
For the case of earthquake excitations, a closed-form solution is not usually available, since there are a large number of random variables and different probability density functions associated with these events. To resolve this problem, the reliability of the structures under these complex phenomena can be represented using a probabilistic methodology incorporating fragility curves. These curves define the probability of exceeding a specific damage state subjected to a hazard by a suitable IM. In this study, following the structural reliability theory by Ditlevsen and Madsen (1996), the fragility functions are formulated as in Eq. (3): Based on the inherent randomness and uncertainty in the capacity (C), the seismic demand (D) and the limit states, a closed-form approximation using a log-normal distribution was proposed by Wen et al. (2004). This fragility formulation is given as: where U(.) denotes the standard normal CDF, and v C and v D|IM are the natural logarithm of the median capacity and demand of the system, respectively, while b C , b D|IM and b M represent the uncertainty in estimating the capacity, demand and structural modeling. Fragility estimates for the studied MRFs are obtained using the IDA results for the optimized HM values and are presented in Figs. 18 and 19. Both modeling approaches are included in the fragility estimation and the results show that two approaches match very well and less than 10 % dissimilarity was observed in all performance levels. For a better demonstration of the collapse mode under concurrent horizontal-vertical excitations, the collapse fragility surfaces are plotted in a 3D space as well (Fig. 20).
Based on Figs. 18 and 19, when the vertical component is coupled with the horizontal excitations, the fragility results can be changed extensively. The maximum increase in the structure's fragility appeared in the CP and Collapse damage modes. Figure 20 shows the developed fragility surface based on a vector-valued IM and can be visualized as fragility curves by projecting the surface onto the planes. These figures demonstrate the wide variation between fragility curves based on scalar-valued intensity measure.
It can be seen in Fig. 20 that there is a discrepancy of up to 30 % between the curves calculated for various HM values. The main advantage of these fragility surfaces is that the variability of structural fragility due to a second parameter can be accounted for in contrast to when fragility curves are used. This means that which records should be used depends on the seismic hazard at the site when scalarvalued IM [S a (T 1 )] is used to evaluate the structural fragility.  Adv Struct Eng (2016) 8:169-192 185 Ignoring the effect of HM will bias the final results. For example, if the seismic hazard disaggregation suggests that extreme motions are associated with records having a mean value of HP of about 0.75, but records are selected with a mean value of about 0.50, then the S a (T 1 )-based result will underestimate the seismic fragility. In other words, the evaluation of structural fragility by means of vector-valued IMs reduces the complexity of record selection procedure based on the seismic environment (i.e., magnitude, distance, site conditions, etc.).

Collapse performance evaluation
Prior to the development of incremental dynamic analysis (IDA) and its use in the FEMA P-695 methodology (2009), accurate modeling of buildings near collapse, in the negative post-peak response range, was not a high priority of research. In recent FEMA guidelines (e.g., FEMA P-440A 2009; FEMA P-695 2009), sideway collapse (where collapse is defined based on unrestrained lateral deformations with an increase in ground motion intensity) is typically assumed to be the governing collapse mechanism.
In the current study, important metrics for quantifying collapse resistance of structures are defined in the previous section and illustrated in Fig. 21, including the median collapse capacity and the conditional probability of collapse at an intensity level of interest, the code-defined maximum considered earthquake (MCE). As the differences in two modeling approaches were negligible, results from the ZEUS-NL models which had higher accuracies are utilized in this section.
The MCE intensities are obtained from the response spectrum of MCE ground motions (Fig. 3) at the fundamental period (T 1 ). The ratio between the median collapse intensity and the MCE intensity is defined as the collapse margin ratio (CMR), which is the primary parameter used to characterize the collapse safety of the structure (Eq. 5): Comparing the intermediate and special MRFs located in the same hazard region, Fig. 21 indicates that the 20and 15-story buildings have better collapse behavior compared to the 7-and 12-story buildings. The reason for this behavior is due to the higher ductility level in SMRFs. As discussed in Haselton et al. (2008) and FEMA-P695 (2009), the collapse capacity of structures and the calculated CMR could be strongly influenced by the frequency content (spectral shape) of the ground motions. Thus, in FEMA-P695, a spectral shape factor (SSF) is proposed to modify the CMR values. To this end, adjusted CMR (ACMR) value for each structure is assessed (Eq. 6): The SSF value for each structure is extracted based on the fundamental period (T 1 ), period-based ductility (l T ) and the seismic design category from FEMA P-695 documentation. For seismic performance evaluation of the studied structures under concurrent horizontal-vertical excitations, on the basis of the 5 % probability of collapse under MCE and the composite uncertainty (b TOT = 0.5) in collapse capacity, the acceptable ACMR for all building models takes the value of 2.28 according to FEMA-P695 recommendation. The final acceptance criteria for the reference structures are reported in Table 6. The results show the CMR, SSF and ACMR values for each archetypical design. Later on, the acceptable ACMRs are compared with the calculated values to see whether the structures could either pass or fail the criterion.
Focusing on the four reference structures, Table 6 shows that SMRFs have acceptable ACMR, while a disturbing trend becomes evident for the IMRF buildings. The results show that in terms of coupled horizontal-vertical excitation, both 7S3B and 12S3B buildings have unacceptable ACMR, and surprisingly 7S3B would also fail for the horizontal excitation case while the 12-story frame passed the criteria marginally. As a result, the IMRFs do not attain the collapse performance required by FEMA P-695 methodology, and additional design requirement adjustments would be needed to improve the overall performance. It means that even the code-conforming design structures with acceptable level of ductility can be vulnerable to seismic excitations.
Comparison of the calculated ACMRs in terms of (H) and (H ? V) excitations demonstrates that generally in all models, the safety margin against collapse reduces by including the vertical component of earthquake and the reductions are very remarkable and pronounced. Fortunately, the SMRF models could pass the collapse performance criteria for both types of excitation. It shows that their elements, being controlled by many detailing and capacity design requirements of the building code, limit possible failure modes.

Mean annual frequency (MAF) of collapse
Ground motion hazard curves for a typical highly seismic and populated region in the Middle-East is Fig. 19 Fragility curves for special RC-MRFs at the optimal HM Int J Adv Struct Eng (2016) 8: 169-192 187 illustrated in Fig. 22 at the periods of 0.2, 1 and 3 s, and the required data at the fundamental period of structures are interpolated from these values. The buildings are assumed to be located at a site in Tehran, for which the hazard curve has been defined through probabilistic seismic hazard analysis (PSHA) by Gholipour et al. (2008). Another collapse metric, called the mean annual frequency (MAF) of collapse, is defined by integrating the seismic hazard curve and the fragility results. This metric describes how likely it is for collapses to occur, considering both the structural collapse capacity and the ground-shaking hazard in the studied region (Krawinkler et al. 2004;Liel 2008). The mean annual frequency of collapse (k Collapse ) is computed using Eq. (7) (Haselton and Deierlein 2007): where P[Sa Collapse B x] is the probability that x exceeds the collapse capacity (i.e., the probability that the building collapses when the ground motion intensity is x), and k IM (x) is the mean annual frequency of the ground motion intensity exceeding x (i.e., a point on the seismic hazard curve). There are many ways to approximate Eq. (7). A closed-form solution is used to fit an exponential function to the hazard curve. To avoid error induced by fitting an exponential function, the PCHIP (piecewise cubic hermite interpolating polynomial) procedure is used to interpolate between the points (MathWorks 2015). In the next stage, the interpolated curves are implemented to complete the numerical integration required to evaluate Eq. (7) (Fig. 23). Based on the extracted results in Fig. 23, the k Collapse(max) in the intermediate MRFs are 2.1 9 10 -5 and 1.86 9 10 -4 collapse/year for the horizontal and combined H ? V, respectively. However, these values are relatively smaller for the special MRFs under both types of excitations. The maximum values for the SMRFs are 2.3 9 10 -6 and 4.5 9 10 -5 collapse/year for the horizontal and combined H ? V, respectively. It is worth mentioning that these results correspond to collapse return periods of 2475 years. Based on MAF calculation, SMRFs are in a higher confidence bounds of safety compared to the intermediate moment frames, while the effect of vertical ground motion is very significant for both groups and cannot be neglected.

Conclusions and recommendations
In this study, the effect of vertical excitations on the seismic performance of intermediate and special RC-MRFs has been evaluated. It is important to emphasize that the main objective of this study was not to quantify numerical design values. The objective was rather to focus on the importance or otherwise of including vertical ground motion in design of RC buildings and its impact on the member and the structure levels. Hence, a large number of NL dynamic analyses were performed using fiber-based and concentrated plasticity approaches. The computational models were utilized in ZEUS-NL and OpenSees platforms to compare the results. The VGM was shown to be significant and should be included in the analysis when the proposed structure is located within 25 km of a seismic source.
The most important findings are summarized as follows: • The vertical component of an earthquake tends to concentrate all its energy content in a narrow band, unlike the horizontal counterpart. This energy  • Based on the frequency content of the vertical ground motion, it can be concluded that structural failure modes from past earthquakes might be attributed to underestimating the effect of vertical acceleration in the design procedures, and there is an urgent need for the adoption of more realistic vertical spectra in future version of seismic design codes. • Although the effect of axial force on shear capacity of the structural elements is an accepted fact and is proved in the current study, current seismic codes do not have a consensus on this effect, and different code equations might lead to different shear capacity estimations. Both the ACI-318 and ASCE-41 equations captured the shear strength degradation due to axial force. ASCE equation predictions could be considered as accurate, because the strength reduction caused by ductility demand could be more significant than that by tension. • As the V/H ratio increases, more fluctuations can be observed in the columns axial force. This phenomenon leads to a significant reduction in the shear capacity in the range of (15-30) %. This reduction in shear capacity of the vertical members increases the potential for shear failure. • Geometric nonlinearities, in terms of the deformed configuration of the system, do not come into play in either IO or LS damaged states of the system. However, P-D effects due to higher interstory drifts of combined H ? V excitations do influence the response of the building in the region near collapse and must be taken into account. • A new vector IM is proposed in the current study to predict accurate fragility results. One of the main advantages of the proposed IM is its hazard compatibility, in which a GM prediction model can be easily developed for the second parameter (HM), implementing the existing attenuation models with an arbitrary set of periods. In case a probabilistic seismic hazard analysis (PSHA) would be required, the calculations can be performed using Ln (HM) as IM in the same way as any single spectral acceleration value. • The results presented in this study indicate that coupling horizontal and vertical ground excitations increases the ductility demand. Therefore, it is highly suggested that the conventional response spectrum (RS) be replaced with a multi-component RS in the next version of seismic design codes. Doing this will consider the effect of vertical ground excitations on the enhanced seismic demand and will provide a better understanding to structural designers.
Taking into account the above observations, the authors would like to recommend for the next version of seismic codes to make sure that the structures locate within 25 km Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.