Hazard‐compatible modification of stochastic ground motion models

A computationally efficient framework is presented for modification of stochastic ground motion models to establish compatibility with the seismic hazard for specific seismicity scenarios and a given structure/site. The modification pertains to the probabilistic predictive models that relate the parameters of the ground motion model to seismicity/site characteristics. These predictive models are defined through a mean prediction and an associated variance, and both these properties are modified in the proposed framework. For a given seismicity scenario, defined for example by the moment magnitude and source‐to‐site distance, the conditional hazard is described through the mean and the dispersion of some structure‐specific intensity measure(s). Therefore, for both the predictive models and the seismic hazard, a probabilistic description is considered, extending previous work of the authors that had examined description only through mean value characteristics. The proposed modification is defined as a bi‐objective optimization. The first objective corresponds to comparison for a chosen seismicity scenario between the target hazard and the predictions established through the stochastic ground motion model. The second objective corresponds to comparison of the modified predictive relationships to the pre‐existing ones that were developed considering regional data, and guarantees that the resultant ground motions will have features compatible with observed trends. The relative entropy is adopted to quantify both objectives, and a computational framework relying on kriging surrogate modeling is established for an efficient optimization. Computational discussions focus on the estimation of the various statistics of the stochastic ground motion model output needed for the entropy calculation.


| INTRODUCTION
The relevance of techniques that model acceleration time-series of seismic events has increased during the past decades due to the growing popularity of simulation-based probabilistic seismic risk assessment [1][2][3] and performance-based earthquake engineering. [4][5][6] Though the most popular methodology for performing this task is the selection of real (ie, recorded from past events) ground motions, 7-10 potentially scaled based on a target intensity measure (IM), an alternative philosophy is the use of simulated ground motions. 11,12 A specific modeling approach for the latter which has been steadily gaining increasing attention by the structural engineering community, [13][14][15] is the use of stochastic ground motion models. [16][17][18][19][20][21][22] These models are based on a parametric description of the spectral and temporal characteristics of the excitation, with synthetic time-histories obtained by filtering a stochastic sequence through the resultant frequency and time domain modulating functions. The parameters involved in this description are related to seismicity (for example moment magnitude and rupture distance) and site characteristics (for example shear wave velocity for soil profile) through predictive models/relationships. 16,20 Sample ground motions for a specific seismicity scenario and site can be generated by determining the parameters of the stochastic ground motion model through these predictive relationships and by utilizing a sample stochastic sequence. This approach may ultimately support a comprehensive description of the seismic hazard, 6,23 and its essential component is the predictive models relating seismicity/site characteristics to ground motion model parameters.
The two main methodologies for establishing such stochastic ground motion models are record-based and physicsbased approaches. Record-based models (also known as site-based) are developed by fitting a preselected "waveform" to a suite of recorded regional ground motions. Regression analysis is used for establishing the predictive models, which leads to a probabilistic characterization described by mean predictions along with an associated variability. 16,19 On the other hand, stochastic physics-based models rely on physical modeling of the rupture and wave propagation mechanisms. 20,21 The predictive relationships in this case are typically described by deterministic models that represent the underlying mean physical properties, though approaches exist for addressing variability in these properties. 13,24 Emphasis in this study will be on record-based models, though the techniques discussed can be extended to any type of stochastic ground motion model, with the assumption that the corresponding predictive models are characterized by both a mean prediction and an associated variance.
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 models, compatibility to the seismic hazard for specific structures and sites is not necessarily obtained. This hazard is typically characterized through a Probabilistic Seismic Hazard Analysis, 25 for example through de-aggregation that identifies the seismicity scenarios, described through relevant seismicity characteristics, mainly the moment magnitude and rupture distance, that have the largest contribution to the hazard for a specific structure. Essential part of Probabilistic Seismic Hazard Analysis are ground motion prediction equations (GMPEs). GMPEs provide predictions, as function of seismicity characteristics, for both the median and the dispersion of IMs, determining ultimately the conditional hazard for seismic events corresponding to these characteristics.
Recognizing the importance of matching stochastic ground motion models to some target IM the modification of the predictive relationships for accommodating such a match was recently examined. 26,27 This was posed as an optimization problem with objective to minimize the discrepancy between the median predictions of the ground motion model and the spectral acceleration estimates of GMPEs for different period ranges (or, more broadly, a range of target IMs). Appropriate constraints were set in this optimization to preserve desired ground motion physical characteristics, so that unrealistic time-histories are not obtained. This formulation was first introduced in Scherbaum et al, 26 with Vetter et al 27 offering a computationally efficient approach for performing the associated optimization, leveraging surrogate modeling principles. Recent work by the authors 28 extended this framework by addressing a critical shortcoming, the fact that physical descriptors of the resulting acceleration time-series were incorporated in the optimization merely as constraints, something that required significant experience in ground motion characterization for proper definition of the optimization problem. The shortcoming was addressed 28 by introducing a bi-objective optimization problem, transforming the aforementioned constraint for the physical characteristics of the resultant ground acceleration time-series to an explicit objective. The proposed tuning was performed for specific seismicity scenarios and identified the ground motion model that achieves the minimum modification of the existing predictive relationships that will yield the desired compatibility with the target IM.
All the aforementioned studies focused, however, on the mean model characteristics and associated hazard. Optimization utilized only the mean of the predictive relationships of the stochastic ground motion model, whereas, more importantly, match only to a target IM vector (taken to represent the mean hazard) was investigated, ignoring any variability in the IM predictions. The latter is an important constraint because for seismic risk assessment applications hazard compatibility is expressed in terms of both the mean and dispersion of the target IMs. 25 From a practical standpoint, capturing the actual variability of target IMs is essential to capture extreme structural response values and, therefore, in properly assessing the likelihood of consequences due to such extreme seismic demand values.
The current study extends approach 28 to (1) match the prescribed conditional hazard (not simply mean IMs) for a specific site and structure (or range of structures) while (2) preserving desired trends and correlations in the physical characteristics of the resultant ground acceleration time-series, including consideration of the variability of these characteristics. This is again formulated as a bi-objective optimization problem. The first objective is to minimize the discrepancy between the statistics (mean and dispersion) of the outputs for a suite of acceleration time-histories obtained from the ground motion model and the target IM statistics for a given seismicity scenario. The second objective is to establish the smallest deviation from the model characteristics suggested by existing predictive models, examining both the mean and the variance of the model parameters. Both objectives are expressed as comparison between probability distributions, and the relative entropy is adopted to quantify them. This setting creates a fundamental difference to the previous study 28 with respect to both the optimization characteristics (alter both mean and variance of predictive models) as well as the goal (match to mean and dispersion for conditional hazard). In Tsioulou et al, 28 the ground motion model was tuned so that outputs from a single parametric description of the ground motion model match a target IM vector for each seismicity scenario. The goal of the current study is to produce an ensemble of ground motion models whose output statistics yield the desired compatibility with the hazard (IM mean and dispersion) for that scenario. For efficiently solving the multi-objective optimization problem, a surrogate modeling approach is adopted 27,28 for approximating the desired IMs for specific values of the ground motion model parameters. Emphasis is placed here on the efficient estimation of the response statistics for the modified ground motion model output, leveraging Monte Carlo techniques. This requires further extension of the surrogate modeling framework, compared with the approach adopted in Tsioulou et al, 28 for facilitating this estimation. Different assumptions are also examined for the evaluation of the entropy for the first objective. The corresponding bi-objective optimization is finally solved using a random search approach. The novelty of the current work stems from both the fundamentally different theoretical framework for formulating the stochastic ground motion modification as well as the computational advances required for efficiently calculating the response statistics for the first objective and performing the associated optimization.
In the next section, the general problem of developing simulated ground motions compatible with target IM distributions is defined, and then specific aspects of the framework are discussed in detail.

