Modification of stochastic ground motion models for matching target intensity measures

Stochastic ground motion models produce synthetic time‐histories by modulating a white noise sequence through functions that address spectral and temporal properties of the excitation. The resultant ground motions can be then used in simulation‐based seismic risk assessment applications. This is established by relating the parameters of the aforementioned functions to earthquake and site characteristics through predictive relationships. An important concern related to the use of these models is the fact that through current approaches in selecting these predictive relationships, compatibility to the seismic hazard is not guaranteed. This work offers a computationally efficient framework for the modification of stochastic ground motion models to match target intensity measures (IMs) for a specific site and structure of interest. This is set as an optimization problem with a dual objective. The first objective minimizes the discrepancy between the target IMs and the predictions established through the stochastic ground motion model for a chosen earthquake scenario. The second objective constraints the deviation from the model characteristics suggested by existing predictive relationships, guaranteeing that the resultant ground motions not only match the target IMs but are also compatible with regional trends. A framework leveraging kriging surrogate modeling is formulated for performing the resultant multi‐objective optimization, and different computational aspects related to this optimization are discussed in detail. The illustrative implementation shows that the proposed framework can provide ground motions with high compatibility to target IMs with small only deviation from existing predictive relationships and discusses approaches for selecting a final compromise between these two competing objectives.


Summary
Stochastic ground motion models produce synthetic time-histories by modulating a white noise sequence through functions that address spectral and temporal properties of the excitation. The resultant ground motions can be then used in simulation-based seismic risk assessment applications. This is established by relating the parameters of the aforementioned functions to earthquake and site characteristics through predictive relationships. An important concern related to the use of these models is the fact that through current approaches in selecting these predictive relationships, compatibility to the seismic hazard is not guaranteed. This work offers a computationally efficient framework for the modification of stochastic ground motion models to match target intensity measures (IMs) for a specific site and structure of interest. This is set as an optimization problem with a dual objective. The first objective minimizes the discrepancy between the target IMs and the predictions established through the stochastic ground motion model for a chosen earthquake scenario. The second objective constraints the deviation from the model characteristics suggested by existing predictive relationships, guaranteeing that the resultant ground motions not only match the target IMs but are also compatible with regional trends. A framework leveraging kriging surrogate modeling is formulated for performing the resultant multi-objective optimization, and different computational aspects related to this optimization are discussed in detail. The illustrative implementation shows that the proposed framework can provide ground motions with high compatibility to target IMs with small only deviation from existing predictive relationships and discusses approaches for selecting a final compromise between these two competing objectives. the most popular methodology for performing this task for seismic risk assessment (or seismic design) applications is the selection and modification of real (ie, recorded from past events) ground motions based on a target intensity measure (IM) level, 5,6 eg, an elastic pseudo-acceleration response spectrum. For seismic risk assessment, such modification is performed for specific seismicity scenarios (typically defined through moment magnitude and source-to-site distance) contributing to the seismic hazard for the chosen site, with the target IM commonly 7 derived through ground motion predictive equations (GMPEs). 8,9 An alternative philosophy for describing seismic excitations is to use simulated ground motions. 10,11 A specific modeling approach for the latter, which has been steadily gaining increasing attention by the structural engineering community, [12][13][14] is the use of stochastic ground motion models. [15][16][17][18][19][20] These models are based on modulation of a stochastic sequence, through functions (filters) that address spectral and temporal characteristics of the excitation. The parameters of these filters are related to seismicity (eg, moment magnitude and rupture distance) and site characteristics (eg, shear wave velocity for soil profile) through predictive relationships. 15,19 Sample ground motions for any desired seismicity scenario can be generated by determining the parameters of the stochastic ground motion model through these predictive relationships and by using a sample stochastic sequence. This approach may ultimately support a comprehensive description of the seismic hazard. 21 The essential component of stochastic ground motion models is the development of the associated predictive relationships, and various approaches have been established to accomplish this, with main representatives being record-based and physicsbased models. Record-based models (also known as site-based) are developed by fitting a preselected "waveform" to a suite of recorded regional ground motions. 15,18,22 On the other hand, stochastic physics-based models rely on physical modeling of the rupture and wave propagation mechanisms. 19,20 Emphasis in this study will be on the former models, although the techniques discussed can be extended to any type of stochastic ground motion model.
An important concern related to the use of stochastic ground motion models for structural engineering applications is the fact that through current approaches in selecting their predictive relationships compatibility to the seismic hazard for specific structures and sites is not necessarily obtained 15 (this is also shown in the illustrative example considered later). Although validation of these models is frequently performed by comparison of their spectral acceleration outputs to GMPEs, 15,17 the match to GMPEs is not explicitly incorporated in the predictive relationships development. Such a match to some desired GMPE (or target IMs in general) is though important for subsequent use of the stochastic ground motion models to describe the seismic hazard. Take for example the recent FEMA P-58 23 guidelines for seismic performance assessment of structures; the scenario-based description of the seismic hazard requires match of the median response to the one described by a GMPE (for the specific seismicity scenario examined). Similarly, the intensity-and time-based descriptions in FEMA P-58 require compatibility with seismic hazard curves, which are ultimately defined through use of GMPEs. 24 This realization has motivated researchers to investigate the selection of predictive relationships for stochastic ground motion models so that compatibility with GMPEs is explicitly established. 25 The formulation introduces an explicit optimization for matching the median predictions of the ground motion model to the spectral acceleration estimates of GMPEs, while maintaining physics-based principles or the matching to trends from real ground motions as an optimization constraint, in an attempt to preserve desired ground motion characteristics. Vetter et al 26 recently extended the work of 25 by providing a versatile and computationally efficient approach, leveraging surrogate modeling principles, for tuning stochastic ground motion models to establish compatibility with the median GMPE predictions for a range of structural periods and seismicity scenarios of interest. One of the main drawbacks of this tuning approach, though, is that the physical characteristics of the resulting acceleration time-series are incorporated in the optimization merely as constraints, something that requires significant experience in ground motion characterization for proper definition of the optimization problem and can furthermore lead to synthetic time histories with unrealistic properties for some seismicity scenarios.
The current study addresses this critical shortcoming and looks at the modification of stochastic ground motion models for specific seismicity scenarios with a dual goal of (i) matching a target IM for a specific structure (or range of structures) while (ii) preserving desired trends and correlations in the physical characteristics of the resultant ground acceleration time series. This is ultimately formulated as a multi-objective optimization problem. The first objective is to minimize the discrepancy between the median ground motion output and the target IM for a given seismicity scenario. Any desired IM can be used for this purpose with only requirement to have a corresponding seismicity scenario. For instance, if the target IM is derived through probabilistic seismic hazard analysis, a corresponding seismicity scenario (or 'design earthquake') can be derived through the disaggregation of seismic hazard 7 for a given hazard level. The second objective is to establish the smallest deviation from the model characteristics suggested by existing predictive relationships. This second objective aims at maintaining regional physical characteristics and parameter correlations with respect to existing predictive relationships. Approach differs significantly from Vetter et al 26 ; rather than tuning the ground motion for hazard compatibility ignoring any existing predictive relationships, goal here is the minimum modification of the existing relationships that will yield the desired compatibility. This is ultimately posed as a multi-objective problem to better investigate the compromise between the two different objectives, and for efficiently solving it a surrogate modeling approach is adopted, similar to that of Vetter et al. 26 A surrogate model (ie, metamodel) is trained based on an initial database of ground motion simulations, and ultimately provides a highly efficient approximation for the spectral acceleration predictions of the stochastic ground motion model. The surrogate model is then leveraged to solve the optimization problem. Emphasis is also placed here on the selection of the database to inform the metamodel development, which constitutes a significant advancement over the approach by Vetter et al. 26 Blind search and gradient-based approaches are considered for the multi-objective optimization, and the relative computational benefits of each are explored.
In the next section, the general problem of developing simulated ground motions compatible with target IMs is defined, and then specific aspects of the framework are discussed in detail.

