Hazard‐consistent simulated earthquake ground motions for PBEE applications on stiff soil and rock sites

Ground motion record selection is a standard step in state‐of‐the‐art performance‐based earthquake engineering (PBEE) applications. It links the structural response to seismic hazard of the site of interest. In this process, suites of hazard‐consistent ground motion recordings of a wide range of intensity levels are selected (and often scaled) from a database of ground motions to be used as input to nonlinear dynamic response analysis. The ideal practice would be to select ground motions from a database of recorded accelerograms in such a way that they are consistent with the seismic hazard at as many intensity measure (IM) levels as possible and for all important causative scenarios. However, the available databases are not heavily populated, especially for the large‐magnitude short‐distance scenarios. Therefore, synthetic ground motions are the natural candidates to enlarge the database of recordings for scenarios that are important to the site hazard but lack, completely or partially, recorded accelerograms. The usage of such recordings, however, can be recommended only after they are proven to be “realistic” through a battery of seismological and statistical tests. In this study, we propose and implement a multi‐fold testing framework for such an evaluation. We first generate a ground motion database of simulated ground motions (SDB) whose values of the parameters of the causative earthquakes are, to the extent possible, the same of those that caused the ground motions in a mirror database of recorded ground motions (RDB). We then utilize the conditional spectrum (CS) approach computed for a hypothetical rock site in Perugia, Italy, to select hazard consistent records separately from the RDB and SDB databases. A comprehensive set of comparisons of distributions of several IMs, computed from the selected sets, and of distributions of engineering demand parameters (EDPs), that such record sets induce on simple SDoF structures suggests the adequacy of using the proposed simulated ground motions for structural response analyses of structures on stiff and rock sites.


NOVELTY
• Proposing a framework for calibration of the spectral content and time modulating functions of a 3D stochastic ground motion simulation technique, based on the iterative comparison of reference earthquake recordings and synthetic replicas.• Proposing a comprehensive assessment methodology to evaluate the use of simulated ground motions for PBEEcompliant seismic risk assessment of structures based on the comparison of recorded and simulated ground motions at the stages of: Ground Motion Databases, CS-selected sets of records, and response of structural systems.

INTRODUCTION
State-of-the-art building-specific risk assessment is now routinely conducted following the performance-based earthquake engineering, PBEE, 1 approach.In this context, the performance of a structure is assessed via its response to ground motions consistent with seismic hazard at the site of interest.The seismic hazard, which is estimated through Probabilistic Seismic Hazard Analysis, or PSHA, 2 is defined as the rate of exceeding (or equaling), over a defined period of time, different levels of a ground motion intensity measure (IM), or intensity measures (IMs), at the site of interest due to all considered earthquake scenarios.Structural response, on the other hand, is typically estimated by means of nonlinear dynamic analyses that use sets of ground motions selected in such a way that the distributions of the IM, or IMs, considered in the PSHA are consistent with those that may be experienced at the site as per the PSHA results.The response of the structure can be measured in different ways.Here we represent it through fragility curves describing the probability of exceeding certain limit states whose onset depends on reaching given values of engineering demand parameters (EDPs), such as maximum interstory drift.For the discussion that follows, we consider only one conditioning IM, for example the pseudo spectral acceleration Sa(T) at a given oscillator period T, but the concepts are valid also for a vector of conditioning IMs.The PBEE procedure, as cast in the seminal work, 1 implicitly assumes that the EDPs can be accurately predicted by the conditioning IM (called IM* here, for clarity), and that all the other IMs have, conditionally speaking, a much lower importance in affecting the EDPs that measure the performance of the structure.However, this approximation is marginally tenable for most types of structures, as structural response is indeed impacted by several IMs. 3 Moreover, different IMs may be good predictors for different EDPs in the same structure. 4,5Therefore, the choice of an appropriate conditioning IM directly affects the record selection procedure, the development of the ensuing fragility curves and, in the end, the accuracy of the seismic risk assessment.Under this reasoning, the selection of sets of records that are, to the extent possible, consistent with those that a structure may experience during its lifetime (denoted as hazard consistency here) is of paramount importance.
Recently, several methodologies have been developed for the selection of hazard-consistent records.The Conditional Spectrum method (CS) proposed by Jayaram et al., 6 for example, enforces hazard consistency by selecting records coherent with the distribution of spectral accelerations at all vibration periods conditional on the spectral acceleration at a period T* of interest, Sa(T*).The generalized conditional intensity measure (GCIM) approach, developed by Bradley 7 provides a generalization of the CS to include IMs other than spectral quantities.
In most cases, unfortunately, databases of recorded ground motions are not sufficiently populated to provide sets of records that match the target conditional distribution of all considered IMs in a satisfactory manner.This limitation has an impact on the record selection scheme.This is especially true for rock sites, for which ground motions are often scarce, particularly those representing severe shaking from events of high magnitudes close to the site.In the past, this limitation has been circumvented by using recordings on both rock and soil sites appropriately scaled to the amplitude levels desired by the application at hand.While this approach is practical, it might potentially bias the structural responses.This is because (i) time histories recorded at soil sites usually present distributions of IMs inherently different than those recorded at rock sites due to the presence of site-effects, [8][9][10] and (ii) weak but scaled ground motions have arguably different characteristics than strong ground motions naturally at high amplitude levels (e.g., 11 ).
To avoid or, better, limit any possible bias, analysts usually screen the database before record selection, for instance, by retaining records only from events of magnitude, distance consistent with the scenarios controlling the site hazard and with soil shear wave velocity like that of the specific site. 12Furthermore, limitations are also usually imposed to the scaling factors to avoid using ground motions deviating excessively from reality.While some level of screening of the initial database is, of course, highly recommended, it comes at a price.Strictly adhering to the site conditions or to the magnitude and source-to site distance events pertinent to the case at hand can cause a significant reduction of the number of records. 13Hence, strictly enforcing those conditions would generally lead to the undesirable consequence of achieving a poor fit to the target spectrum. 14Thus, the analyst often faces a trade-off between the number of adequate records made available for selection and the desired level of hazard consistency achievable from this limited data set.
Although the findings of recent work 15 have decreased the legitimate concern of mixing scaled soil and rock ground motions for structural response assessment (provided that they are selected to be consistent with the site hazard via the CS method) it is incontrovertible that an injection of simulated ground motions to the pool can be seen as an appealing alternative option.Ground motion simulation aims to produce accelerograms for earthquake scenarios that are either missing or are barely present in the databases of recorded ground motions.Before using them for structural response assessment, however, these simulated ground motions must be proven to be realistic when compared to those recorded from earthquakes having the same event (and site) causative parameters such as magnitude (M w ), hypocentral distance (R hyp ) and the average shear wave velocity in the upper 30 m (Vs 30 ).Furthermore, these ground motions should be validated by comparing pairs of reference and simulated ground motions in terms of distributions of (i) IMs, (ii) EDPs, or (iii) a combination of both.
0][21] Due to its wide-spread use by engineers, the most employed validation IM is the spectral acceleration at several oscillator periods of interest.Several other IMs, such as peak ground velocity (PGV), peak ground displacement (PGD), Arias intensity (AI) and significant duration (Ds 5-95 ) have also been tested in the past. 22,23Complementary, several studies have based their validation proposals on the comparison of EDPs obtained from structural response analysis.4][25][26][27][28][29][30][31] For this purpose, single degree of freedom systems (SDoFs) are often used as a proxy for complex structures, however some studies make use of multi degree of freedom systems (MDoF) instead. 32his study investigates the implications, for risk assessment of structures at rock sites, of relaxing the criteria for record selection by injecting into the available databases of recorded ground motion a pool of synthetic ground motions.The proposed assessment compares the distributions of IMs, EDPs, and fragility curves obtained from sets of real records and simulated ground motions selected with the CS approach.We first introduce the database of recorded ground motions used as reference.Next, we introduce the stochastic ground motion simulation methodology and the calibration technique used for the generation of the synthetic ground motions populating the simulated database.Finally, we present the case study, results obtained for a set of nonlinear SDoFs characterized by different periods of vibration and hysteretical behavior.