| Preliminaries and baseline predictive relationships formulation
The foundation of the problem formulation is the same as in Tsioulou et al. 28 A stochastic ground motion model is considered that provides acceleration time-histories € a tjθ; w ð Þby modulating a discretized Gaussian white-noise sequence, w, through appropriate time/frequency functions that are parameterized through the n θ -dimensional model parameter vector θ ¼ θ 1 θ 2 … θ n θ ½ ∈ℜ n θ . θ completely defines the parametric description of the model, ie, along with w facilitates the generation of a time-history € a tjθ; w ð Þ, 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. It should be noted that θ corresponds typically to a low-dimensional vector and w to a high-dimensional sequence.
Predictive models/relationships are utilized to relate θ to seismicity and local site characteristics, such as the fault type F, the moment magnitude, M, the rupture distance, R, and the shear wave velocity in the upper 30 m of soil, V s30. 16,20 The vector of these characteristics is denoted as z. For record-based models, the standard approach for development of these predictive relationships 16,29 relies on first matching the waveform of recorded ground motions (ie, identify first θ for each of the recorded ground motions in a given database) and then carrying out a regression to relate θ to z. Typically, this is performed by first transforming problem to the standard Gaussian space through a nonlinear mapping for each component θ i. 16 The transformed Gaussian vector is denoted v. Approach ultimately leads to a Gaussian probability model v~N(μ r (z),Σ r ) with mean μ r (z) and covariance matrix Σ r . Note that the latter is independent of z. The notation N(a, b) stands for Normal distribution with mean a and covariance b whereas notation c~d stands for random variable c following distribution d. The resultant probability model for θ is denoted p(θ|μ r (z),Σ r ) and determines the predictive model for the ground motion model parameters. Note that a similar description can be readily established for physics-based models, 28 with the uncertainty characterization stemming from an explicit treatment of the epistemic uncertainties associated with the physics-based formulation. 13,24 The predictive model for θ will be denoted herein as p(θ|μ r (z),Σ r ), with the understanding that it is not necessarily constrained to models established through regression analysis, rather simply can be parameterized by quantities μ r (z) and Σ r , with μ r (z) representing the mean predictions, and Σ r the variability of these predictions.