| PROBLEM FORMULATION
Consider a stochastic ground motion model that provides acceleration time-histories € a tjθ; W ð Þby modulating a Gaussian white noise sequence, W, through appropriate time/frequency functions that are parameterized through the n θ -dimensional model parameter vector θ ¼ θ 1 θ 2 … θ n θ ½ ∈ ℜ n θ . This vector completely defines the model and is typically composed of various excitation properties such as Arias intensity, strong ground motion duration or parameters related to frequency characteristics of the ground motion. A specific example for such a model, the one used in the illustrative example later, is provided in Appendix A. This particular record-based model efficiently addressed both temporal and spectral nonstationarities. The former is established through a time-domain modulating envelope function, whereas the latter is achieved by filtering a white noise process by a filter with characteristics that vary in time.
Synthetic time-histories can be created by relating θ to seismicity and local site properties through predictive relationships. The vector of these properties, called seismological parameters, is denoted as z. Common characteristics used for z 15,19 include the fault type F, the moment magnitude, M, the rupture distance, R, and the shear wave velocity in the upper 30 meters of soil, V s30 . For record-based models, the standard approach for development of these predictive relationships 15,27 relies on first matching the waveform characteristics to recorded ground motions (ie, identify first θ for each of the recorded ground motions in a given database) and then performing a regression to relate θ to z. This leads ultimately to a probabilistic regression model for θ with mean predictive relationship μ θ (z) that is dependent on z and some associated uncertainty characterization U, identified from the residuals of the regression, that is independent of z. Typically, this is performed by first transforming the problem to the standard Gaussian space through a nonlinear mapping for each component θ i . The transformed Gaussian vector is denoted v(θ) herein. Approach ultimately leads to a Gaussian probability model v~N(μ(z),Σ) with mean μ(z) and covariance matrix Σ. In this case, the uncertainty characterization U corresponds to the covariance matrix Σ and to the fact that probability model for v is identified as Gaussian. Appendix A includes more details for a specific ground motion model. 15 Note that a similar description can be established for physics-based models. In this case, the predictive relationships μ θ (z) are obtained through rupture and wave propagation principles, 19,28 whereas the uncertainty characterization U can be established by assigning probability models for θ through an epistemic uncertainty treatment. 12,29 As discussed in Section 1, this formulation for the predictive relationships of stochastic ground motion models, prioritizing a match to regional trends, provides synthetic ground motions whose output IMs do not necessarily match hazardcompatible IMs (eg, as derived from GMPEs). For this purpose, a modification of the model parameter vector θ is proposed for specific seismicity scenarios defined by z with objective to (i) match a target IM vector, while (ii) maintaining similarity to the predictive relationships already established for the model. Equivalently, this can be viewed as identifying the model characteristics θ that are closest to the established model μ θ (z) (considering, when available, any additional information provided by U) and also match the intended hazard (described through some IM). The IM vector may include different response quantities of interest, for example, (i) direct characteristics of the ground motion, such as peak ground acceleration, velocity, and displacement; (ii) elastic and inelastic spectral responses for different periods of a single degree of freedom (SDoF) oscillator; or (iii) more complex spectral or ground motion related quantities proposed by different researchers. [30][31][32] The target for most of these IMs can be described through a GMPE. 8,33 However, this is not necessary; any IM description can be used, with only requirement that a corresponding seismicity scenario is defined. Note that if match to spectral responses (ie, a spectral plot) is of interest, then a range of structural periods for which the match is established needs to be determined.
To formalize these concepts mathematically, let, Y i (z) ; i = 1,.…n y denote the target response quantities of interest and Y m i θ ð Þ the median predictions for the same quantities provided through the stochastic ground motion model. The median predictions are obtained through the following process: Step 1. Generate n w sample acceleration time-histories for different white noise sequences € a κ tjθ; W κ ð Þ; κ ¼ 1; …; n w . Step 2. For each sample evaluate the responses of interest. For spectral quantities this will entail numerical simulation of SDoF responses.
Step 3. Estimate the statistics (median) over the established sample-set.
The modification problem is ultimately formulated as bi-objective optimization problem: The first objective F 1 corresponds to a measure of the discrepancy from the chosen target. One choice for this measure is the average weighted square error given by with γ i corresponding to the weights prioritizing the match to different IM components. A typical selection for γ i is 1/Y i (z), then the quantity in Equation 2 corresponds to the average squared relative error. Alternative formulations for this first objective, facilitating perhaps a better physical intuition, are the average absolute error or the maximum absolute error over the different IMs, given, respectively by, Both these objectives lead, though, to more challenging optimization problems as they correspond to nondifferentiable functions. Thus, preference will be here for the objective given by Equation 2. The alternative measures will be used to evaluate the suitability of different solutions. The second objective F 2 measures the discrepancy of θ from the established predictive relationships. One choice for F 2 (θ|z) could be [θ-μ θ (z)] Τ [θ-μ θ (z)], ie, the discrepancy from the mean predictive relationships. If there is no available uncertainty characterization U, then this is the only option that can be made. If this characterization is available, as is the case typically with record-based models, 15 correlation and variability information for the model parameters can be additionally incorporated. This is established by considering the maximization of the likelihood of θ based on the probability model described through U. For the stochastic ground motion model described in Appendix A, this corresponds to the maximization of the probability density for v(θ), leading to The covariance matrix Σ incorporates in the formulation the correlation between the model characteristics as well as the fact that variability for each of these characteristics is different. Objective function F 1 enforces the match to the target IMs. Objective F 2 guarantees compatibility of the physical characteristics of the resultant ground motions with the regional trends observed in recorded ground motions. Solution of the multiobjective optimization of Equation 1 ultimately leads to a Pareto set of dominant solutions {θ p ; p = 1,…,n p } that express a different compromise between the competing objectives F 1 and F 2 . A solution is characterized as dominant (and belongs in the Pareto set) if there is no other solution that simultaneously improves both objectives F 1 and F 2 . The representation of the Pareto set in the performance objective [F 1 , F 2 ] space, {[F 1 (θ p |z), F 2 (θ p |z)]; p = 1,…,n p } is termed as the Pareto front. Illustrations of such Pareto fronts are included in the example discussed later. One extreme point of this front will always correspond to the unmodified model with θ = μ θ (z), representing the minimum of objective F 2 = 0. Unless this point also yields a match to the targeted hazard (ie, corresponds to F 1 = 0), optimization of Equation 1 will identify points that improve upon F 1 (μ θ (z)| z) while deviating from the unmodified model (F 2 > 0). One can eventually select a model configuration from the identified Pareto set that yields the desired IM-compatibility without deviating significantly from regional ground motion characteristics. This will be further discussed in the illustrative implementation in Section 5.4.
Identifying the Pareto set for this problem is challenging because the computational burden in evaluation of objective F 1 is significant, requiring n w ⋅ n y time-history analyses for each objective function estimation. To facilitate an efficient optimization that can be repeated for any desired seismicity scenario z, a surrogate modeling approach is adopted here, similar to the one used by Vetter et al. 26 Specifically, kriging is used as metamodel since it has a proven capability to approximate highly complex functions, 34 while simultaneously providing gradient information that will be directly exploited in the optimization and allowing to explicitly consider the local metamodel approximation error within the optimization formulation. These aspects of the optimization problem will be discussed in Section 4. The details of the kriging metamodel development are discussed first in the next section.