RECORDED GROUND MOTIONS DATABASE (RDB)
To assemble the RDB, we considered a subset of the Engineering Strong-Motion (ESM) database of ground motions, 33 which contains accelerograms recorded in Europe and Middle East from a total of 2179 earthquakes and 2080 stations.ESM is complemented by a flatfile that includes verified metadata and IMs of the manually processed waveforms.The RDB subset considered in this study contains the information from earthquakes with M w in the range of [4, 8), R hyp ≤ 100 km and Vs 30 ≥ 400 m/s.These constraints are deemed adequate for the case study and limit the comparisons to the cases where we expect the simulation technique to perform well.As discussed later, this is the case of ground motions from short-distance scenarios less impacted by surface waves, which are not explicitly modeled by the herein considered simulated method, and generated for sites characterized by stiff soil or rock where severe non-linearities of the soil response are not expected.These constraints result in an RDB of 6912 three-component recordings (see 33 for additional details).

SIMULATED GROUND MOTIONS DATABASE (SDB)
The database of simulated ground motions, called here the SDB, was made consistent by design with the characteristics of the records in the RDB described in the previous section.The consistency was enforced in terms of the joint probability distributions of the main causative parameters of the scenarios contained in the RDB, that is, M w , R hyp , Z hyp and of the characteristics of the local soil, Vs 30 (see glossary for the definition of such quantities).Here, we decided to explicitly consider Z hyp (even if it is already reflected through R hyp ) because of its importance for the propagation of the rupture within the kinematic source model considered in the ground motion simulation methodology discussed next.In the following subsections, we first provide general information about the simulation technique used to generate the synthetic ground motions.Next, we discuss the calibration of the simulation technique to reproduce different IMs from the scenarios contained in the RDB, and finally, we describe the process of assembling the SDB by means of Monte Carlo sampling of the causative parameters of the scenarios of interest.

Ground motion simulation methodology
The simulated ground motions presented in this study were generated using the stochastic ground motion simulation technique introduced by Alvarez, 34 herein referred to as SGMSM.This model is based on the methodology originally presented by Otarola et al. and Ruiz et al., 17,18 but includes a series of modifications that improve the physical representation of the simulated earthquake ground motions.For example, the technique adopted here uses the kinematic source model proposed in Gallovič, 35 which allows for a distribution of the asperities on the fault plane consistent with the properties of the crust, and a modeling of the rupture process superior to the radial rupture models typically coupled to the original representation of the source in earlier application of the SGMSM, for example, 17,18 because it considers a depth-dependent rupture propagation velocity (thus the importance of Z hyp ).Furthermore, the representation of the source is also updated to account for the evolution of the corner frequency during the rupture phenomenon, as per the proposal of Dang et al. 36 Finally, the SGMSM incorporates a modified version of the post processing procedure recommended by Wang et al., 37 to include the inter-frequency correlation structure of spectral accelerations empirically derived from recorded ground motions.
The SGMSM constructs the acceleration Fourier Amplitude Spectrum (FAS) of the ground motion as a function of M w , R hyp , and Vs 30 : Where U is the mean simulated ground motion spectrum, f is the frequency, S 0 , is the modulated white noise with normalized power spectral density and M 0 is the source seismic moment.The ground motion spectrum may be modeled considering a point source model or a finite fault source model.Small earthquakes, that is, those with negligible rupture dimensions compared to the propagation distance, are typically represented by a point source.Larger magnitude earthquakes, with propagation distances like (or smaller than) the dimensions of the rupture, are described through finite-fault models composed of arrays of sub-sources.Each of these sub-sources radiates seismic waves following the omega square model, as defined by Aki, 38 where the energy emitted from the source is directly determined by M 0 and the stress drop parameter (Δσ).The path travelled by the seismic waves is approximated considering ray theory while the effect of the attenuation is modeled through semi-empirical filters representing geometrical spreading and anelastic attenuation.The technique also considers the (de) amplification of the FAS due to phenomena independent from the path of the waves, such as crustal amplification and high-frequency decay phenomena present in recorded earthquake ground motion.Finally, site-specific effects may also be considered by means of convolution with transfer functions representing the response of the stratified profile of the soil.
The phase spectrum of the simulated ground motions is modeled as a random variable in the form of the band limited Gaussian white-noise with finite duration, S 0 .A modulating function is applied to the white noise with the intention to include the characteristic temporal non-stationary nature of ground motion time series.The Saragoni and Hart modulating function, ω(t, e, η, t η ), 39 shown in Equation ( 2), is often considered in stochastic simulation methods, for example 17,18,40,41 and was also considered for our study.