| Modification of predictive models for hazard matching
As discussed in the introduction, the formulation of the predictive model for θ provides synthetic ground motions whose statistics (mean and dispersion) of output IMs do not necessarily match the intended hazard for specific structures and sites. For accommodating such a match, a modification of the existing predictive model for θ is proposed for specific seismicity scenarios defined through z, with objective to get a suite of acceleration time-series for that scenario whose (1) mean and dispersion match a target IM mean and dispersion vectors, while (2) maintaining similarity to the predictive relationships already established for the model. Equivalently, this can be viewed as identifying the updated probabilistic model p(θ|μ,Σ) that is closest to the established model p(θ|μ r (z),Σ r ) and also matches the intended conditional hazard. μ and Σ represent the updated parametric description for the probability model of θ. In the context of recordbased models, these correspond, respectively, to the mean vector and covariance matrix for v. The IM vector may include different response quantities of interest, 28 for example direct characteristics of the ground motion, such as peak ground acceleration or elastic and inelastic spectral responses for different periods of a single degree of freedom oscillator. The conditional hazard for most of these IMs may be described through a GMPE. 30,31 Because the proposed modification refers to the conditional hazard, for simplifying terminology for the remainder of the paper, the description as "conditional" will be removed: term hazard corresponds to conditional hazard.
To formalize these concepts mathematically, let, Y i (z); i = 1, . …n y denote the IMs of interest. The target hazard for them is quantified through a probabilistic description, denoted by p t (ln(Y i )| z). As is common in earthquake engineering and without loss of generality, the statistics are assumed here to be determined for the logarithm of the IM. To better align approach with current GMPE standards a lognormal underlying model is assumed, leading to and σ 2 i z ð Þ corresponding to the mean and variance, respectively, of the logarithmic IM. Note, however, that the computational framework can support any probabilistic IM description p t (ln(Y i )| z), not constrained to one provided by GMPEs or defined through a lognormal probabilistic model. Also, the superscripts/subscripts t and g (the latter defined in the next paragraph) are used herein to distinguish between target prediction and prediction facilitated through the stochastic ground motion model.
To quantify the hazard predictions through the stochastic ground motion model, let Y g i θ; w ð Þdenote the estimate for Y i established through this model for specific values of the model parameter vector θ and a specific white noise sequence w (ie for a specific ground motion time-history € a tjθ; w ð Þ). Y g i θ; w ð Þ will be referenced herein as response output of the ground motion model. The statistical characterization for ln(Y i ) through the stochastic ground motion model for the updated parametric description is denoted by p g (ln(Y i )| μ, Σ) and using the total probability theorem is equal to: where p(w) is the probability distribution for the stochastic sequence w, which by definition is independent of p(θ|μ,Σ).
In deriving the second equality in Equation 1 the fact that the stochastic ground motion model is completely defined by pair θ, w was also used. This simplifies p ln Y g Þbecause knowledge of μ and Σ is redundant if θ is also known. Note that Y g i is itself a random variable with randomness in its description stemming from both w and θ. In this context, Y g i θ; w ð Þ may be viewed as a realization of the random variable. Distribution p g (ln(Y i )| μ, Σ) can be approximated through kernel density estimation (KDE) as discussed in Appendix A. To further simplify framework, the same hypothesis as for the hazard can be adopted, assuming a lognormal distribution, leading to ln , where the mean and variance for ln(Y i ) g , considering the variability in both the low-dimensional θ and high-dimensional w, are: with E[.], Var[.] denoting the expectation and variance operators, respectively. The functional dependence of the statistics of ln (Y i (z)) g on μ and Σ is explicitly noted herein to facilitate an easier understanding of the ground motion model modification framework. The log-normal assumption will be primarily used for p g (ln(Y i )| μ, Σ), and its validity will be examined in the illustrative example later by comparing to the KDE-based estimation. The hazard compatible modeling corresponds to modification of the probability model for θ, ultimately of the parametric description defined through μ and Σ, and is formulated as multi-objective optimization problem with two competing objectives The first objective F p1 corresponds to the weighted discrepancy of the target seismic hazard to the hazard predicted through the ground motion model, ie, to a comparison between p g (ln(Y i (z)| μ(z), Σ) and p t (ln(Y i )| z). The relative entropy, a popular measure to quantify differences between distributions, 32 is utilized as metric for this discrepancy, given by This leads to definition of F p1 as with γ i corresponding to the weights prioritizing the match to different IM components (eg, spectral accelerations at different structural periods). For the assumption of lognormal distribution for Y g i , this simplifies to with closed-form solution 32 : Objective F p2 measures that discrepancy between the initial predictive model for θ, p(θ|μ r (z),Σ r ), and the modified one, p(θ|μ,Σ). The relative entropy is utilized again as measure to quantify differences: Because the relative entropy is invariant under a coordinate transformation for θ, 32 the comparison can be established in the transformed Gaussian space, leading to with the latter expression readily evaluated 32 as where tr[.] and det [.] stand for trace and determinant, respectively. Objective F p1 enforces hazard compatibility (mean and dispersion), whereas objective F p2 guarantees compatibility of the physical characteristics of the resultant ground motions with the regional trends observed in recorded ground motions. Solution of the multi-objective optimization of Equation 4 ultimately leads to a Pareto set of dominant solutions {(μ p, Σ p ); p = 1, …,n p } that expresses a different compromise between the competing objectives F p1 and F p2 . A solution is characterized as dominant (or Pareto optimal) and belongs in the Pareto set if there is no other solution that simultaneously improves both objectives F p1 and F p2 . The representation of the Pareto set in the performance objective [F p1 ,F p2 ] space, {[F p1 (μ p, Σ p ),F p2 (μ p, Σ p )]; p = 1, …,n p } is termed as the Pareto front. Illustrations of such Pareto fronts are included in the example discussed later. One can eventually select a model configuration from the identified Pareto set that yields the desired hazard compatibility without deviating significantly from regional ground motion characteristics. This will be further discussed in the illustrative implementation. Note that a Pareto set for optimization of Equation 4 always exist apart from the extreme case that the original stochastic ground motion model provides an exact match to the hazard, ie, F p1 (μ r (z), Σ r | z) = 0. In that case, the entire set corresponds to a single point (μ r (z), Σ r ).
Identifying the Pareto set is, however, challenging because the computational burden in evaluation of objective F p1 is significant, requiring calculation of the multidimensional integrals of Equations 2 and 3 (or samples needed for the KDE approximation discussed in Appendix A) which can be only performed numerically and entails hundreds estimates of the response output of the ground motion model. To facilitate an efficient optimization that can be repeated for any desired seismicity scenario z, a surrogate modeling (ie, metamodeling) approach is adopted, specifically selecting kriging as metamodel due to its proven capability to approximate well even complex functions. 33 As input for the metamodel, the low-dimensional vector θ is chosen. This choice corresponds to the smallest possible dimension for the metamodel input, something that can greatly enhance accuracy, 33 a very important consideration because metamodel will be eventually used for optimization, and larger metamodel errors can lead to identification of suboptimal solutions. Alternative choices were to additionally include w in the metamodel input definition which is completely impractical (due to high-dimensionality of w), or to use directly pair {μ,Σ} as the metamodel inputs (and statistics of Equations 2 and 3 as outputs) which significantly increases, however, the input dimension. Under this selection of θ as metamodel input, the metamodel output corresponds to the response output statistics considering variability with respect to w only (ie, conditional on θ statistics). Here, statistics refers to all quantities needed to eventually calculate F p1 . The variability with respect to θ in quantifying the hazard (and ultimately objective F p1 ) will be addressed through the approach discussed in Section 4.
This creates a similar setting as in Tsioulou et al 28 for the metamodel formulation, with the additional requirement, however, to eventually incorporate variability stemming from θ in the overall framework, so that variability in both θ and w is explicitly considered. The kriging metamodel development is briefly reviewed in the next section, and the details of the kriging-aided optimization problem are discussed in Section 4. A schematic of the overall optimization approach is provided in Figure 1 on the next page.

| KRIGING METAMODEL DEVELOPMENT
As discussed in the previous section, θ is chosen as metamodel input. A modification of this input should be further adopted if the relationship between some components of θ and the response output Y g i θ; w ð Þ is explicitly known. This is true for stochastic ground motion models that include a scaling parameter, denoted θ s herein, that directly controls the amplitude of the excitation. For example, for the model that will be used in the illustrative example 16 θ s corresponds to the Arias intensity. This means that Þ with x corresponding to the remaining model parameters excluding θ s , s i (x, w) representing the output Y g i θ; w ð Þ for θ s = 1 and f(θ s ) being a simple function of θ s . Without loss of generality, we will adopt here this assumption, ie, representation The choice of square root for function f(.) is made here simply to match the model used in the illustrative example. Numerical approach discussed herein can accommodate any other function. This setting leads to modification of the metamodel input to x, and, respectively, of metamodel output to the statistics of s i (x, w) conditional on x. Once the latter statistics are known the statistics for Y g i θ; w ð Þconditional on θ can be easily obtained as will be shown in Section 4, because the relationship to the remaining component of vector θ, θ s , is explicitly known. This modification ultimately reduces the dimension of the metamodel input, which as discussed earlier facilitates improved accuracy. Note that similar to Y g i , s i is also a random variable with randomness in its description stemming from w and x.
The approximated statistics, the ones required for calculating F p1 , are the conditional on x mean and variance of ln(s i ), given by where independence of w from x was used in deriving the last equalities in both expressions. These statistics define the metamodel output. It should be pointed out that this output is independent of the predictive relationships; rather is simply functions of x (reason for introducing the functional dependence on x notation). This is what contributes to the efficiency of the approach: the surrogate model is established with respect to the low dimensional vector x and can be then leveraged to evaluate the required statistics for different selections of the predictive models (ie, different μ and Σ) as detailed in the next section. For developing the metamodel, a database with n observations is initially obtained that provides information for the 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 to avoid extrapolations, ie, domain covered by p(θ|μ,Σ), approximated in this case as domain for p(θ|μ r (z),Σ r ) (for further details for definition of X please refer to Tsioulou et al 28 ). The predictions provided through the ground motion model for each x j are then established through the following process: (Step 1) generate n w sample acceleration time-histories for different white-noise sequences {w k , k = 1, …, n w } using θ s =1 for all samples; (Step 2) for each sample evaluate the response outputs of interest {s i (x j , w k ); i = 1, …, n y } using responsehistory analysis for spectral IMs; (Step 3) estimate the statistics (logarithmic mean and variance) over the sample-set to obtain ln s i x j À Á À Á and σ s i x j ð Þ. Using this database, the surrogate model can be formulated following the approach presented in Tsioulou et al. 28 Only difference is that output for metamodel also includes σ s i x ð Þ. The computationally intensive aspect of the formulation is the development of the database which requires response-history analysis. This needs to be performed, though, only once. As soon as the metamodel is established using this database, it can predict ln s i x ð Þ ð Þ and σ s i x ð Þ for any other x desired. Metamodel predictions can be also vectorized, 34 something that will be leveraged in the numerical FIGURE 1 Overview schematic of the hazard compatible stochastic ground motion modeling approach optimization discussed in the next section. The accuracy of the metamodel depends on the number of experiments n as well as the exact selection of these experiments. Details for both these tasks, including a sequential, adaptive metamodel formulation that gradually increases n, stopping when pre-specified accuracy criteria are satisfied, are provided in Tsioulou et al. 28

| Calculation of statistics of interest
Evaluation of objective F p1 boils down to estimation of Equations 2 and 3. Using representation where Cov[a,b] stands for the covariance between random variables a and b. The statistics with respect to θ s , that is the mean E[ln(θ s )] and variance Var[ln(θ s )] can be readily calculated using the marginal distribution p(θ s |μ,Σ). The statistics that involve s i may be calculated using the metamodel to approximate variability with respect to w. Using the laws of total expectation and variance we have for these statistics where the expectation (integrals) with respect to x or θ have been explicitly expressed in all these equations. These integrals address the variability stemming from the predictive relationships and need to be calculated through Monte Carlo simulation (MCS). Rather than performing two different MCS: one for the integrals involved in the expectation E[ln(s i )] and variance Var[ln(s i )], which require samples from the marginal distribution p(x|μ,Σ), and one for the covariance Cov[ln(θ s ), ln(s i )], which requires samples from joint distribution p(θ|μ,Σ), a single MCS is performed, utilizing a common set of samples for all these statistics. This leads to the following approximations for the quantities needed for Equations 14 and 15 where θ j s ; x j Â Ã correspond to samples from p(θ|μ,Σ), N s is the total number of samples used, and ln s i x ð Þ ð Þand σ s i x ð Þ are approximated through the kriging metamodel for each one of these samples. Utilizing vectorized manipulations for the metamodel predictions, both these quantities can be calculated with very small computational effort, meaning that the MCS-based estimation of Equations 19 to 21 can be performed very efficiently. Further numerical details (computational times) are provided in the illustrative example.
For reducing the relative importance of the MCS estimation error within the multi-objective problem of Equation 4, common random numbers (CRN) are utilized. 35 This is facilitated by getting first samples in the standard Gaussian space (which is independent of μ and Σ) and then transforming them to the desired samples for θ. Approach is equivalent to transforming integrals in Equations 16 to 18 to standard Gaussian space and for record-based model is easily performed, as these models typically entail 16 a Gaussian predictive model for v. The CRN is implemented by using the same sample set for the standard Gaussian samples across the entire optimization as also shown earlier in Figure 1. Adoption of CRN creates a consistent estimation error in the MCS application for the different examined μ and Σ values and therefore improves the optimization accuracy for identifying the correct Pareto front. 35 In other words, it allows use of smaller value for N s because it reduces the relative importance of the MCS estimation error within the optimization.

| Multi objective optimization
Calculation of statistics given by Equations 14 and 15, utilizing MCS estimates of Equations 19 to 21, facilitates an efficient approximation for performance objective F p1 given by Equation 8. If the lognormal assumption for the distribution of Y g i is not used, then F p1 can be estimated through the approach discussed in Appendix A, leveraging again MCS principles. Note that calculation of objective F p2 is computationally trivial. Therefore, through the introduction of the metamodel, an efficient estimation of both objectives involved in the optimization described by Equation 4 can be established. It should be noted that in Tsioulou et al 28 the explicit incorporation in the optimization of the metamodel error was also considered, established through appropriate modification of the equivalent objective function F p1 . It was shown, however, that this reduces computational efficiency and is not necessary if the underlying surrogate model has high accuracy. As calculation of objective F p1 requires N s evaluations of the metamodel in a MCS setting (only 1 evaluation was needed in Tsioulou et al 28 ) the inclusion of the metamodel error is not advocated here as the computational burden is expected to be higher. Additionally, the benefits of explicitly including this error in the current formulation are expected to be smaller because objective F p1 represents a statistical quantity over θ. In evaluating the necessary statistic, the response, and therefore associated metamodel errors, over different θ values are averaged. Potential large errors that may exist for specific θ values end up averaged with smaller errors from other θ's and therefore do not impact F p1 estimates as much as they would if these θ values were the only contribution to F p1 (as was the case in Tsioulou et al 28 ).
For solving the resultant multi-objective optimization, a variety of numerical approaches can be utilized. 36 In Tsioulou et al, 28 two such approaches were examined, one gradient-free and one gradient-based, and preference was ultimately given to the former because it can support higher robustness and computational efficiency for identifying the entire Pareto front. This recommendation is adopted here, and a gradient-free, random search approach is implemented as follows. A large number of n bc candidate solutions for μ and Σ are generated that are close to μ r (z) and Σ r . This is established by creating uniform random samples centered around μ r (z) and Σ r . The range for these samples is chosen so that the value of F p2 does not become excessively large, because the latter indicates significant departure of the modified model from the initial one, which might produce ground motions with unrealistic characteristics. Because estimation of F p2 is simple, a pre-screening can be implemented, rejecting any candidate solution with (large) value of F p2 over some desired threshold. Value close to 10 is appropriate for the latter threshold. This choice corresponds to modification of the standard deviation up to 50% and modification of μ within hypersphere 5 standard deviations away from μ r (z). Once all n bc candidate solutions are obtained, objective functions F p1 and F p2 are calculated for each of them. Estimation of objective F p1 in this case leverages the computational efficiency of the metamodel in performing vectorized predictions: the calculations are simultaneously performed for all n bc candidate solutions, or using subsets with a lower number of members depending on the available computational resources (memory can be a problem for vectorized operations depending on n and n bc 28 ). The dominant solutions representing the Pareto front can be then readily identified by comparing the values for the two objectives over all candidate solutions. The only challenge is that the value of n bc needs to be large in order to obtain an adequate representation of the Pareto front. Proper selection for n bs is examined in the illustrative example. Once the Pareto front has been identified, a solution (among the Pareto set) can be adopted using any desired criterion. 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 16 with model parameter vector including the arias intensity I a , the significant duration D 5-95 , the time corresponding to 50% of the intensity t mid and the associated spectral frequency ω mid , the rate of change for that frequency ω ' , and the damping ratio ζ f for the excitation spectrum. The arias intensity simply scales the ground motion, so it corresponds to parameter θ s in representation Y g i θ; w ð Þ¼ ffiffiffiffi θ s p s i x; w ð Þ, with the 5 remaining parameters corresponding to x. For the target hazard, GMPEs used in the Western US are selected here, [37][38][39][40] whereas the suggestions in Kaklamanos et al 41 are adopted to estimate unknown inputs for some of the GMPEs. As target, IM predictions of the logarithmic mean and variance from individual GMPEs as well as the average of their predictions are used. Note that the latter provides single target hazard for each structural period examined, simply that hazard is obtained by averaging information from multiple GMPEs. All computations are performed in a quad-core 3.0-GHz Xeon processor with 16 Gb of RAM. Fault and site characteristics are taken, respectively, as strike-slip fault and shear wave velocity in upper 30 m of soil V s,30 = 800 m/s. For moment magnitude M and rupture distance R, different values will be examined. As IM for the seismic hazard description, the peak pseudo-acceleration (PSA) (S a ) for a single degree of freedom system with 5% damping ratio is utilized. Different ranges of structural periods will be examined for S a and, unless otherwise specified, the objective function F p1 is estimated as the weighted average of the entropies for each scalar IM as in Equation 7, with the weights chosen as γ i = 1, so that no specific structural period is prioritized.

| Details for metamodel development
For the characterization of domain X and the selection of the support points for the metamodel development, the same approach as in Tsioulou et al 28 is adopted, the only difference being that a larger domain is considered here with a relative increase of 70% compared with Tsioulou et al. 28 The reason for the latter is to support higher accuracy in the MCS implementation: in Tsioulou et al 28 evaluation only at the new (modified) predictive relationships was required so proximity to the initial predictive relationships needed to account only for that modification, whereas here the MCS estimation will need to utilize the metamodel predictions for parameters even further away from the modified (and therefore from the initial) mean predictive relationships.
Logarithmic mean and variance for S a for different periods, the ones used by the aforementioned GMPEs, is adopted as the response output for the metamodel development whereas the white noise samples are chosen as n w = 100. Three different accuracy criteria with associated coefficient of determination 0.92, 0.94, and 0.96 are selected for the adaptive metamodel development. 28 This leads to metamodels with 1500, 3000, and 4500 number of support points, respectively. Note that the accuracy of the established metamodels is much higher for prediction of the logarithmic mean rather than the logarithmic variance. For example, for the metamodel with 4500 support points the average coefficient of determination is 0.99 for the logarithmic mean and 0.94 for the logarithmic variance. This trend agrees with past studies 42 that have examined the use of metamodels in approximating seismic hazard when stochastic ground motion models are utilized for the description of the latter. This should be also viewed as a positive attribute of the metamodel approximation for the purposes of this study: the logarithmic mean has a higher importance towards the hazard description for the typical earthquake engineering applications that have logarithmic variance smaller than 1 (the scaling of difference of logarithmic means by 1/ σ 2 i z ð Þ >1 in Equation 8 when estimating discrepancy form target also reveals this), which means that a higher degree of confidence exists for properly approximating this hazard than the average coefficient of determination (averaged over both logarithmic mean and variance) indicates, because the accuracy for the more important component (logarithmic mean) is higher.
Estimation of metamodel response for 10 000 samples requires 3.6, 7.5, and 10.7 seconds for the metamodels with 1500, 3000, and 4500 support points, respectively. Note that adoption of larger value of samples prohibits efficient vectorization of operations for the n = 4500 points due to memory restrictions. Calculation of objective function F p1 for N s = 100 MCS samples requires 0.005, 0.10, and 0.16 seconds for the metamodels with 1500, 3000, and 4500 support points, respectively. It should be stressed that, as advocated earlier, evaluation of the metamodel across multiple candidate solutions is simultaneously performed to better leverage the computational efficiency allotted by the vectorized metamodel evaluation. For example, for N s = 100, 100 different candidate solutions are examined; this leads to a total of 100 × 100 = 10 000 samples for the metamodel evaluation, which avoids memory problems even for the metamodel with 4500 support points.

| Validation of lognormal distribution assumption
First the validation of the lognormal assumption for the response output Y g i distribution is examined. This is performed separately for the response output from the exact stochastic ground motion model as well as for the metamodel approximation. Distribution p g (ln(Y i )| μ, Σ) is approximated though KDE (Equation A1) as discussed in Appendix A using n s = 1000 samples of Y g i . Comparison is primarily expressed though objective F p1 because this is the critical quantity utilized in the proposed framework. The target hazard corresponds to the average of the aforementioned four GMPEs. The predictive model proposed in Rezaeian and Der Kiureghian, 16 referenced herein as unmodified model, is used for defining μ and Σ, but comparison extends over a wide range of seismicity scenarios, M in range [6 8] and R in range [10 100] km. The latter guarantees that validation extends over the wider range that the optimization is going to examine with respect to the predictive relationships. Figure 2 shows F p1 evaluated either through the lognormal assumption or through the KDE approximation for three periods T s = [0.01 0.5 2] seconds (columns of the figure). In this instance F p1 is evaluated for each period separately to assess impact over specific structural characteristics, rather than averaged over the different IMs. Top row corresponds to comparison with respect to the actual stochastic ground motion model and bottom row to comparison with respect to the metamodel approximation. The comparisons in this figure show exceptionally close agreement for objective F p1 between the lognormal assumption and the KDE approximation (compare the two curves in each subplot) over the entire seismicity range and for all structural periods. This agreement holds for both the actual model and the metamodel. The results in the figure allow for some additional observations. First there is a very good agreement between the metamodel and the actual model (compare the results across the rows). This provides a preliminary validation, examined further in the following sections, of the proposed, metamodel-based framework. Secondarily, for certain seismicity ranges, the unmodified model has large discrepancy (large F p1 values) from the target hazard. This validates the claim that motivated this study, that existing approaches for selecting predictive relationships do not necessarily provide a close match to the desired hazard for some structures or seismicity scenarios. For The compatibility with respect to the lognormal assumption is further examined in Figure 3, which shows the cumulative distribution functions (CDF) for Y g i for a particular seismicity scenario, corresponding to M = 6 and R = 20 km. Similarly to Figure 2, each column corresponds to a different period T s . The results in Figure 3 show, again, a very good match not just between the CDFs using the KDE approximation or the lognormal assumption but also between the CDFs using the actual model and the metamodel. Overall, the comparisons in Figures 2 and 3 justify the use of a lognormal assumption for the distribution of Y g i . Additionally, they provide a first validation for the accuracy of the metamodel approximation, when compared with the actual ground motion model, with respect to both the entropy ( Figure 2) but also the exact distribution of Y g i (Figure 3).

| Optimization details and metamodel accuracy
This section examines details of the numerical solution of the optimization problem. For the MCS of Equations 19 to 21, Latin hypercube sampling was adopted, and for the sample size N s different values in range [20 150] were examined. The selection of the exact value of N s is a compromise between numerical efficiency (it proportionally impacts the computational burden) and robustness of the multi-objective optimization (identification of correct Pareto front). After sensitivity analysis, a value of N s = 70 was chosen. The coefficient of variation for F p1 using MCS estimates of Equations 19 to 21 with N s = 70 was found to be in range 3% to 4% which should be deemed sufficient considering the fact that CRN are further utilized 35 to reduce the importance of the associated estimation error in the comparisons established across the optimization. The discussion moves next to the impact of the number of points used in the random search n bc . The considerations are similar to the selection of N s (efficiency versus robustness) although in this case the selection is not as straightforward because there are no simple statistics like the coefficient of variation to compare. The sensitivity analysis is performed instead by solving the optimization problem for different values of n bc . The target used in this section corresponds to structural periods T s = [0.4 0.5 0.75 1.0 1.5 2.0] seconds and hazard described by the average of the considered GMPEs. Three different seismicity scenarios M = 6-R = 20 km, M = 7-R = 30 km, and M = 8-R = 50 km are considered. The first scenario represents a case that the unmodified model provides a poor match to the target hazard, the second one achieves a good match, and the last case lies in between the other two. Figure 4 presents the Pareto fronts identified through random search for different number of points n bc. The metamodel with 4500 support points is used in this case. Only a few representative solutions and not the entire front are shown for clarity. Results are reported in this figure and in the remaining of the manuscript with respect to ffiffiffiffiffiffiffi F p1 p and ffiffiffiffiffiffiffi F p2 p to facilitate an easier comparison (differences of extreme values easier to discern). Figure 4 shows that the differences between the identified fronts for different n bc values are only minor, and they occur primarily for small ffiffiffiffiffiffiffi F p1 p large ffiffiffiffiffiffiffi F p2 p combinations. Note that the objective function ffiffiffiffiffiffiffi F p1 p values for the M = 7-R = 30 km scenario are very low (due to an already good match of the unmodified model to the target hazard); therefore, any identified differences between the Pareto fronts for different number of points n bc are of smaller relevance. This comparison shows that a value of n bc around 100 000 to 250 000 should be considered as sufficient for identifying an adequate representation of the Pareto front. An n bc value equal to 150 000 is adopted for the random search results presented in the remaining of the manuscript.
The impact of the metamodel accuracy is examined next. This is established by considering additionally the results obtained by using the exact stochastic ground motion model, which represents the measure for evaluating the actual hazard compatibility of the identified ground motion model. Comparison is performed across different Pareto optimal solutions. The Pareto fronts identified by using the metamodels with the three different number of support points are presented in Figure 5 for seismicity scenario M = 6-R = 20 km and M = 8-R = 50 km. Note that for seismicity scenario M = 7-R = 30 km which was also presented in Figure 4, the results are of limited interest because the unmodified ground motion model provides a good compatibility to the target hazard. This is the reason that this seismicity scenario is not presented here. In all cases, the random search is implemented with the same candidate solutions, to facilitate a consistency in the corresponding comparisons. Then, Figure 6 presents spectral plot comparisons for the solution (among the Pareto set identified in each case) corresponding to the minimum of F p1 for the seismicity scenario M = 6-R = 20 km. Similar trends hold for the M = 8-R = 50 km seismicity scenario (not reported here due to space constraints). In all plots, the predictions using the metamodel and the actual ground motion model are reported. Figure 6 offers comparisons in context of both average response (top row) using curves corresponding to median and 14 th to 86 th percentiles (denoted as median ± σ log herein), and logarithmic standard deviation (bottom row). The former assesses the hazard compatibility with respect to different IM statistics and the latter explicitly with respect to the IM dispersion.
The results show that good agreement is established between the metamodel predictions and the actual model predictions along the Pareto front for all metamodel cases: 1500, 3000, and 4500 support points. This is evident in both the Pareto fronts ( Figure 5) and the corresponding spectral plots (Figure 6), the latter indicating a good match in terms of both the mean and variability of the hazard. For the lower accuracy metamodel utilizing only 1500 support points greater differences exist, especially for smaller F p1 values, but the discrepancies even for it are overall quite small, significantly smaller than the ones reported in Tsioulou et al 28 where this metamodel was shown to lead to erroneous results, identifying suboptimal solutions when the metamodel error was not explicitly considered in the optimization. This should be attributed to the fact, discussed also in Section 4.2, that objective F p1 represents a statistical quantity over θ, and averaging over different θ values for evaluating F p1 reduces the potential influence of larger errors for specific θ values. This discussion and the good agreement reported in Figures 5 and 6 also further validate the choice to avoid the explicit consideration of the metamodel prediction error in the problem formulation (see discussion in Section 4.2). This consideration would increase the computational burden while providing negligible benefits in terms of the quality of the identified Pareto fronts. Overall, all examined here metamodels provide adequate accuracy in the identification of the Pareto front, with no need to explicitly consider the metamodel prediction error, with preference towards the metamodel with 3000 or 4500 support points. The latter will be utilized in all remaining comparisons in this manuscript. Some final remarks are warranted with respect to the overall computational cost. The primary computational burden of optimization of Equation 4 stems from the N s ·n bc evaluations of the metamodel required for the MCS of Equations 19 to 21 across the search for the Pareto front. For the scenario advocated earlier with N s = 70, n bc = 150 000 and use of metamodel with 4500 support points the time is 180 minutes per seismicity scenario. If the metamodel with 3000 support points is used instead the required time reduces to 120 minutes. Compared with the study Tsioulou et al, 28 this represents an important increase of the computational cost: there is ultimately a N s -fold increase of this cost for same n bc value, stemming from the MCS step.

| Implementation for different seismicity scenarios
With the computational details ironed out, the discussion moves next to the hazard compatibility established by the proposed modification of the ground motion model. Figure 7 shows results for 6 seismicity scenarios targeting seismic hazard given by the average of the aforementioned GMPEs for two different ranges of T s : T s = [0.4 0.5 0.75 1.0 1.5 2.0] seconds and T s = [0.4 0.5 0.75] seconds. These two different cases are referenced herein as long and short, respectively, period ranges. The proposed framework identifies 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 hazard. Choosing a shorter period range for the target IMs facilitates an overall better FIGURE 6 Spectral plots for the solutions corresponding to minimum of F p1 in the Pareto fronts identified in Figure 5 for the M = 6-R = 20 km seismicity scenario. Top row shows curves corresponding to median and median ± σ log for the response. Bottom row shows logarithmic standard deviation (σ log ) match (smaller F p1 values); this is anticipated because objective F p1 imposes less strict requirements (fewer number of components to match) in this case.
To select a solution within the identified Pareto set, one of the most common approaches within the multi-objective optimization literature 43 is to choose the one that has the smallest normalized distance from the utopia point, defined as the point in the Pareto front that corresponds to the idealized (unachievable) minimum of the two objectives across the front. Following the guidelines in Tsioulou et al 28 for choosing the normalization, the following distance metric is chosen The corresponding point is identified in all cases in Figure 7. Another option would have been to choose the solution that satisfies a pre-determined threshold for the match to the targeted IMs. These selections are finally demonstrated for a wide range of seismicity scenarios, M in range [6 8] and R in range [10 100] km, in Figures 8 to 13. For each scenario, three different Pareto points are selected, the one with smallest distance D p (θ) from the utopia point and the ones with objective ffiffiffiffiffiffiffi F p1 p smaller than 0.15 or 0.075. 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 8 shows plots for (i-first row) ffiffiffiffiffiffiffi F p1 p for U t and U n ( ffiffiffiffiffiffiffi F p1 p is constrained for the other two cases) and for (ii-second and third rows) ffiffiffiffiffiffiffi F p2 p for U t , C l , and C s ( ffiffiffiffiffiffiffi F p2 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 hazard given by the average of the aforementioned four GMPEs for both the (A) long and (B) short period ranges as well as (C) target hazard given only by GMPE 38 for the long period range. These scenarios are denoted herein as SC 1 , SC 2 , and SC 3 , respectively. Figures 9 to 11 show spectral plots for a 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 spectral curves corresponding to the target hazard, the unmodified model and the predictions by the three aforementioned model modifications are shown to facilitate comparisons. More specifically, Figure 9 shows curves corresponding to different statistics of the response FIGURE 7 Pareto fronts for different seismicity scenarios considering match to long (black) or short (gray) period range IMs for defining the seismic hazard. In each plot, the Pareto point with minimum distance from utopia point is shown with x and the utopia point with a □ (median and median ± σ log ) for the unmodified model and the model corresponding to the Pareto point with smallest distance from the utopia point (U t case). Figure 10 presents the same curves for the models corresponding to the Pareto points with average relative entropy thresholds C l and C s . Figure 11 presents same comparison directly in terms of logarithmic standard deviation (ie, IM dispersion). In all these figures, the curves corresponding to the target are also shown.
Finally, Figures 12 and 13 show the predictive relationships for all examined seismicity scenarios for the modified model corresponding to the Pareto point with smallest distance from the utopia point (U t case) for SC 1 . In particular, Figure 12 shows the mean model parameters θ of the modified predictive model (corresponding to μ for the aforementioned Pareto point). The unmodified ground motion model parameters (corresponding to μ r ) as well as another case that will be discussed in the next section are also included. Then, Figure 13 shows the ratio of standard deviation for the modified and unmodified model (comparison of square root of the diagonal points of Σ for the aforementioned Pareto point and Σ r ). Note that some of the curves shown in these figures have non-smooth characteristics. As discussed in Tsioulou et al, 28 this should be attributed to the fact that a discrete representation of the Pareto front is obtained and to the random search characteristics of the adopted optimization algorithm for identifying the Pareto front.
The baseline trends observed in the figures are similar to the ones in Tsioulou et al. 28 These trends are enhanced in this study with additional considerations with respect to the variability associated with the seismic hazard. The unmodified model does not provide a good match to the target hazard for the entire seismicity range as observed in first row of Figure 8 and also in the spectral plots in Figures 9 and 11. This is particularly true for the logarithmic mean, ie, median (Figure 9), with logarithmic variance, ie, dispersion (Figure 11), showing better compatibility to start with. The proposed modification (cases U t , C l , and C s ) significantly improve the match to the hazard (Figure 8), establishing a balance between F p1 and F p2 , with the characteristics of the balance depending on the criteria for selection of the final model among the Pareto optimal solutions. Similar to Tsioulou et al 28 when the unmodified model has larger discrepancies from the target hazard, then the modifications lead to larger values for F p2 , but still successfully identify models, independent of the implementation case, that provide an improved match to the IM target mean and variance. This is clearly observed for C l and C s cases; C s imposes a smaller discrepancy between the modified model and the target hazard and the modification leads to identification of a model with bigger differences from the original one (larger values for F p2 ) . The U t modification identifies a model with moderate discrepancy from the unmodified one, corresponding to values of ffiffiffiffiffiffiffi F p2 p in the range of 0.2 to 0.8, 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 hazard, and therefore modification of it provides limited benefits. This is perhaps better captured by the C l case, which corresponds to low ffiffiffiffiffiffiffi F p2 p values for this seismicity range, significantly lower than the U t modification. This indicates that a satisfactory match to the target hazard, ie, a match satisfying the predefined compatibility threshold, is established with no need to greatly modify the initial predictive models. Thus, selection of the Pareto optimal model based on a targeted accuracy to the GMPEs, ie, value for F p1 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. Selection of a small threshold (C s case), however, results in identification of a model with large discrepancies from the unmodified model (large F p2 values). This model will typically be far away from the U t case and belong in a steep part of the Pareto front (check Figure 7 earlier), meaning that small improvements in F p1 come at a large increase of F p2 . A multi-level selection criterion similar to the one proposed in Tsioulou et al 28 may be therefore advocated: select the Pareto optimal solution that satisfies a certain accuracy threshold for F p1 unless this solution leads to a greater F p2 value than the Pareto optimal solution with minimum distance

FIGURE 9
Spectral plots of average response for seismicity scenarios (different subplots) corresponding to combinations of M [6. 2, 6.8, 7.4, 8] and R [30,60,90] km, for the target hazard (target), the unmodified ground motion model (U n ), and the modified ground motion model corresponding to the minimum distance from utopia point (U t ). Implementation scenario shown corresponds to matching to the average considered GMPEs and long period range. Curves corresponding to median and median ± σ log shown from the utopia point. If the latter happens, then select the Pareto optimal solution with minimum distance from the utopia point.
From the comparison of the different implementation cases shown in Figure 8, it is evident that selection of a shorter period range [compare cases (B) and (A)] establishes an easier match to the target, as it results in smaller overall values for F p1 and F p2 . The selection of a different target [compare cases (A) and (C)], although not imposing any additional constraint in the optimization implementation, leads to different results. Overall, the discrepancy between the unmodified model and the target hazard follows the same trends in all implementation cases. This demonstrates the versatility of the proposed framework as the modification facilitates an enhanced hazard compatibility for all IMs and hazard targets considered.
The spectral plots in Figure 9 to 11 provide the decomposition of the overall F p1 match to the different structural periods and the statistics of the IM distribution (median and dispersion). The comparison between the unmodified and modified cases shows that the match to both the target median and dispersion are separately addressed within the proposed optimization. This is also supported by the trends in Figures 12 and 13; not only the mean ( Figure 12) but also the standard deviation ( Figure 13) of the predictive model is adjusted. Overall for the ground motion model examined here, the match to the target dispersion offered by the unmodified model, and subsequently by the modified one, is quite good (compare relative discrepancies for U n in Figures 8 and 11). Although the proposed modifications also impact this dispersion (Figure 11), the impact on the median (Figures 9 and 10) is more substantial. This should be attributed to ability to influence value of objective F p1 more by adjustments in ln Y g i μ; Σ ð Þ À Á rather than σ g i μ; Σ ð Þ [formulation of Equation 8 and scaling of differences from the target for the former by σ 2 i z ð Þ also reveals that], something automatically leveraged by the optimization. With respect to the different modifications, the C s always provides a better match compared with C l one, although the spectral curves of the two modifications are very close to each other especially for lower magnitude ranges M = 6.2-6.8. The U t modification also provides a good

Cs Cl
Target Cs Cl  FIGURE 10 Spectral plots of average response for same seismicity scenarios as in Figure 9 (different subplots) for the target hazard (target), and the modified ground motion model corresponding to value ffiffiffiffiffiffiffi F p1 p smaller than 0.15 (C l ) or 0.075 (C s ). Implementation scenario shown corresponds to matching to the average considered GMPEs and long period range. Curves corresponding to median and median ± σ log shown match to the target, frequently very close to the C l and C s cases, depending on the original match of the unmodified model (ie, how much is the modification truly needed). This observation further supports the multi-level selection criterion discussed earlier.
The results in Figures 12 and 13 show that the model modification leads to similar characteristics as observed for the unmodified model. This guarantees that the proposed ground motion model modification does not deviate significantly from observed regional trends and was achieved by incorporating this deviation as an objective in the problem formulation (objective F p2 ). Parameters I a , ζ f , and ω ' show bigger variability compared with their initial mean values ( Figure 12). This should be attributed to a greater sensitivity with respect to them of the resultant ground motions and agrees with the trends reported in Tsiouou et al. 28 With respect to the adjustment of the variability of the predictive model, the ratio of standard deviations remains close to 1 ( Figure 13) with values in the range of 0.8 to 1.05, indicating small (but not negligible) overall adjustment.
The overall discussion shows the importance of the established framework: once the initial metamodel is developed, it can support the efficient identification of ground motion models that (1) match conditional hazard for any desired IMs and chosen period range while (2) maintaining a small deviation from the initial predictive models. This can be seamlessly repeated for any seismicity scenario. The final ground motion model modification can be chosen based on the criteria discussed earlier.

| Comparison to modification of mean value characteristics only
As discussed in Section 5.3, the computational cost of the proposed modification to match the probabilistic seismic hazard is significantly higher than previous efforts to match only the IM predictions corresponding to the mean predictive relationships for the ground motion model, 28 ie, completely ignoring the variability stemming from Σ in the predictive model for θ (using Σ = 0). This increase in computational burden stems ultimately from the need to estimate the statistics of the response when the variability of the predictive model p(θ|μ,Σ) is considered. It is therefore of interest to examine whether the computationally less demanding problem of matching only the IM predictions corresponding to the mean of the predictive relationships 28 can be adopted as a surrogate for the problem of interest here. This is    The quality of the solution obtained from this approximate problem may be then assessed by evaluating objectives F p1 and F p2 assuming distribution p(θ|μ m ,Σ r ), ie, adopting the initial variability of the predictive models Σ r , or p(θ|μ m ,0), ie, ignoring any variability in the predictive models and utilizing only the variability stemming from the white noise to calculate response statistics. The comparison can be performed with respect to the entire Pareto front for limited number of seismicity scenarios or with respect to a specific solution along the front over a wider range of scenarios. Due to space constraints, the latter is only reported here, with specific solution chosen as the point with minimum distance from the Utopia point. The solution available from Tsioulou et al 28 is directly utilized for μ m . The case corresponding to p(θ|μ m ,Σ r ) is denoted as U tm and case corresponding to p(θ|μ m ,0) as U tmn . Results are shown in Figures 14 and 15, following same guidelines as the study reported in Figures 8 and 9. Comparison with respect to the dispersion only (similar to Figure 10) is not reported due to space constraints. Figure 14 compares U n , U t , and U tm across both objectives over a range of seismicity scenarios (ie, adds U tm curve in the results reported in Figure 8), and Figure 15 presents spectral plots for the target seismic hazard, U tm and U tmn . Latter figure should be compared directly to Figure 9 to evaluate the relative advantages of U t (adding that curve in this plot is avoided to improve clarity of the presentation). Solution for the physical parameters corresponding to μ m has been also reported earlier in Figure 12. Note that this solution is same for both U tm and U tmn .
With respect to U tmn first (Figure 15), the results show the importance of considering the variability of the predictive relationships. Reliance only on the variability stemming from the white noise provides significantly lower variability for the seismic hazard than the variability prescribed by GMPEs (compare median and median ± σ log curves for target and U tm cases). This ultimately leads to large values for F p1 for U tmn which is the reason that results for it are not reported in Figure 14. Comparison now between U t and U tm in Figure 14 shows that the explicit optimization for the probabilistic hazard provides for some scenarios a noticeably better match (smaller F p1 values for U t ) for the same level of modification of the initial probability model (similar F p2 values for U t and U tm ). U tm modifications even underperform the unmodified model U n for some scenarios, corresponding to cases for which U n provides an adequate match to the target hazard to start with. Same trends are observed in the spectral plots in Figure 15. Overall U tm is shown to provide an improvement over U n when the initial match to the hazard is not adequate, although it might underperform U t due to its inability to explicitly accommodate the hazard variability. It should be also pointed out that based on the results in Figure 12, even the mean value vectors for the predictive relationships identified by different modification approaches (U t and U tm ) are different. A final interesting comparison can be established with respect to the mean FIGURE 14 Results for ffiffiffiffiffiffiffi F p1 p and ffiffiffiffiffiffiffi F p2 p for unmodified ground motion model (U n ) and modified ground motion model with minimum distance from utopia point for matching the complete probabilistic hazard (U t ) or the hazard corresponding to mean predictive relationships (U tm ). For the latter estimation of objectives adopts variability of the initial (U n ) predictive model but with the corresponding updated predictive mean. 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 38 [Colour figure can be viewed at wileyonlinelibrary.com] spectral response for U tm and U tmn (black curves in Figure 15). Consideration of the variability in the predictive model after the optimization (U tm case) impacts even the mean of the predicted response, not only the variability of that response, and overall moves that response further away from the desired target, which, recall, was the objective in the underlying modification for the mean predictive relationships μ m . This discussion shows that even though the surrogate optimization for only μ m may provide an improvement over the unmodified model for scenarios of high initial discrepancy from the target probabilistic hazard, the post-consideration of variability in the predictive models is problematic. Setting initially Σ = 0 to identify μ m and then calculating the hazard for predictive model (μ m ,Σ r ) provides ultimately a lower quality fit to the target hazard. Despite the higher computational burden associated with it, the simultaneous modification of the entire predictive model for θ (both μ and Σ) is therefore advocated.

| CONCLUSIONS
The modification of stochastic ground motion models to establish hazard compatibility for specific seismicity scenarios was discussed in this paper. The hazard for each scenario was described with respect to some IM of interest and a probabilistic description was adopted for it, for example defined through mean and dispersion characteristics. The modification of the ground motion model was defined as an adjustment of the probabilistic predictive models/relationships that relate the parameters of the ground motion model to seismicity characteristics. Both the mean of the predictive model and the associated variance were adjusted. The proposed modification was defined as a bi-objective optimization with dual objective of minimizing the discrepancy between the hazard for a given structure/site and the predictions established through the stochastic ground motion model, while maintaining a small deviation from the original predictive relationships, assumed to facilitate similarity to observed regional trends. This setting extends previous efforts of the FIGURE 15 Spectral plots for seismicity scenarios (different subplots) corresponding to combinations of M [6. 2, 6.8, 7.4, 8] and R [30, 60, 90] km, for the target hazard (target), and modified ground motion model with minimum distance from utopia point for matching the hazard corresponding to mean predictive relationships. For the latter, the cases adopting variability of the initial predictive model (U tm ) or no variability (U tmn ) are presented. Implementation scenario shown corresponds to matching to the average considered GMPEs and long period range. Curves corresponding to median and median ± σ log shown authors 28 that examined modification of only the mean predictive relationships with goal to match the corresponding IM predictions, ignoring any variability in either of these two components. The relative entropy was adopted as metric to quantify the objectives considered in the multi-objective optimization problem whereas the same surrogate modeling framework as in Tsioulou et al 28 was utilized for an efficient optimization. Emphasis was placed on the estimation of the various response statistics needed for the entropy calculation, and a MCS approach was advocated for it coupled with assumption for lognormal distribution of the response when considering the variability of the predictive models/relationships. The computational approach explicitly considered the fact that most ground motion models involve a separate parameter that impacts their scaling. This parameter was separately treated with respect to both the surrogate model development and the MCS. Different statistical assumption for the distribution of the ground motion model output were also examined for the evaluation of the entropy for the first objective.
In the illustrative example, the proposed framework was applied using a recently developed record-based stochastic ground motion model. It was shown that lognormal distribution assumption for calculating the first objective provides an adequate approximation for the application at hand, whereas the metamodel-aided optimization can facilitate an accurate identification of the Pareto front even when lower accuracy metamodels are utilized, a feature that did not hold in study. 28 The necessity to calculate response statistics through MCS increases, though, the computational burden of the approach. Application to wide range of seismicity scenarios and different approaches for determining the seismic hazard (different IMs or sources for the target values) were examined, illustrating the advantages for the proposed framework: it allows significant improvements to the target hazard match, especially for seismicity scenarios for which the unmodified model provides a poor initial fit, with minor only modifications to the original predictive model, something that can guarantee a good agreement with observed regional trends. With respect to selection of the final model across the identified Pareto front, same recommendation as in Tsioulou et al 28 is made: select the Pareto optimal solution that satisfies a certain accuracy threshold for match to the target hazard (F p1 constraint) unless this solution leads to a greater modification for the predictive model (F p2 value) than the Pareto optimal solution with minimum distance from the Utopia point. Finally, the approach presented in this paper was compared with the modification of only the mean predictive relationships to match the corresponding hazard for these mean predictions. 28 It was shown in this case that the latter approach may provide an adequate surrogate for seismicity scenarios with high initial discrepancy to the target hazard, though overall it is better to explicitly consider the impact of the variability in the predictive models (ie, calculate response statistics) and simultaneously modify the entire predictive model (ie, not only focus on the mean of the predictive relationships).
The main limitation of the approach is the significant computational burden for performing the multi-objective optimization to identify the Pareto front, a burden stemming from the MCS step. Considering the fact that the proposed modification needs to be repeated for each seismicity scenario of interest, further reduction of this burden, which will have to come from a more computational efficient implementation of the surrogate model predictions, is an important extension of this work, currently under investigation by the authors.