| KRIGING METAMODEL DEVELOPMENT
The kriging metamodel is developed to provide an efficient approximation to the input-output relationship θ −Y m i θ ð Þ considering every potential response quantity of interest that can be eventually used for the definition of objective F 1 . A further simplification can be established if the relationship between some components of θ, and the response output Y m i θ ð Þ is explicitly known. This is true for stochastic ground motion models that include a parameter, denoted θ s herein, that directly controls the amplitude of the excitation. This means that Y m i θ ð Þ ¼ θ s s m i x ð Þ with x corresponding to the remaining model parameters Without loss of generality, we will adopt here this assumption, ie, representation In this case, the metamodel needs to be established to approximate only relationship x−s m i x ð Þ. For developing the metamodel, a database with n observations is initially obtained that provides information for the x−s m i x ð Þ pair. For this purpose, n samples for {x j , j = 1, … , n}, also known as support points or experiments, are obtained over the domain of interest for x. This domain, denoted X, should encompass the anticipated range that the metamodel will be implemented in and its determination is discussed later in this section. The median predictions provided through the ground motion model are then established through the 3-step process discussed in Section 2 considering n w white noise samples. Using the data set É , the kriging model can be formulated. Details for this formulation may be found in the studies of Sacks et al 35 or Vetter et al, 26 with the latter reference focusing on a similar application as the one considered here, looking at approximating the predictions of stochastic ground motion models.
This approach ultimately leads to a kriging predictor that has a Gaussian nature with predictive mean b s m i x ð Þ and local prediction variance, which is also a function of x, σ 2 i x ð Þ . 35 Each response output can be approximated through this predictor leading to where ε i is a standard Gaussian variable. This facilitates a computationally efficient approximation to Y m i θ ð Þ for each θ. This efficiency can be improved by ignoring the metamodel prediction error (ie, setting ε i = 0) since calculation of predictive variance σ 2 i x ð Þ entails a significant higher computational cost than estimation of predictive mean. The computationally intensive aspect of the entire formulation is the development of the database , which requires response-history analysis for a large number of model parameters to populate X and a sufficient number of white noise samples to address the resultant variability in the response. However, this needs to be performed only once. As soon as the kriging metamodel is established based on this database, it can be then used to efficiently predict the responses for any other θ desired. Calculation of b s m i x ð Þ and σ 2 i x ð Þ can be also vectorized, 36 something that will be leveraged in the numerical optimization discussed in the next section.
The accuracy of this metamodel depends on the number of experiments n used as well as the exact selection of these experiments. A larger value of n can improve accuracy but at the same time can reduce significantly computational efficiency. Because the latter is important for solving the challenging multi-objective optimization problem discussed here, the value of n needs to be kept moderately low. Therefore, the improvement of metamodel accuracy is primarily sought after through an adaptive design of experiments (DoE). The adaptive DoE strategy gradually increases the number of support points, leveraging the metamodel developed through the existing support points to guide the selection of the new experiments. This leads to an iterative identification of support points, whereas the specific strategy adopted here for selecting the support points in each of the iterations corresponds to a sample-based DoE. 37,38 In the first iteration, since no metamodel is yet available, the initial n 1 experiments are obtained using Latin hypercube sampling in X, ie, a space-filling DoE. Subsequent iterations adopt the adaptive DoE strategy. At the kth iteration, a surrogate model is developed using the available n av (k) support points. The prediction error of this metamodel is then leveraged to identify new experiments in regions with low metamodel accuracy. This is accomplished through a sample-based implementation: a large number of candidate experiments is first sampled within X that are distributed proportional to σ 2 i x ð Þ (for example, through rejection sampling 39 ), and these experiments are then clustered (for example, using k-means clustering 40 ) to the desired number n a of additional experiments, which are added to the existing experiments for a total of n av (k+1) = n av (k) + n a support points. The clustering is important for avoiding close-proximity support points that ultimately provide information over the same domain in X. A new metamodel is then developed using the n av (k+1) experiments, and its accuracy is assessed, for example, by calculating error statistics through cross-validation. 38 If sufficient accuracy is achieved, then adaptive DoE is terminated and n = n av (k+1) , else algorithm proceeds to the k + 1st iteration.
Another important feature of the metamodel development is the definition of domain X. This domain needs to (i) cover the entire range of values that the metamodel will be eventually used for (to avoid extrapolations that unavoidably have reduced accuracy) while (ii) avoiding unnecessarily broad definitions that lead to computational effort spent in subdomains within X of no practical interest. Here the definition of X is established by ignoring subdomains that correspond to high values for F 2 (θ|z) for any potential seismicity scenario z, ie, x values that are away from the current predictive relationships. Solutions within such subdomains will not be selected in the multi-objective optimization because they correspond to large values for one of the objectives. The X definition is established through these steps: Step 1. Create a range of seismicity scenarios z l c ; l ¼ 1; …; n c È É that are representative of the scenarios that could be eventually considered in the ground motion model tuning.
Step 2. Define a box-bounded domain X d that is expected to be a superset of X, and create a large number of samples for for each seismicity scenario z l c where θ j d is the sample corresponding to x j d and the mean value for θ s (for the given z l c ). Evaluate which corresponds to the smallest distance for x j d from the current predictive relationships for any potential seismicity scenario.
Step 4. Set a threshold δ d and then classify each sample x j d as belonging in set X if D j is smaller than δ d (sample is given membership classification 1) or not otherwise (sample is given membership classification 0). Using this classification information, domain X can be characterized through support vector machine (SVM). 41 This approach ultimately leads to an SVM characterization of X. Figure 1 demonstrates some of the steps of this process. Finally, samples within X, as needed for the adaptive DoE, can be generated with negligible computational effort by creating first samples within box-bounded domain X d and then using the SVM classifier to maintain only the samples belonging in X.