𝛚
( , , ,   ) = Where e and η are parameters defining the shape of the function, a = exp(1.0)/e,b = -e⋅ln(η)/(1.0+e⋅(ln(e)−1.0),c = b/e, t η = f Tgm⋅ T gm .T gm is the duration of the signal for the computation of the window function, and f Tgm is a factor modeling the elongation of the windows.Previous implementations of this modulating function have considered different combinations of parameters.These are often adjusted for matching reference signals, for example. 17,18

Calibration of the ground motion simulation methodology
We calibrated the SGMSM based on an iterative process with the goal of minimizing the differences between the distributions of IMs of interest, extracted from the ground motions in the RDB, and those extracted from synthetic ground motions "replicas".These so-called replicas were constructed by generating simulated ground motions with the same values of the parameters M w , R hyp , Z hyp , and Vs 30 as those of the records in the RDB.Acknowledging the fundamental differences introduced in ground motions by these parameters, we classified the recorded scenarios into 12 M w -R hyp -Vs 30 groups based on the following criteria: 4.0 ≤ M w < 5.0, 5.0 ≤ M w < 6.0, M w ≥ 6.0; R hyp ≤ 50 km, 50 < R hyp ≤ 100 km; 400 ≤ Vs 30 < 600, and Vs 30 ≥ 600 m/s.Note that we did not consider Z hyp amongst the causative parameters used for the classification of the scenarios.This because of its negligible effect in the IMs considered in the calibration procedure, when compared to the other parameters.

Generation of the synthetic earthquake scenarios
The simulation of the earthquake scenarios requires the description of the source geometry, slip distribution and hypocenter location.Complementary, the model requires the characterization of the crustal structure necessary for the simplified estimation of the path travelled by the seismic waves and, if site-specific effects are considered, also the transfer functions used in the convolution of the ground motion to include the effect of the superficial soil layers.Finally, the SGMSM requires the input parameters describing the FAS.Due to the scarcity of site-specific information, and for our implementation of the SGMSM, we recurred to generic descriptions of the source and regional velocity models.We used the regression model proposed by Wells et al., 42 to define the geometry of the seismic source.This model allows the estimation of the dimensions of the rupture, and its associated uncertainty, based on the M w of the earthquake scenario.We considered a uniform spatial density function to quantify the slip rates of a seismic source.This assumption results in the allocation of the scenario-consistent slip, based only on the properties of the crust at the depth of each sub-source, and thus not guided by any pre-disposed distribution based on inversions of the source distribution of previous events.For this exercise, all finite-source models were discretized into 100 sub-sources.The geometry of the finite-source was completed by defining the orientation of the ruptured plane.Due to the lack of specific information for each one of the scenarios in the RDB; dip, strike and rake of the seismic sources were sampled as random variables with uniform distributions.The ranges of variation of these variables were considered as 10 • −90 • for the dip, 0 • −360 • for the strike and 0 • −90 • for the rake.Finally, we randomly located the location of the hypocenter within the length and the deepest (down-dip) half of the ruptured plane.
We considered a crustal structure, or regional velocity model, following the models suggested by Boore and Joyner 43 for hard rock.This model provides a generic S-wave propagation velocity profile (β) for the crust as a function of depth.To complete our model, we approximated the P-wave velocity profile (α) considering a constant velocity ratio of β/α = 0.67, which corresponds to a Poisson ratio of 0.25.The density of the crust was estimated with the generic function for hard rock proposed in. 43Site effects were considered by sampling site-specific elastic transfer functions as a function of the reference Vs 30 .These transfer functions were computed with the Haskell-Thompson method, 44 considering the soil profiles of stations taken from the Kik-Net database, 45 and contemplating them to be directly on top of the velocity model defining the crustal structure.In other words, the computation of the transfer functions considers the uppermost layer of the crustal structure as the infinite half-space over which the soil column lies.
Finally, the input variables defining the ground motion spectrum were taken from Bindi and Kotha, 46 where the authors used the Generalized Inversion technique to parameterize the FAS of earthquake ground motions contained in the complete ESM database. 33

Calibration of the simulation technique
Currently, there are no available studies providing guidance on the modeling of the time modulating window functions within a finite-fault context.The lagged summation of sub-sources contributions does not allow the definition of a single set of parameters suiting all simulation possibilities.Additionally, the consideration of stress drop estimates obtained from the spectral inversion of large number of events results in estimates heavily biased toward the source physics of lower magnitude events.This is simply due to the scarcity of recorded earthquakes with medium or large magnitudes.In the preliminary phase of this work, these discrepancies resulted in simulated ground motions with unrealistic durations and overly large high-frequency content.We addressed this deficiency using an iterative calibration process to optimize our target modeling parameters including: (i) the time modulating window parameters, that is, e, η, and f Tgm , modeling the location of the peak, attenuation, and elongation of the window; and (ii) a correction factor, CF Δσ , applied to the stress drop.The goal of this optimization process was to minimize the differences between the IMs of interest extracted from the RDB ground motions and those extracted from sets of synthetic ground motion replicas computed with different combinations of the target modeling parameters.The replicas were constructed by generating one simulated ground motion per record in the RDB using the same causative parameters, that is, M w , R hyp , and Z hyp , of the recorded ground motion.The difference in the recorded and simulated ground motions in terms of the j th IM can be defined by Equation (3).
In Equation ( 3), IM ref,ij and IM sim,ij show the corresponding values of the IMj in i th pair of recorded and simulated ground motions.To find the total error of a set of simulated ground motions, one can compute the sum of the square of the error terms obtained by Equation (3) for multiple IMs (i.e., total of N IM ) and for all the records in the database, N rec , according to Equation (4).
In Equation ( 4), w j represents a weighting factor allowing the prioritization of the match of IM j , j = 1,. . ., N IM considered in this procedure.Here, we considered a vector of six target IMs including spectral accelerations, Sa, at periods of vibration of PGA, 0.1, 0.3, 0.5, and 1.0 s, and Arias Intensity, AI, with weights of 0.4, 0.1, 0.1, 0.1, 0.1, and 0.2, respectively.These IMs were selected because of their representativity of the optimal range of performance of the SGMSM, 17,18 that is, the highfrequency range.More specifically, the IMs were considered for the Orientation-Independent Geometric Mean, using period-independent rotation angle component GMRotI50, computed as per the procedure detailed in. 47We selected this specific component for consistency with the record selection scheme used in the coming sections.Note that herein we performed the optimization process separately for each one of the 12 M w -R hyp -Vs 30 groups that we defined earlier because of the inherent differences in the IMs of ground motions with different causative parameters.For each one of the groups, we generated a total of 25 sets of synthetic ground motion, each with different time modulating window parameters and stress drop correction factor selected by means of Latin Hypercube Sampling. 48For each of the groups, we defined the optimal set of parameters by selecting the set resulting in the lowest error and we consider this set to simulate the ground motions populating our SDB.The number of considered sets was defined based on acceptable computational times and errors.
Figure 1 shows the results of the calibration procedure for the case of R hyp ≤ 50 km.Figure 1A,B show the timemodulating functions for simulations obtained for soil sites (400 ≤ Vs 30 < 600 m/s), and rock sites (Vs 30 ≥ 600 m/s).These results show the importance of the site-effects on the duration and overall distribution of the amplitudes in the simulated ground motions, here the changes in the phase produced by softer soil layers disguise the differences in the arrival times of the energy of events with different magnitudes.For example, see how the modulating functions are fairly similar for all ranges of magnitude in Figure 1A, and how these change in Figure 1B, where site-effects are less relevant.
Figure 1C shows the correction factors for the stress drop, CF Δσ .These were noticed to consistently decrease with magnitude, except for the simulations at soil sites and for the largest magnitude bin.The decreasing correction factor as a function of magnitude is related to the influence of the size of the considered finite-source model.The point-source model concentrates the source at a single point, thus resulting in attenuation effects computed for a single R hyp .Finite-source models, on the other hand, are subject to different values of attenuation for every sub-source composing the model, and this results in more energy release from sub-sources closer to the observation point and, therefore, requiring a smaller correction factor.The case of the correction factor related to simulations for soil sites and for the highest magnitude bin is influenced by the non-linear response of the soil in the recorded events.This is a feature that is not captured in our simulated ground motions as site effects were computed considering linear elastic transfer functions.These functions are not capable of modeling the dissipation of energy produced by non-linear soil response.

Generation of the ground motion simulation database
After conducting the calibration of the SGMSM, for each one of the 12 predefined RDB groups, we populated the SDB with 7000 3D acceleration time histories generated by Monte Carlo sampling of the joint probability density function of the causative parameters of interest, that is, M w , R hyp , and Z hyp .The consideration of site-effects, source geometry and crustal structure model coincide with the methodology described for the calibration of the SGMSM.Figure 2A,B show the M w -R hyp scatter plot of the recorded and simulated ground motions of the RDB and SDB, respectively.Figure 2C compares the distribution of the Vs 30 in both databases.The minor differences observed in Figure 2C are mainly due to the limited number of sites contained in the considered database of transfer functions.As a final validation step in verifying the coherence of the ground motions in the synthetic versus recorded databases, Figure 3 compares IMs for two extreme M w -R hyp -Vs 30 groups.These groups are of special interest due to the number of ground motions considered in the record selection procedure discussed next.The comparisons of the response spectra, shown in the first row of Figure 3, indicate an overall similarity between the median spectra of the SDB and RDB ground motions, especially for the lowest magnitude group (5 ≤ M w < 6).Overall, we observed differences in the dispersion of the group with the largest magnitude range, M w ≥ 6, here represented by the shaded area that identifies the region between the 16th and 84th percentiles of the distribution of spectral accelerations.We attributed the lower dispersion in the RDB case to the lack of diversity in the earthquake scenarios populating this group, that is, most of the records were generated by the same earthquake, and thus share the same features related to the earthquake source.This represents a fundamental difference with our simulation scheme, where each simulated record comes, by definition, from a different earthquake scenario defined by the Monte Carlo sampling of all the input parameters.Similar distributions were observed for AI shown in the third row of Figure 3.The two IMs related to duration of the ground motion, that is, cumulative absolute velocity (CAV) and significant duration of the strong ground motion phase (Ds 5-95 ), however, have more distinct distributions.This is expected because these two quantities were not explicitly considered in our calibration process.In the next section we will also compare the distributions of the IMs of the subsets of ground motions extracted from these databases using the CS-based selection scheme.

ASSESSING THE ADEQUACY OF SYNTHETIC GROUND MOTIONS FOR PBEE APPLICATIONS
This section intends to investigate the suitability of the proposed simulated ground motions for utilization in the framework of the state-of-the-art PBEE seismic risk assessment approach.The proposed methodology considers three types of comparisons involving the hazard-consistent selected ground motions in terms of (a) parameters of the causative earthquakes of the selected accelerograms, (b) their IMs, and (c) EDPs they cause to SDoF representations of structures.These comparisons are made to identify the ranges of applicability of these synthetic ground motions without introducing any significant bias into the structural responses when compared with those obtained from recorded ground motions.

Case study structures, site hazard, and record selection
We studied the responses of an ensemble of nonlinear SDoFs with five fundamental periods of vibrations of 0.2, 0.5, 1.0, 1.5, and 2.0 s representing a range of different type of structures, with elastoplastic, elastic-hardening and degradation (pinching) hysteretic behaviors (shown in Figure 4).The structures are assumed to be located at a fictitious rock site in Perugia, Italy, with coordinates of 43.11 • N and 12.39 • E and Vs 30 = 800 m/s.The SDoFs were designed for a yield base shear coefficient, C y , computed as the lateral strength equivalent to that of the site-specific PSHA-based spectral acceleration value of 10% in 50 years at the fundamental period of each SDoF.C y is the yield base shear, V y, normalized by the weight W, and is numerically equivalent to the yield spectral acceleration Sa y .Sa y /g = V y /W is obtained by C y = Sa des (T 1 )⋅Ω/q where Sa des (T 1 ) is the design spectral acceleration at T 1 , q is the behavior factor assumed to be equal to 4.0 for new ductile buildings, and Ω = 2.0 is the over-strength factor.Accordingly, the corresponding yield displacement δ y , of the SDoF is obtained by δ y = Sa y ⋅[T 1 /(2⋅π)] 2 .Herein, we conventionally assume that collapse of a SDoF occurs once a ductility value of 8 is exceeded.As for the dynamic properties, the SDoFs were assumed to have a 5% mass proportional Rayleigh damping.
We performed the PSHA calculations using the OpenQuake 49 based on the area source model of SHARE 50 for spectral accelerations at T 1 = {0.2,0.5, 1.0, 1.5, 2.0}s using the ground motion prediction equation of Boore and Atkinson (2008). 51To evaluate the structural responses, we use the multiple stripe analysis (MSA) approach 52 using suites of hazard consistent ground motions separately extracted from the RDB and the SDB.To do so, we adopted the standard CS-based record selection approach of Jayaram et al., 6 with T 1 as conditioning IM* to select suites of records for the 10 IM levels (IMLs) corresponding to 70.0, 50.0, 30.0, 10.0, 5.0, 2.0, 1.5, 1.0, 0.6 and 0.2 % probability of exceedance (PoE) in 50 years.At each IML, we selected 40 pairs of records based on GMRotI50 of Sa that collectively match the target CS.Note that we limited the number of records generated by the same earthquake to a maximum of three to preserve, to the extent possible, the independence of ground motions utilized to match the target spectrum.A maximum scaling factor (SF) of 10 was allowed in the selection to give the priority to the matching of the record set with the target CS.Note that the selection is based on the "approximate" CS approach, 53 and that this method uses the mean of the causative scenarios to select hazard-consistent sets of records.The misfit of the sets of selected records is evaluated using the sum of squared errors SSEs 12 shown in Equation (5).
In Equation ( 5), m lnIMk is the sample mean of lnIM k and s lnIMk is the sample standard deviation of lnIM k values of the selected ground motions.The quantities µ lnIMk and σ lnIMk are the target conditional means and standard deviations, p is the number of periods of interest, and w is a weight factor here assumed equal to 2.0.This value assigns a higher degree of importance to the mismatches in the standard deviation rather than those in the target mean.In this study, SSE s = 0.12 is considered as the acceptable threshold based on the discussion presented in Garcia et al. 15 Figure 5 shows the overall fit of the selected suites of records by comparing the SSE s for all conditioning periods with the acceptable ratio.Overall, we observe that the quality of the match reduces as the IML increases, specifically for the RDB shown in Figure 5A, and for the highest periods of vibration.This indicates that the RDB lacks the diversity of spectral shapes required to match the CS at the highest IMLs.This observation is not as severe for the SDB, shown in Figure 5B, because we generated a different scenario for each M w -R hyp pair for the population of the database.In contrast, for the case of the RDB, many of the M w -R hyp pairs come from the same earthquake scenario, thus limiting the diversity of the spectra of the available records. Figure 6 shows the hazard consistency verification, following the procedure detailed in, 54 comparing the mean annual rates of exceedances from the PSHA-obtained hazard curves for Sa(0.2 s), Sa(1.0 s) and Sa(2.0 s), with the mean annual rates of exceeding the same values extracted from the CS selected sets of records with IM* = Sa(1.0s).Our comparison showed an acceptable match, and thus hazard consistency, for all sets of records and conditioning IM* considered in our study.

Comparison of earthquake and site parameters
As a first step in the assessment of the adequacy of our simulated ground motions for PBEE applications, we compared the distributions of the earthquake causative parameters and site parameter of the CS-based selected records from the RDB and SDB. Figure 7 shows for three of the five SDoFs, with T 1 = {0.2,1.0, 2.0}s, the comparison of the empirical distributions Comparison of the distributions of the causative parameters M w , Vs 30 , R hyp and of SF from first to the last row, respectively.Note: The boundaries of each box correspond to the lower and upper quartiles (16th and 84th, respectively), the line within the box corresponds to the median, and the whiskers extend to the minimum and maximum observed values.
of M w , Vs 30 , R hyp and, although not a causative or site parameter, also of the SF for the ground motions selected for all 10 IMLs.The boundaries of each box correspond to the lower and upper quartiles (16th and 84th, respectively), the line within the box corresponds to the median and the whiskers extend to the minimum and maximum observed values.Similar plots for the other two SDoFs can be found in the Supplementary material.Overall, there is an acceptable match of the distributions of the studied parameters for all conditioning periods.The similarity of these distributions confirms the adequate representation of the spectral form of the simulated ground motions.We did notice a larger dispersion in the distribution of the SF for T = 2.0 s, which indicates that spectral accelerations for this period of vibration are different in both sets of selected ground motions.Note that this specific vibration period was not considered in the calibration procedure of the simulation technique and, therefore, some discrepancies were expected.Nonetheless, we are presenting these results to show the robustness of the performance of the simulated ground motions in reproducing the IMs even for this limiting case.

Comparison of intensity measures (IMs)
Figure 8 compares the distributions of the IMs: PGA, AI, Ds 5-95 , CAV, and spectral intensity (SI), computed from the sets of records selected from the RDB and the SDB.Overall, there is a quite close fit for all the distributions of the considered Comparison of the distributions of PGA, Ds 5-95 , AI, CAV, and SI, from first to the last row, respectively.The boundaries of each box correspond to the lower and upper quartiles (16th and 84th, respectively), the line within the box corresponds to the median, and the whiskers extend to the minimum and maximum of the observed values.
IMs. Specifically, the IMs related to the spectral content of the ground motions, that is, PGA and SI, which were explicitly and implicitly considered in the calibration procedure, respectively, showed the best observed match.The match of the distributions is poorer for AI, Ds 5-95 and CAV that are related (directly or indirectly) to the ground motion duration, which was not accounted for in the calibration procedure.These differences were most noticeable in the median of Ds 5-95 , and in the overall ranges of variation of AI and CAV.We attributed the discrepancies in the distributions of IMs to different reasons.Firstly, the duration was not directly targeted in the calibration procedure.This is because we focused our attention on a CS-based record selection scheme.Duration, nevertheless, was indirectly accounted for in the calibration procedure by means of the AI.From another point of view, the number of events considered for the calibration of the group containing the largest magnitudes was the smallest within the reference ESM subset, thus limiting the available information for its characterization.Finally, the herein considered site effects were included by means of elastic transfer functions, therefore, differences in the amplitude of the simulated ground motions may be expected for seismic scenarios prone to result in the nonlinear response of the considered soil columns.Furthermore, the sampling of the site effects based only on Vs 30 may not guarantee similar changes in the phase of the ground motions.This is because different distributions of soil layers may result in comparable Vs 30 but completely different phase and spectral amplitude effects.In other words, Vs 30 is not a sufficient variable to describe site effects, as stated earlier in many articles available in the literature.

Comparison of structural response
This section presents the results obtained from MSA for multiple EDPs including ductility ratio (µ), peak acceleration and velocity of the system, and the dissipated energy.These quantities were computed from the response of the SDoF to the rotated geometrical mean of the horizontal components resulting in the GMRotI50 spectral acceleration considered for record selection.This procedure maintains the consistency between the CS-based selection scheme and the input to the nonlinear dynamic analyses.Figure 9 shows the comparison of the response of these EDPs, for the same three SDoFs presented earlier, and for the pinching material model.The EDP results for the remaining two SDoFs and hysteresis models, which show similar trends, are included in the Supplementary material.In Figure 9, the runs leading to ductility ratios above eight were treated as collapses, and, for illustration, they were lumped at the dashed horizontal grey line defining the maximum limit for each EDP and marked with an 'x' marker.When more than 50% of data points were collapse cases, the median ductility is infinite and, therefore, is omitted in the figure.The continuous lines describe the median response, the dashed lines refer to the 16th and 84th percentiles, and finally the individual points identify the results obtained for single analyses.For the case of dissipated energy, a value of 1.0 was fixed for records resulting in a linear response of the SDoF, that is, µ ≤ 1.0.In general, we observed that the distributions of responses obtained for the records selected from the RDB and the SDB are in good agreement for all considered SDoFs and IMLs.We noticed slight differences when the SDoF systems accumulate damage due to non-linear incursions, that is, for IMLs above the design level.The accumulated damage results in the elongation of the period of vibration of the SDoF, consequently shifting its response to one more influenced by lower frequencies in earthquake ground motion.It is well known that the SGMSM does not accurately represent this frequency range, for example, 17,18,34,40 mainly because of the importance source and wave propagation phenomena (oversimplified in stochastic techniques).The monitoring of the dissipated energy allowed us to investigate the effect of the differences in the distributions of the significant duration of the CS-based selected sets of records.Herein, statistically speaking, we do not see any considerable effects due to the differences in duration, which would be captured by differences in the dissipated energy.We argue that the differences in duration are not enough to cause systematically larger EDPs in the considered structures.
To complete our assessment of comparison of SDoF responses, we generated fragility curves, which provide the probabilities of exceeding a given limit state for different levels of IM, as well as the response hazard curves in terms of ductility factors.The limit states considered here refer to minor damage, moderate damage, and ultimate damage (collapse) states whose onset are defined here by the ductility ratios µ = {2, 4, 8}.The curves for the three SDoFs discussed so far with the pinching hysteresis model are compared in Figure 10.Again, the omitted cases are included in the Supplementary material.Overall, we observed a very good agreement between the sets of fragility curves derived from records selected from the RDB and SDB.In agreement with the trends noticed throughout our analysis, some slight differences were noticed for the largest IMLs, and for the longest vibration period SDoF.These differences cause an offset of the fragility curves for the collapse limit state for the SDoF systems with the longest periods of vibration.These differences, however, are not there for lower damage states, which suggests that the discrepancies appear due to the shifting of the response of the structure due to the period elongation phenomenon discussed earlier.In addition, we noticed slight differences in the dispersion of the moderate and collapse fragility curves for the SDoF with T 1 = 0.2 s.We attributed these discrepancies to the overall larger dispersion associated to the simulated records selected for the highest IMLs and for this SDoF, that is, 6 < M w and R hyp < 50 km.These differences may be seen, for example, in the comparison of the response spectra shown in Figure 3, where we observed that our simulations resulted in spectral accelerations characterized by a larger dispersion when compared to those from the RDB.As mentioned earlier, we attributed this divergence in dispersion, at least partially, to the independence of the M w -R hyp pairs considered for the simulation of the earthquake scenarios F I G U R E 9 Comparison of the distributions of the Ductility ratio (µ), peak acceleration, peak velocity, and dissipated energy, from first to the last row, respectively, computed for three SDoFs with pinching hysteresis model.The values of these EDPs were obtained via nonlinear dynamic analysis using the rotated geometrical mean of the horizontal components resulting in the GMRotI50 of 40 earthquake scenarios at each IML selected from the RDB and SDB.Median (solid line) and 16th and 84th percentiles (dashed lines).The individual analysis marked with an 'x' refer to the records leading to collapse.populating the SDB.In fact, this is a fundamental difference in the construction of our SDB with respect to the RDB and its impact becomes more evident when comparing groups with fewer scenarios, where records often share the same seismic source.

CONCLUSIONS AND RECOMMENDATIONS
This paper investigated the adequacy of using synthetic ground motions simulated based on the methodology introduced by Otarola et al. and Ruiz et al., 17,18 and modified by Alvarez 34 in a PBEE framework for applications on stiff and rock sites.For this purpose, the appeal of using simulated ground motions to supplement missing or scarce recorded ground motions especially for those earthquake scenarios controlling the site hazard is unquestionable.The responses, crafted in this study in terms of both fragility curves and response hazard curves for three different damage states ranging from minor to collapse, were computed for an ensemble of five structures modeled as nonlinear SDoF systems with fundamental period of vibrations of 0.2, 0.5, 1.0, 1.5, and 2.0 s with three alternative hysteretic models, elastoplastic, hardening and pinching.These structures were assumed to be located at a rock site in Central Italy for which we evaluated the seismic hazard via PSHA.The fragility curves for these structures were analytically computed via nonlinear dynamic analysis using a MSA approach fed by either simulated or recorded ground motions selected in both cases to match site-specific conditional spectra at different intensity levels.The response hazard curves were obtained by convolving the fragility curves for the three damage states with the corresponding hazard curves for the site.
To allow for a fair comparison of these two sets of fragility curves, we assembled a database of simulated ground motions (SDB) that were generated from earthquake scenarios with causative parameters (M w , R hyp , and Z hyp ) and soil conditions (Vs 30 ) consistent, to the extent possible, with those of the recorded ground motions included in the reference database (RDB).This action ensures that the discrepancy in the distributions of IMs that matter for the prediction of the EDPs is kept to a minimum (or at least not produced by differences in the distribution of the causative parameters of the earthquake scenarios populating the databases).
The ground motions, either recorded or simulated, selected by the CS method for different intensity levels were able to satisfactorily match the target spectra.This finding ensures the capacity of the proposed simulation technique to adequately represent the spectral content of the motions for all earthquake scenarios populating the reference database.However, when comparing the distributions of other IMs related to ground motion duration, such as Ds 5-95 , and CAV, which were not included in the calibration process, we observed some differences in the median and standard deviation.These differences point to a limitation of the adopted simulation technique to reproduce duration-related IMs based solely on the matching of spectral content and energy.We intend to overcome this limitation in the future by enforcing a broader calibration of the simulation technique to include these IMs.
Moving to responses of these SDoF systems, we found no significant statistical discrepancies, for most of the studied cases, between the distribution of several considered EDPs, namely ductility, peak acceleration and velocity, and dissipated energy, generated by the recorded and simulated ground motions selected to be hazard-consistent via the CS approach.Consequently, the fragility curves for the three ductility-based damage states considered herein showed little differences for all SDoF systems and damage states, except for the case of collapse of SDoF systems with long elastic vibration periods.We attributed the observed discrepancies to different factors including the difficulty to match the target spectrum, and the limitations of the SGMSM to model the low-frequency content of earthquake ground motion (exacerbated by the elongation of the period of vibration due to non-linear response).The differences observed in the fragility curves computed for collapse are mostly noticeable when considering the pinching material model, where the systematic differences in duration of the records populating the databases play a role in the degradation of the response.These results suggest that if carefully calibrated to match the distributions of relevant IMs in tandem with careful rock hazard consistent record selection enforced via a CS method, there is no detectable statistical difference between the distributions of responses of structures on stiff and rock sites to real ground motions and simulated ones with this technique, at least for the set of structural systems considered herein.
To conclude, we do not advocate in real life PBEE applications to replace recorded ground motions with synthetic ones simulated for specific earthquake scenarios and soil conditions (here stiff soil and rock).However, the very good agreement of the response distributions showed here suggests that adding carefully calibrated and validated simulated ground motions to the mix of available accelerograms may be a good practice ready for prime-time consideration.The scarcity of recorded ground motions may hinder the match of the CS target, potentially leading to biased fragility curves and, in turn, of response hazard curves.This deficiency can be mitigated, if not eliminated, by adding simulated motions into the pool available for properly ensuring hazard consistency.F y yield base shear g acceleration of the gravity, (9.8 m/s 2 ) GCIM generalized conditional intensity measure GMRotI50 Orientation-Independent Geometric Mean, using period-independent rotation angle component IM(s) intensity measure(s) IML intensity measure level M 0 Seismic moment MDoF multi-degree of freedom system MSA multi-stripe analysis M w moment magnitude N rec number of records PBEE performance-based earthquake engineering PGA peak ground acceleration computed for component GMRotI50 PGD peak ground displacement computed for component GMRotI50 PGV peak ground velocity computed for component GMRotI50 PoE probability of exceedance PSHA probabilistic seismic hazard analysis q behavior factor RDB recorded database R hyp source-to-site hypocentral distance S 0 band-limited, finite duration and power spectral normalized gaussian white noised Sa(T) spectral acceleration, for component GMRotI50 and for period of vibration T Sa des (T 1 ) spectral acceleration, computed for component GMRotI50 at the fundamental period of vibration T 1 SCEC Southern California Earthquake Center SDB simulated database SDoF single degree of freedom system SF earthquake spectrum scaling factor SGMSM stochastic ground motion simulation method SI spectral intensity (), which is computed as  = We would also like to acknowledge Dr. Irmela Zentner for the constructive discussions that improved this article and Prof. Fabian Bonilla for his key advising role in developing the ground motion simulation technique adopted here.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The simulated ground motion data base and its correspondent flat file, produced for the purpose of this study, is available upon request to the first author.

F I G U R E 1
Optimal modeling parameters obtained based on the calibration procedure.(A) Time modulating functions for R hyp ≤ 50 km and 400 ≤ Vs 30 < 600 m/s; (B) time modulating functions for R hyp ≤ 50 and Vs 30 ≥ 600; and (C) stress drop correction factors obtained for R hyp ≤ 50 km.F I G U R E 2 Distribution of M w -R hyp in (A) RDB and (B) SDB; and (C) comparison of the distribution of Vs 30 for the ground motions in the RDB and SDB.

F I G U R E 3
Comparison of the IMs distributions for the ground motion records populating different groups considered in the calibration: Acceleration response spectra (with shaded area representing the region between the 16th and 84th percentiles of the distribution), and empirical probabilistic mass functions for (EPMF) for Significant duration (Ds 5-95 ), Arias intensity, AI, and cumulative absolute velocity and CAV from first to the last row, respectively.

F I G U R E 4
Hysteresis behaviors considered in this study, namely (A) elastoplastic, (B) elastic with hardening (hardening ratio of 1.03), and (C) pinching.Note: panel (C) includes the pinching material modeled through the IMKPeakOriented model in Opensees.A detail characterization of the degradation parameters as implemented in Opensees is provided in the Supplementary material.

F I G U R E 5
SSEs of the records selected from the (A) RDB and (B) SDB to match the CS at 10 IMLs for the five conditioning periods specified in the legend.F I G U R E 6 Hazard consistency verification for spectral ordinates of Sa(0.2 s), Sa(1.0), and Sa(2.0) using the IM* = Sa(1.0)as conditioning period for CS-based selection of records from: (A) RDB and (B) SDB.The dashed lines show the PSHA-based hazard curves while the solid lines show the values of the mean annual frequency of exceedance of each spectral ordinate as estimated from CS-based selected records.

F I G U R E 1 0
Comparison between the fragility curves (first row) and response hazard curves (second row) obtained from RDB (solid line) and SDB (dashed line) groups for three different ductility levels and SDoFs with vibration period of T 1 = 0.2 s, T 1 = 1.0 s, and T 1 = 2.0 s.Note: These results pertain to SDoFs modeled with the pinching hysteresis model.The vertical dashed lines in the second row show the ductility values at the onset of the three damage states.

AI
Arias intensity computed for component GMRotI50 CAV cumulative absolute velocity, computed as  =   ∫ 0 |()| , where   is the total duration of the ground motion a(t).Computed for component GMRotI50 CF Δσ stress drop correction factor CS conditional spectrum Cy yield base shear coefficient of a SDoF Ds 5-95 duration of the strong ground motion phase, computed as the duration between 5%−95% of the total Arias intensity.Computed for component GMRotI50 e Saragoni-Hart window parameter locating the peak of the function EDP(s) engineering demand parameter(s) EPMF empirical probabilistic mass function ESM European strong motion database FAS Fourier Amplitude Spectrum f Tgm Saragoni-Hart window parameter modeling the elongation of the function ), where  is the oscillator period and PSA is the pseudo-spectral acceleration.Computed for component GMRotI50 SSEs sum of squared errors T oscillation period T* conditioning period U simulated ground motion spectrum V s30 average shear wave velocity of the uppermost 30 m W structure's weight Z hyp depth of the hypocenter α P-wave velocity profile β S-wave propagation velocity δ y yield displacement Δσ stress drop ε logarithmic difference between recorded and simulated IMs η Saragoni-Hart window parameter modeling the decay of the function Ω over strength factor µ Ductility ratio ω Saragoni-Hart window function A C K N O W L E D G M E N T S This study was funded under ANRT CIFRE thesis contract 8610-5920074779 and by the Horizon 2020 program under grant agreement n • 945121.