| MULTI-OBJECTIVE OPTIMIZATION TO MATCH TARGET IMS SUPPORTED BY KRIGING METAMODELING
The multi-objective optimization of Equation 1 can be efficiently performed by using the kriging approximation given by Equation 6 when evaluating performance objective F 1 . Note that calculating objective F 2 is computationally trivial. In addition, the approximation error of the metamodel can be incorporated in the objective function definition, leading to the following modification 26 : The motivation for incorporating this error is to improve the robustness of the optimization and to avoid convergence to erroneous solutions due to poor quality of the metamodel. This feature will be further explored in the illustrative implementation.
For solving the multi-objective optimization, a variety of numerical approaches can be used. 42 Here two are considered, one gradient-free and one gradient-based, offering different advantages. The first approach adopts an exhaustive search. 43 A very large number of n bc samples for θ are generated that are close to μ θ (z), and objective functions F 1 and F 2 are calculated. Estimation of objective F 1 in this case leverages the computational efficiency of the metamodel in performing vectorized predictions: the calculations are simultaneously performed for all n bc samples, or using subsets with a lower number of samples depending on the available computational resources (memory can be a problem for vectorizing operation). This greatly reduces computational time for estimating F 1 . The dominant solutions representing the Pareto front can be then readily identified by comparing the values for the two objectives. The challenge in this case is that the value of n bc needs to be large to obtain an adequate representation of the Pareto front. The advantage is that vectorized calculations can be used for the metamodel predictions. For the stochastic ground motion model described in Appendix A, the samples for θ can be generated by obtaining samples for v from Gaussian distribution N(μ(z),Σ) and then transforming these to samples for θ through the inverse of Equation A6. This guarantees that samples will correspond to lower values for F 2 [ie, are close to μ θ (z)] and can therefore emerge as dominant solutions. Note that since this approach does not use gradients, it can seamlessly accommodate the alternative objective functions for F 1 given by Equations 3 and 4.
The second optimization approach is a gradient-based one. The epsilon constraint approach 44 is specifically adopted because of its ability to explicitly define the value for one of the objectives. This method converts the multi-objective optimization problem to a set of single-objective constraint optimization problems with different constraint bounds ε r . Through systematic variations of this constrain, different Pareto optimums can be obtained. Here objective function F 2 is adopted as objective and F 1 as constraint. This allows identification of the stochastic ground motion model that provides a specific compatibility with the target hazard (the prescribed constraint). The range of values ε r of interest is determined by identifying first the anchor points of the Pareto front, corresponding to the minimum of objective functions F 1 (θ | z) and F 2 (θ | z) (unconstrained single-objective optimizations). Evidently, optimization for F 2 (θ | z) yields solution μ θ (z). The range for ε r corresponds then to [min F 1 F 1 (μ θ (z) | z)], and a number of different constraint values can be considered, with the exact number depending on the desired resolution of the front. For each such value, the single-objective, constrained optimization is solved This optimization problem is not convex, and a gradient-based approach appropriate for constrained global optimization problems needs to be adopted. This is accomplished through an approximate multistart approach that addressed both the potential existence of multiple local minima as well as challenges associated with identifying feasible starting points for the gradientbased approach. Initially, a large value n init of trial solutions for θ is examined, then the solutions corresponding to lower values of F 2 (θ | z) while satisfying constraint ε r for F 1 (θ | z) are taken as initial points for a gradient-based optimization. The latter is achieved through a SOL solver implemented through the TOMLAB optimization environment. 45 The same candidate solutions are used for all values of ε r (no need to repeat this step). Evaluation of F 1 (θ | z) over the n init candidate solutions is vectorized, so has small burden, whereas the efficiency of the gradient-based optimization is improved by obtaining analytically the gradients for both objectives. For objective F 1 , the components of the gradient vector are where the partial derivatives inside the brackets can be readily obtained through the metamodel. 26 For the specific objective function F 2 that will be used in the illustrative example later, given by Equation 5, the components of the gradient vector are calculated by the chain rule as where φ(.) stands for the standard Gaussian pdf, the partial derivative ∂v i /∂θ i was calculated by differentiating Equation A6, and the partial derivatives ∂F 2 /∂v i correspond to the components of the gradient row vector The challenge for this optimization approach is that a gradient-based step needs to be repeated multiple times (for each different value of ε r ) and cannot leverage vectorized calculations for the metamodel because the metamodel is separately utilized for each sequential objective function evaluation. The advantage is that gradient information can improve computational efficiency and that the optimization can be performed for specific values of ε r . This allows the identification of a specific part for the Pareto front if desired, for example, the front that corresponds to specific levels of compatibility to the chosen IMs.
Once the Pareto front has been identified, a dominant solution can be adopted using any desired criterion, for example, the solution that provides a specific compromise between the two objective functions. This will be further discussed in the illustrative implementation.

| ILLUSTRATIVE IMPLEMENTATION
The illustrative implementation considers the stochastic ground motion model developed by Rezaeian and Der Kiureghian 15 and reviewed in Appendix A. For the target IMs, GMPEs used in the Western US are considered here, 9,[46][47][48] whereas the suggestions by Kaklamanos et al 49 were adopted to estimate unknown inputs for some of the GMPEs. As target, IM predictions from individual GMPEs as well as the average of their predictions will be adopted later. Note that the latter still provides a single target IM for each structural period examined. All computations are performed in a quad-core 3.0 GHz Xeon processor with 16 GB of RAM.  30 and discrete {0,1} for F. These ranges correspond to the ones for which the initial predictive relationships 15 were developed. The samples for the SVM-based characterization of X are chosen as n s = 10000 with δ d = 9 (latter corresponding to radius of 3 standard deviations away from mean). The adaptive DoE discussed in Section 3 is implemented with n 1 = 600 and n a = 300. Three different accuracy criteria are selected, with associated coefficient of determination (averaged over all outputs) 0.95, 0.97, and 0.985. This leads to number of support points 1500, 3000, and 4500, respectively.

| Details for metamodel development
For generating a total of 900,000 time histories and performing the required 19,800,000 simulations to develop the database for the metamodel for the elastic responses, close to 600 CPU hours were required. Although this computational burden is significant, it should be stressed that it corresponds to an initial only overhead of the approach. Once the metamodel is developed, it can be then used for any required predictions since the established accuracy is high. This large burden should be also attributed to the large number of white noise sequences (200), periods (22), and the wide seismicity range examined. The former provides high accuracy for the calculation of the relevant statistics, whereas the latter two support a wide applicability of the developed metamodel, as it can provide accurate predictions for all the responses of interest for the considered stochastic ground motion model (covers its range of applicability) and GMPEs (covers all the periods addressed). The burden for a metamodel that considered smaller number of stochastic sequences or constrained seismicity or period ranges would be drastically decreased.
Estimation of metamodel response for 10000 samples requires 4.9 s, 8.8 s, and 14.05 s for the metamodels with 1500, 3000, and 4500 support points, respectively. When the metamodel prediction error variance is not computed, the corresponding times are 3.6 s, 7.5 s, and 10.7 s, respectively. Note that adoption of larger value of samples prohibits efficient vectorization of operations for the n = 4500 points because of memory restrictions. The calculation of objective function F 1 along with its gradient requires 0.031 s, 0.052 s, and 0.091 s for the metamodels with 1500, 3000, and 4500 support points, respectively. When the metamodel prediction error variance is not included in the calculation, the corresponding times are 0.008 s, 0.0107 s, and 0.015 s, respectively. Comparison of these computational times shows that (i) increase of the support points has a considerable effect on computational efficiency, (ii) vectorization of calculations provides significant benefits, and (iii) inclusion of the prediction error variance in the calculations increases the computational burden, especially when gradient information needs to be obtained. All these aspects should be taken into account when choosing computational details for the optimization problem.

| Comparison of optimization approaches
The focus is first placed on the numerical solution of the optimization problem. The target used in this subsection, and also in the next one 5.3, corresponds to structural periods T s = [0.4 0.5 0.75 1.0 1.5 2.0]s and IM described by the average of the considered GMPEs. 9,46-48 Weights are chosen as γ i = 1/Y i (z), so that objective function is expressed in terms of the relative error. In other words, no specific structural period is prioritized in evaluating the match to the target IMs. Three reference seismicity scenarios are examined in this section, corresponding to M = 6-R = 20 km, M = 7.8-R = 30 km, and M = 7-R = 40 km, for a strike-slip fault (F = 0) and V s,30 = 800 m/s. The first two scenarios correspond to cases where the unmodified stochastic ground motion model does not provide an adequate match to the target GMPEs 15 (significantly over predicts for the former, moderately under predicts for the latter) and the last to a case where the unmodified model facilitates a good match. Because of space limitations, discussion will focus around specific cases of interest. The metamodel with 4500 support points is used in this section. Comparison across metamodels with different number of support points (and therefore accuracy) will be discussed in the next section. Figure 2 presents illustrative swarms of candidate solutions in the objective space from the exhaustive search using n bc = 50000 samples for seismicity M = 6-R = 20 km and M = 7-R = 40 km, in the former case examining the case with and without metamodel error in the calculation of objective F 1 . For the second objective, results are reported with respect to ffiffiffiffiffi ffi F 2 p , which corresponds to the distance between the initial and the modified predictive relationships (not the squared distance) and offers a better normalization for the results in the comparison. The solutions located at the left boundary of the swarms in Figure 2 correspond ultimately to Pareto optimal solutions since there is no other solution that can simultaneously improve both performance objectives. It is evident from these swarms that modification of the predictive relationships can indeed facilitate a difference in the IM match (check the range of F 1 values obtained), whereas for smaller F 1 values, the candidate solutions deviate more from the rest of the swarm. This means, ultimately, that there are fewer model configurations that can provide a good match to the target IMs (small F 1 values). When the unmodified model is closer to the target IMs (M = 7-R = 40 km seismicity scenario), higher compatibility to these IMs can be obtained through the modification of the predictive relationships [compare case (C) to the other two], whereas addition of the metamodel error [compare case (B) to (A)] increases function F 1 , especially for smaller F 1 values. Figure 3 examines different approaches for the solution of the multi-objective problem for the reference seismicity case of M = 6-R = 20 km (similar trends hold for the other cases). Part (A) shows the Pareto fronts identified through the exhaustive search using three different n bc values. Results are reported herein with respect to ffiffiffiffiffi ffi F 1 p and ffiffiffiffiffi ffi F 2 p because this facilitates an easier comparison (differences of extreme values easier to discern). Only 10 representative solutions are shown, and not the entire front. It is evident that minor differences only exist between the identified fronts for different n bc values, and such differences occur primarily for small F 1 -large F 2 combinations. This comparison shows that a value of n bc around 200,000 to 400,000 should be considered as sufficient for efficiently identifying the front. For the remaining of the manuscript, results for the exhaustive search will be presented for a value of n bc equal to 400,000. Figure 3B then compares the Pareto fronts obtained by the exhaustive search of the gradient-based (epsilon constraint) approaches for the case that the metamodel error is considered or not in the objective function F 1 . For the epsilon constraint, n init is taken as 10000, whereas the gradient-based optimization is performed using a couple only different initial points. For small F 1 constraints, some challenges were encountered in converging to an admissible solution. This should be attributed to the trend identified in Figure 2; for such values, a smaller number of candidate solutions exist that can satisfy the constraint, and a potential increase in n init might be needed to identify feasible initial points for the gradient-based optimization. The comparison in Figure 3B shows that the two approaches identify similar Pareto fronts, with the gradient-based optimization converging to suboptimal solutions for lower F 1 values, evidently because of the existence of local minima with greater differences in achieved performance. For such performance ranges (ie, with an already good match to the target IM, as indicated by the lower F 1 value), these differences might be unimportant. With respect to the computational burden, the exhaustive search requires 15 s per 10000 candidate solutions examined (11 s if metamodel error is disregarded). For the epsilon constraint, the computational cost is 15 s for the initial 10000 trials and 20 s (2 s if metamodel error is disregarded) for each different constraint examined. These comparisons show that both optimization approaches may be considered as adequate and preferred, depending on the application context. Overall, some preference exists for the exhaustive search due to the fact that epsilon constraint method needs identification of an appropriate starting point to avoid convergence to suboptimal local minima. The epsilon constraint approach, though, might be beneficial when a single solution is sought after, the one that satisfies a desired match to the target IM, rather than the entire front. Consideration or not of the metamodel error has no effect on the differences between the optimization approaches. The computational efficiency for the gradient-based search is reduced, as discussed in the previous section, when this error is included. Table 1 provides summary of these results, including some details discussed in the next section.  Figure 4 presents a comparison between the alternative objective function selections for quantifying the discrepancy from the target IM, ie, comparison between F 1 , F 1r , and F 1m . Figure 4A discusses F 1r and Figure 4B F 1m . The exact Pareto front, ie, optimization using F 1r (or F 1m ) and F 2 as objectives, as well as an approximate front are compared. The approximate front is obtained by using F 1 and F 2 as objectives, identifying the Pareto set for θ and then evaluating objective F 1r (or F 1m ) over that Pareto set. As discussed earlier, this optimization is less challenging but evidently identifies a sub-optimal solution. The results in Figure 4 show that there is overall very good correlation between the different objective functions and that the approximate solutions have only small deviations from the optimal front and those only for small values of the first objective. Therefore, use of F 1 as objective may be considered as an adequate surrogate even when the interest is in objectives F 1r or F 1m for quantifying discrepancy from the target IMs.

| Impact of metamodel accuracy
The discussion moves next to the examination of the effect of the metamodel accuracy. This is established by considering additionally the results obtained by using the exact stochastic ground motion model (ie, not relying on metamodel predictions), which represents the measure for evaluating the actual hazard compatibility of the identified ground motion model. The   Figure 5 shows the Pareto fronts identified by using the metamodels with the three different number of support points. Cases with or without the metamodel error in the estimation of objective function F 1 are separately shown. This leads to two different Pareto sets, one without error and one with error, and for each two different fronts are reported, one corresponding to metamodel predictions and one to the use of the exact stochastic ground motion model. Then Figure 6 shows the spectral plot comparisons for the solution (among the Pareto set identified in each case) corresponding to the minimum of F 1 . The period range used in Figure 6 corresponds to the target structural periods T s . In all cases, the exhaustive search is implemented with the same candidate solutions to facilitate a consistency in the corresponding comparisons. Figure 7 then shows Pareto front results for different seismicity scenario, M = 7.8-R = 30 km. Note that for seismicity scenario, M = 7-R = 40 km results (another case discussed in the previous section) are of limited interest since the unmodified ground motion model provides a good compatibility to target IMs. The results show that for the higher accuracy metamodel (4500 support points), good agreement is established between the metamodel predictions and the actual model predictions along the Pareto front, whereas the inclusion of the metamodel error has only a small effect on the identified Pareto front. This is observed in both Pareto fronts (Figures 5 and 7) as well as in the corresponding spectral plots (Figure 6). Different trends are observed, though, for the lower accuracy metamodel (1500 support points). For lower F 1 values, and therefore large F 2 values, the Pareto optimal solutions identified when metamodel error is not included in the problem formulation lead to erroneously modified ground motion models. Based on the metamodel predictions, these models provide a very good match to the target IMs (small predicted F 1 value), but when response is evaluated with the exact model, larger differences are observed from the target IMs. This is particularly evident in the spectral plots shown in Figure 6. Evidently, the respective solutions identified correspond to parameters θ for which the metamodel accuracy is low.  The moderate accuracy metamodel (3000 support points) falls in between the two aforementioned cases, with characteristics that resemble more closely the ones for the high accuracy metamodel.
Note that the seismicity case examined here is ideal for exploring vulnerabilities in the optimization associated with lower metamodel accuracy, as it corresponds to a case at the boundary of the scenarios used to define domain X. Therefore, for larger values of F 2 , the corresponding parameters θ are expected to be close to the boundary of X, where metamodel accuracy is lower. However, even for this challenging case, metamodels with higher accuracy (3000 or 4500 support points) face small challenges, whereas the inclusion of the metamodel error for the lower accuracy metamodel (1500 support points) greatly improves the robustness of the optimization, leading to identified solutions with good agreement between metamodel and actual model. Figure 8 sheds further light into this topic. It focuses on the case examined in Figure 6A, but besides the mean metamodel predictions, it includes the predictions that are 1.5 standard deviations (σ) from the mean based on the estimated error variance. The solution identified when metamodel error is not included in the evaluation of F 1 is associated with a larger anticipated error ( Figure 8A). Note that the actual model is actually even further away than the plotted 1.5 σ. When the metamodel error is included in the evaluation of F 1 , such θ values with large associated σ are avoided because the large σ contributes to larger values for the objective function. This ultimately contributes to identification of solutions with greater robustness, ie, better agreement between metamodel and actual model ( Figure 8B).
Overall, the above discussion shows that metamodels with higher accuracy (coefficient of determination 98%) can be considered as a good surrogate for the proposed optimization, whereas the inclusion of the prediction error greatly improves the   Figure 5 for the metamodel with 1500 support points. For the metamodel predictions, the mean predictions and the predictions within 1.5 standard deviations from the mean are shown robustness of this optimization, avoiding identification of erroneous solutions, even when metamodels with lower accuracy are adopted. Consideration of this error is not necessary, though, when higher accuracy metamodels are used.

| Implementation for different seismicity scenarios
With the computational details ironed out, the discussion moves finally to the IM compatibility established by the proposed modification of the ground motion model. Figure 9 shows results for three seismicity scenarios targeting PSA given by the average of the aforementioned GMPEs for two different ranges for T s : T s = [0.4 0.5 0.75 1.0 1.5 2.0] s and T s = [0.4 0.5 0.75] s. These two different cases are referenced herein as long and short, respectively, period ranges. The proposed approach identified in each case a Pareto front that clearly demonstrates the compromise between the two objectives, with different characteristics in each case, depending on how close the unmodified ground motion model was to the target IM. Choosing a shorter period range for this target facilitates an overall better match; this is anticipated because objective F 1 imposes less strict requirements in terms of IM compatibility (fewer number of components to match).
The question finally arises which point should be selected within the identified Pareto set. Various approaches have been proposed in the greater multi-objective optimization literature for making this choice. 50 Perhaps the most common one is to select the solution that has the smallest normalized distance from the utopia point, defined as the point in the Pareto front (ie, objective function space) that corresponds to the minimum of the two objectives across the front. This utopia point represents the best but unachievable performance [minF 1 minF 2 ]. Normalization is typically established with respect to the maximum of each performance objective across the Pareto front, leading to distance metric Instead of F 1 and F 2 , this can be implemented for ffiffiffiffiffi ffi F 1 p and ffiffiffiffiffi ffi F 2 p due to the better normalization properties. This point is identified in all cases in Figure 9. Another choice would have been to choose the solution that satisfies a predetermined threshold for the match to the targeted IMs. Perhaps this is better set with respect to objectives F 1r or F 1m rather than F 1 . ffiffiffiffiffi ffi F 1 p smaller than 0.15 or 0.05. These three cases are denoted, U t , C l , and C s , respectively. The thresholds for C l and C s modifications were chosen so that to reflect medium and small, respectively, incompatibility to the target hazard. In addition, results for the unmodified model are presented denoted U n . Figure 10 shows plots for (i, first row) ffiffiffiffiffi ffi F 1 p for U t and U n ( ffiffiffiffiffi ffi F 1 p is constrained for the other two cases) and for (ii, second and third rows) ffiffiffiffiffi ffi F 2 p for U t , C l , and C s ( ffiffiffiffiffi ffi F 2 p is zero for U n ). To better demonstrate the differences, results are presented separately for U t (second row) and for the pair C l and C s (third row) in the latter case. The three different columns in the figure correspond to three different implementation cases: target IM given by the average of the aforementioned 4 GMPEs for both (A) the long and (B) the short period ranges as well as (C) target IM given only by one only GMPE 47 for the long period range. These scenarios are denoted herein as SC 1 , SC 2 , and SC 3 , respectively. Figure 11 shows spectral plots for a smaller (due to space limitations) selection of seismicity scenarios, defined by combinations of M [6. 2, 6.8, 7.4, 8] and R [30,60,90] km, for SC 1 . For each of the 12 M-R combinations, the curves corresponding to the target IM, the unmodified model, and the predictions by the three aforementioned model modifications are shown to facilitate comparisons. Finally, Figure 12 shows for all examined seismicity scenarios the model parameters θ for the unmodified ground motion model as well as for the modified model corresponding to the Pareto point with smallest distance from the utopia point (U t case) for SC 1 . Note that some of the curves shown in these figures have nonsmooth characteristics. This should be attributed to multiple facts: (i) a discrete representation of the Pareto front is identified, rather than the actual Pareto front; (ii) problem has multiple local minima as discussed earlier, especially for larger F 2 FIGURE 10 Results for ffiffiffiffiffi ffi F 1 p and ffiffiffiffiffi ffi F 2 p for unmodified ground motion model (U n ) and modified ground motion model corresponding to three different selection criteria: Minimum distance form utopia point (U t ) and value ffiffiffiffiffi ffi F 1 p smaller than 0.15 (C l ) or 0.05 (C s ). Implementations in the different columns correspond to A, long and B, short period ranges for matching to the average considered GMPEs and C, long period ranges for match to GMPE 47 [Colour figure can be viewed at wileyonlinelibrary.com] values; and (iii) algorithms used for the optimization have characteristics of stochastic search, well known to lead to discontinuous results.
The unmodified model U n does not provide a good match to the target IMs for the entire seismicity range, with trends observed in first row of Figure 10 and in the spectral plots in Figure 11 being similar to the ones reported by Rezaeian and Der Kiureghian, 15 ie, greater challenges in lower moment magnitudes and combination of higher magnitudes and larger distances. Figure 10 shows that the proposed modification (cases U t , C l and C s ) improves this match, establishing a balance between F 1 and F 2 , with the characteristics of the balance depending on the criteria for selection of the final model among the Pareto optimal solutions (ie, which specific case is chosen). When the unmodified model has larger discrepancies from the target IMs, then the modifications lead to larger values for F 2 , but still successfully identify models, independent of the implementation case, that provide an improved match to the IM target. The U t modification identifies a model with moderate discrepancy from the unmodified one, corresponding to values of F 2 in the range of 0.3-1, whereas the two other modification approaches, C l and C s , identify models with greater variability across the different seismicity scenarios. For scenarios in the range of M = 7-7.5, the unmodified model provides a good match to the target IMs and therefore modification of it provides limited advantages. This is perhaps better captured by the C l case. This discussion shows that selection of the Pareto optimal model based on a targeted accuracy to the GMPEs, ie, value for F 1 below a certain threshold as in the C s and C l cases, provides a more rational selection for the final model as it allows a more direct recognition of the seismicity ranges where modification is not truly required. However when that threshold is selected small (C s case) and the unmodified model has larger discrepancy from the target IM, then the modification leads to identification of a model with big differences from the original one (large F 2 values). This model will typically be far away from the U t case and will belong in a steep part of the Pareto front, meaning that small improvements in F 1 come at a large increase of F 2 (check the Pareto fronts shown in previous figures). A multilevel selection criterion seems therefore more appropriate: select the Pareto optimal solution that satisfies a certain accuracy threshold for the target IM (target), the unmodified ground motion model (U n ) and modified ground motion model corresponding to three different selections criteria: minimum distance form utopia point (U t ) and value ffiffiffiffiffi ffi F 1 p smaller than 0.15 (C l ) or 0.05 (C s ). Implementation scenario shown corresponds to matching to the average considered GMPEs and long period range for F 1 unless this solution leads to a greater F 2 value than the Pareto optimal solution with minimum distance from the utopia point. If the latter happens, then select the Pareto optimal solution with minimum distance from the utopia point.
With respect to the different implementation cases shown in Figure 10, selection of a shorter period range [compare cases (B) and (A)] promotes, as also identified earlier, an easier match to the target (smaller overall values for F 1 and F 2 ). This demonstrates the importance of carefully selecting the target IMs, as this selection affects the ability to match this target, and of the availability of a framework, as the one developed here, that allows you to do so. The exact selection of the target [compare cases (A) and (C)] does not impose any additional constraint in the optimization implementation, although the results evidently change.
The spectral plots in Figure 11 provide further validation for the observed accuracy of the different examined cases, showing the decomposition of the overall F 1 match to the different structural periods. The ground motion model modification greatly improves the match to the target spectral curves, with C s modification always providing a better match compared to the C l one. The U t modification exhibits a varying behavior, frequently in between the C l and C s cases, in other instances outside their envelope, with characteristics that ultimately depend as discussed above (based on Figure 10) on the ability of the unmodified model to provide an adequate match to the target IMs (ie, how much is the modification truly needed). Of course the proper evaluation of the proposed modifications comes from examining both F 1 and F 2 values, as detailed above (discussion focusing on Figure 10).
Finally, the results for the ground motion model parameters in Figure 12 show that the model modification leads to similar trends as observed for the unmodified model. This is the direct results of incorporating the difference between these parameters as an objective in the problem formulation (objective F 2 ) rather than merely as a constraint. Parameters I a , ζ f , and ω ' show bigger variability compared with their initial values. This should be attributed to a greater sensitivity with respect to them of the resultant ground motions.
The overall discussion shows the importance of the established framework: once the initial metamodel is developed, through the adaptive guidelines established here, it can support the efficient identification of ground motion models that (i) match any desired IMs for any chosen period range while (ii) maintaining a small deviation from the initial predictive relationships. This can be seamlessly repeated for any seismicity scenario. The final ground motion model modification can be chosen based on the criteria discussed earlier.

| CONCLUSIONS
The modification of stochastic ground motion models to match target IMs for specific seismicity scenarios was presented in this article. This was formulated as a multi-objective optimization problem with first objective (F 1 ) quantifying the discrepancy between the ground motion predictions and the target IMs, and the second objective (F 2 ) corresponding to the discrepancy between the new model characteristics and the model characteristics suggested by existing predictive relationships (ie, the unmodified model). The second objective explicitly incorporates in the modification process physical characteristics and parameter correlations described in the initial ground motion model. The developed framework facilitates a match to any desired IM or to a collection of them, for example, spectral accelerations over a period range, for any chosen seismicity scenario. Repeating process for different seismicity scenarios can then facilitate the development of a suite of models that can support comprehensive seismic risk assessment. Computational efficiency was achieved by adopting a metamodel for approximating the median ground motion model predictions for the targeted IMs. Although the upfront cost for development of this metamodel is significant, once established, it can be used to support a highly efficient multi-objective optimization. Gradient-based and gradient-free approaches were discussed for the latter, whereas to reduce the computational burden for the metamodel development, an adaptive design of experiments was proposed for selecting the database informing the metamodel.
The framework was demonstrated in an illustrative example considering a recently developed record-based stochastic ground motion model and IMs described through ground motion prediction equations. It was shown in the context of this example that the metamodel-aided optimization can support an accurate identification of the Pareto front of dominant solutions, provided that the metamodel accuracy is significant high, and that inclusion of the metamodel error in the optimization formulation greatly improves the robustness of this optimization, avoiding the identification of erroneous solutions. Comparisons between the two optimization approaches showed that the gradient-free one demonstrates overall preferable attributes because the gradient-based one might converge to suboptimal local minima, especially for lower F 1 values. However, the gradient-based one provides greater relative efficiency when identification of a single solution, rather than of the entire front, is desired. Still, an adequate representation of the overall Pareto front can be obtained in as little as 2 minutes using the blind-search, gradient-free optimization, which should be considered as an acceptable computational burden. Different approaches can be then used to select the final ground motion model, eg, the Pareto optimal point that has the minimum distance from the utopia point or the point that satisfies a desired compatibility to the target IMs (ie, value of F 1 smaller than a threshold). Implementation for a range of seismicity scenarios showcased the advantages of the proposed framework: small modifications of the original ground motion model (small to moderate F 2 values) provided significant improvement for the match to the target IM (F 1 value) for seismicity ranges where the unmodified model faces challenges in matching the target IMs. With respect to the selection of the final ground motion model, the following suggestion is provided after carefully examining the various trends observed: select the Pareto optimal solution that satisfies a certain accuracy threshold for F 1 unless this solution leads to a greater F 2 value than the Pareto optimal solution with minimum distance from the utopia point. If the latter happens, then select the Pareto optimal solution with minimum distance from the utopia point.