Engineering-oriented ground-motion model for Israel

This study presents a response-spectral Ground Motion Model (GMM) for Israel, called KYB22 herein and derived with practical applications in mind. This model is based on the former work by Maiti et al. (Bull Seismol Soc Am 111: 2177–2194, 2021), who derived a suite of nine Fourier amplitude spectra GMMs, using both empirical data and calibrated point-source simulations. In this study, a weighted average of the Maiti et al. (Bull Seismol Soc Am 111: 2177–2194, 2021) FAS models is computed and a synthetic database is created. Next, this database is converted to the Response spectral domain using the random vibration theory and a new GMM is regressed, constraining the magnitude and distance scaling on the synthetic response-spectral data. Site scaling is represented by VS30 and is a combination of empirical scaling with other considerations which are discussed in the text. Nonlinear site response is constrained from a global model, as are finite-fault effects, such as hanging-wall, mechanism and top of rupture—which cannot be constrained from the data because it does not contain enough large magnitude data. State-wide hazard is then computed using KYB22, comparing results with other GMM combinations. It is found that the hazard results obtained by using KYB22 as a backbone model are comparable to results obtained using other popular combinations of GMMs in the logic tree. Therefore, we recommend using the new GMM as one of the branches within the ground-motion logic tree when conducting seismic hazard calculations for Israel.


Introduction
Ground motion models (GMMs) are a critical component within seismic hazard assessment (SHA). They allow us to compute the expected engineering load at the site of interest, given information about the earthquake source, distance, and site conditions. In regions where seismic data is abundant, GMMs are typically fully empirical, using the available data to constrain the model behavior. In data-poor regions, empirical models cannot be 1 3 fully constrained and hence two main approaches have been traditionally adopted: (a) either a global model is used, with or without regional adjustment factors (e.g. Campbell 2003; Gülerce et al. 2016;Shahjouei and Pezeshk 2016), or (b) simulations are used to synthetically create the data on which the model is then regressed (e.g. Goulet et al. 2021;Yenier and Atkinson 2015).
In Israel, there are currently no models that were developed, or calibrated, to local data and conditions, which are also applicable for standard engineering practice. This is mostly due to lack of empirical data within the magnitude range applicable to hazard. The existing GMMs are either based on too little data (Meirova et al. 2008), do not cover the full spectrum (Lior and Ziv 2018), or do not include all relevant parameters for engineering applications . That is why the Israeli national seismic hazard map, used in the building code, uses the Campbell and Bozorgnia (2008) GMM, which is therefore often used in practice, without any adjustment or validation. Nevertheless, Avital et al. (2018) showed that the uncertainty associated with the GMM can significantly increase the hazard in Israel.
In a recent study, Maiti et al. (2021, MYK21 hereinafter) developed a suite of alternative GMMs for Israel in terms of the Fourier amplitude spectra (FAS). The models are based on different datasets and different functional forms, to ensure that both modeling and parametric epistemic uncertainty are captured. The datasets include an empirical dataset (including only ground motions from Israel, hence very limited in magnitude range and frequency band) and four alternative simulation datasets. The simulations used the pointsource Stochastic Method SIMulation (SMSIM) platform (Boore 2003), with four alternative sets of calibration parameters, all compatible with the empirical data. The models were then regressed using two host-models-the Bayless and Abrahamson (2019) model, which is based on the NGA-West2 database, and the Bora et al. (2015) model, which is based on the European RESORCE database. The regression was performed such that some terms were locked to their original host values and some terms were regressed to the local data-depending on data availability within the regressed data. Eventually, all nine models are functions of moment magnitude (Mw), rupture distance (Rrup) and the time-average shear-wave velocity of the top 30 m (V S30 ). Additional scaling terms, such as mechanism or hanging-wall are not included in the MYK21 models. A short summary of the MYK21 study, including a table with calibration parameters of the four simulation sets, is provided as supplementary material to this manuscript, KYB_SI3.
In this study, we continue the work of Maiti et al. (2021), by developing a practiceoriented response-spectral GMM and testing its effect on the computed seismic hazard. The model is based on a weighted average of the nine alternative MYK21 models, with additional parameters that are significant for forward-analysis engineering applications, as will be explained in the following sections.

Dataset
The dataset used for regression is a synthetic dataset, which has a uniform distribution in the magnitude-distance (M-R) domain (see Fig. 1a). The synthetic dataset was created as follows: first, a weighted mean model in FAS was calculated using the nine models in Maiti et al. (2021). A weighted model is used in lieu of using a compilation of all nine models, so that the separation between aleatory variability and epistemic uncertainty is maintained. In this way, the MYK21 model variations are mapped into the logic tree branches instead of into the variability of the dataset itself. The weights were distributed such that the empirical model received 50%, the low stress-drop simulation models received 25% and the high stress-drop simulation models received the remaining 25%. This was done to ensure equal weight between empirical and simulation-driven models. Next, the FAS synthetic database was converted to response-spectral domain (RSP), by using the open-source software PyRVT (see data and resources section). Using PyRVT requires specification of a peakfactor model, as well as a duration model. In a preliminary stage, alternative peak-factor, as well as duration models, were used, and the sensitivity of results to the choice of either was found to be insignificant. Hence, for the peak-factor model, we use the built-in BJ84 calculator, based on Boore and Joyner (1984). For the duration model, we use the duration model proposed by Lee and Green (2014) for Western North America (WNA). This model was found to be very similar to that proposed more recently by Ashfari and Stewart (2016), only smoother, such that the large-magnitude data did not oversaturate. Note that the synthetic database includes only one value of V S30 (600 m/sec), due to the original FAS simulation-based models being derived for a single value of V S30 .
The empirical GM database presented by Yagoda-Biran et al. (2021), is used herein to validate and calibrate the regressed model. The M-R distribution of our empirical database is presented in Fig. 1b. Because the empirical database is mostly populated by 3 < M < 4 events, finite-fault effects were not accounted for in MYK21 and hence are not represented in the synthetic database.

Regression
The functional form of the new GMM generally follows the Abrahamson et al. (2014) model, referred to as ASK14 hereafter, except for a characteristic break in the distance scaling. The general form of the median model is as follows: (1) in which f1 describes magnitude and distance scaling and is fully regressed on the synthetic database; f5 describes the site-response and is constrained partly by local data and partly globally; and f4, f6, f7 and f8 are finite-fault effects which are fully adopted from ASK14, as will be explained below.

Base form
The basic form of the magnitude and distance dependence, for strike-slip earthquakes, is similar to that of ASK14, but includes an additional break in the distance scaling between 70 and 200 km. As explained and demonstrated in Maiti et al. (2021), the empirical data shows a very characteristic break in the distance scaling, leading to diminished attenuation between 70 and 200 km, typically associated with reflections from the Moho discontinuity. Hence, this feature is also reflected by the associated simulation data and hence in the median FAS model suggested by Maiti et al. (2021).
Due to the additional break in the distance scaling, we formulate the basic form slightly differently than the way it is presented in ASK14, such that the magnitude scaling breaks are defined by three additional magnitude-dependent parameters-Ma, Mb, and Mc, which are required to account for magnitude-dependent distance scaling.
The functional form of the basic model is as follows: where and c 4M is identical to its form in ASK14: where c 4 is period-independent and is c 4 = 6. (2) Finally, the additional magnitude-scaling parameters are defined as follows: where M 1 = 6.75 and M 2 = 5. All model coefficients are provided in Table 1 in the supplementary material KYB_SI1.

Site response
Our site-response term is similar to that of ASK14, with some modifications-it includes a linear and a non-linear term and is based on V S30 as the main site characterization parameter. The linear term is constrained by data and simulations, as will be explained shortly. The nonlinear term, which was constrained in ASK14 by 1D simulations , is adopted in full. It uses the spectral acceleration on hard rock-Sa 1180 -to define the strength of shaking. However, we do not include the second limiting shear-wave velocity term, V1, because surface rocks in Israel are on average stiffer than in California (CA) and we do not see evidence of decreased scaling at higher V S30 in our database. The functional form of the site response term is as follows: in which b, n, c, and V lin are fixed to the ASK14 values (following Kamai et al. 2014), and a 10 is constrained externally, as discussed below, because the synthetic dataset does not include V S30 scaling (is computed for one value only). Table 1 List of model coefficients and method used to constrain them *a 2c constrained to be − 0.5, to allow for a smooth transition from body-wave to surface-wave geometrical spreading, following Chiou and Youngs (2014) **M 1 slightly modified from ASK14, such that it is fixed rather than period-dependent The linear V S30 scaling terms, as obtained from the two empirical models in Maiti et al. (2021), are presented in Fig. 2. The original host model terms are also shown, for comparison. Finally the V S30 scaling as obtained from a set of linear-elastic 1D simulations, using STRATA (Kottke and Rathje 2008), is also shown. In the 1D simulations, the soil velocity profiles are based on a velocity profile database of about 1800 profiles in Israel (Baram et al. 2020) and amplification is computed with respect to an outcropping generic rock profile, for which V S30 = 1100 m/sec (Baram et al. 2019).
The comparison shows that at high frequencies (approx. f > 2 Hz), all models are generally within a similar range. However, at low frequencies, the local amplifications-both analytical and empirical-predict much less amplification (smaller negative value) than the two global empirical models. This reduced amplification is believed to be partly due to true regional differences and partly due to limitations introduced by the 1D calculations and the local profile database. These differences are discussed in a supplementary document KYB_SI2.  Baram et al. (2020). IL empirical refers to the empirical models derived by Maiti et al. (2021) and the other two models are the host models used in that study  Seyhan and Stewart 2014). IL revised is the suggested V S30 scaling for Israel, with a reduced scaling (relative to CA) at long periods, keeping the IL scaling at short periods Finally, a revised linear V S30 scaling term is suggested (IL_revised in Fig. 3). The revised scaling is compared with the analytically derived linear V S30 -scaling terms (Baram et al. 2020) in response-spectral domain (Bea20 in Fig. 3), and with equivalent terms derived empirically for five regions, based on the NGA-West2 GM database . We maintain the IL analytical scaling up to T = 0.5 s, which is equivalent to 2 Hz, due to the general agreement in Fig. 3 at high frequencies and also because the short-period scaling of the IL models are very similar to other regions, specifically CA and Mediterranean (MED). In the long period range, we suggest a revised IL model, which is about half-way between the analytically-derived IL scaling (Bea20) and the CA scaling. This revised model represents a reduced dependence of amplification on V S30 at long-periods, while taking into account the limitations within both our analytical as well as our empirical models, shown in Fig. 2.

Finite-fault effects
Most modern GMMs include scaling of additional factors on top of magnitude, distance and site, such as mechanism, rupture depth, and hanging wall. These effects, referred to here as finite-fault effects, have been shown to affect ground motions, especially for large magnitudes at short distances-which is the range that is typically most important for hazard. Henceeven if these effects are not represented by the available data, they must be accounted for in forward engineering applications, which often extend beyond the range of available data. Because the Maiti et al. (2021) models do not include these effects (and are simply functions of magnitude, distance, and V S30 )-the synthetic database cannot be used to regress these factors because they are not represented within it. Furthermore, our empirical dataset includes only one single finite-fault model (for the 1995 Mw7.2 Nueiba earthquake, for which the closest record is at R rup = 92 km) and hence it cannot be used to constrain or validate these effects. Therefore-all finite fault effects have been fully adopted from ASK14, including their scaling coefficients, under the assumption that they are not sensitive to regional differences and should be quite stable between active crustal regions (Stafford et al. 2008). Because some of the components are correlated (e.g. mechanism and depth to top of rupture-Z TOR ), adopting them from the same host model (ASK14) is important, so that such unavoidable correlations are accounted for.
Hence, the following functions (and the associated ASK14 coefficients) are incorporated into our GMM: As in ASK14-the scaling of reverse faulting events is accounted for by the Z TOR term (a 11 = 0).
The hanging-wall (HW) effect in ASK14 was largely constrained by 3D finite-fault simulations (Donahue and Abrahamson 2014). It includes a factor a 13 , and five tapers to produce a smoothly varying HW effect as a function of the dip, magnitude, location over the rupture, depth, and distance off the ends of the rupture: (20) , h 1 = 0.25, h 2 = 1.5 andh 3 = −0.75. R y0 can be computed from R y0 = R x *|tan(Src2SiteA)|, where Src2SiteA is the azimuth from the source strike to the site, measured clockwise. R y0 can only be zero (for sites located along the rupture) or positive. If the R y0 distance metric is not available, the T 5 taper can be replaced using the following model: The effect of buried ruptures is accounted for through the Z TOR term and is capped at 20 km depth. While the average crustal depth in Israel is larger than that in CA (especially within the Dead-Sea basin), we did not have sufficient data to extend it to larger depths and hence maintain the exact same functional form and scaling. (11) for M ≥ 6.5 1 + a 2HW (M − 6.5) − (1 − a 2HW )(M − 6.5) 2 for 5.5 < M < 6.5 0 for M ≤ 5.5

Median model results
The model was regressed on the synthetic database in one step, due to the completeness and uniformity of the catalog. Regression residuals can be seen in Fig. 4 and as expected, there is no significant trend or bias and the variability is very low (Std < 0.1). Despite the limitation associated with the completeness and quality of the empirical database, it is crucial to confirm that the model is not rejected by the available data. Hence, model residuals were computed with respect to the empirical database (Yagoda-Biran et al. 2021). The residuals were then separated into between-event and within-event components, such that the between-event residuals were simply the median of all residuals associated with a specific event. Within-event residuals with respect to magnitude, distance and V S30 show very large variability but no observable bias or trend (Fig. 5) and hence those effects within the model are considered to agree with the data. However, the between-event residuals were found to be on average higher than zero for all spectral periods (see Figs. 6 and 7). This bias leads to an underestimation of the expected ground motions across the entire spectral range. Therefore, we adjust the constant a 1 with a fixed value, ∆ a 1 , such that the spectral shape is maintained, and the model agrees better with the available empirical data. This is a relatively large adjustment, which may be related to the poor quality of the data used to calibrate the simulations in MYK21, leading to a narrow bandwidth and very low values at high frequencies. Until the local empirical dataset is improved, this should be reflected in a relatively large epistemic uncertainty when using this model for hazard calculations.
While a random-effects regression is generally preferred for computing between-event residuals, in this case a random-effects regression yielded almost identical variability and slightly larger bias, but the resulting spectrum was not smooth and required additional adjustments. Therefore, we use a simple computation of the median residual per event, but limit the data used for computation to ensure that we are using events with higher quality. These limitations include a maximum rupture distance of 300 km, a minimum magnitude of 3.0 and a minimum number of records per event of 5 so that the median is representative. Furthermore, we only include recordings from the 'IS' network and discard our accelerometer data ('AA' network in flatfile), which is less reliable. The final model coefficients, including the a 1 adjustment, are presented in Fig. 8 and provided in Table 1 in the supplementary material KYB_SI1.
Finally, model magnitude and distance scaling are presented in Figs. 9 and 10, respectively. Both figures show the model response for a soft soil site, compared to seven other global models-four NGA-West2 models and three European models. Generally, the model trends are consistent with other global models. Specifically, at short spectral periods (T = 0.01 s), the model, KYB22 from hereon, is lower than other models at short distances and small magnitudes. At long spectral periods (T = 1 s), KYB22 is higher than other models, especially at large magnitudes and short distances. The reason for these differences is the spectral shape, as seen in Fig. 11-the KYB22 spectral shape is somewhat different than ASK14, for example, such that the peak is shifted to longer periods. For example, at a distance of 10 km-at a magnitude 5 KYB22 has a lower PGA value but similar T = 1 s value to ASK14, while at magnitude 7 KYB22 has similar PGA value but higher T = 1 s value than ASK14, due to the spectral shape.

Aleatory variability
The total aleatory standard deviation-σ-is computed as = √ 2 + 2 , where ϕ and τ are the within-event and between-event standard deviations, respectively. The between-and within-event standard deviations of a GMM are typically computed from model residuals obtained directly during random-effects regression. In this study, the model is regressed on a synthetic database, and the variability in the regression residuals is low and does not represent the true aleatory variability.
There are two options for computing σ in this case-one is to use the aleatory variability of the original empirical FAS model in Maiti et al. (2021). However-variability in FAS is consistently larger than in RSP, and strongly depends on smoothing. For example-the aleatory variability in Bayless and Abrahamson (2019) is much larger than that of ASK14 (2014), despite using the same database. The second option is computing the aleatory standard deviation (σ) from residuals between the Israeli RSP database and our RSP model, as explained in the previous section and shown in Fig. 12. The resulting standard deviations are very high, compared to typical values in global models. This is due to two reasons: first, σ has been shown to be magnitude-dependent, with increasing variability at small magnitudes (e.g. Abrahamson et al. 2014;Boore et al. 2014;Campbell and Bozorgnia 2014;Chiou and Youngs 2014). Because our empirical database is almost entirely comprised of small-magnitude data, the resulting variability is high and does not represent the variability at large magnitudes. This overestimation of σ is expected to dramatically increase the estimated hazard, which is driven by large magnitude events. The second reason for bias in σ could be due to errors in magnitude estimates. Such errors are mostly expected for small magnitudes, in which a finite-fault model does not exist and conversion between magnitude scales could lead to relatively large bias (see Yagoda-Biran et al. 2021). Errors in assigned magnitudes could significantly bias the model predictions and show up as increased between-event variability (Kuehn and Abrahamson 2018).
Due to the limitations listed above, we choose to adopt a global model for the aleatory variability which reflects the true magnitude-dependent nature of between-event  Spectral shape at three distances, compared to two other global models. Model names same as in Fig. 9 and within-event aleatory variability. The standard deviation terms are adopted from the ASK14 model, as follows: and As mentioned above, all model coefficients are provided in Table S1 in the supplementary material KYB_SI1. In addition, Table 1 below provides a list of the model coefficients, with an indication of the method used to constrain them-a regression on the synthetic database, the empirical data, the model ASK14 or other.

Hazard
In this chapter, the newly derived GMM, KYB22, is used within a forward engineering seismic hazard application, in order to compare its performance to other common alternatives, within the final framework in which it is to be used. State-wide probabilistic seismic hazard maps are produced, using different alternatives for ground-motion logic-trees, and the results are compared and discussed.
s 4 for M > 7 Fig. 12 Standard deviation of between-event and within-event residuals, using empirical data

Epistemic uncertainty
Epistemic uncertainty in GMM context is associated with both parametric uncertainty, as well as modeling uncertainty, in that the model is a simplified description of much complicated physical processes (Abrahamson 1990;Abrahamson and Bommer 2005). Capturing this uncertainty in SHA is therefore crucial for hazard estimations, as it expresses the limitations of our current knowledge. Epistemic uncertainty is incorporated in SHA using logic trees-each branch represents an alternative model and is assigned a relative weight. There are generally two approaches to building a GMM logic tree: (1) using multiple published GMMs, local or global, and assigning weights according to judgement, and (2) selecting a center model or models, referred to as the backbone model, and scaling them up and down, with weights representing a three-point normal distribution ). There is a growing interest in the backbone approach, which is being used more widely, because of its transparency regarding the level of uncertainty and the meaning of weights in the logic tree (Douglas 2018). We therefore use the new GMM-KYB22-as the center model in the backbone approach, and define the standard deviation for the scaled models, ∆.
The ∆ value derived by Maiti et al. (2021) for their nine GMMs for Israel, which KYB22 is based on, ranges between 0.23 and 0.46 natural log units. However, all the models in Maiti et al. (2021) are based on the same dataset, with very limited empirical data. Therefore, their derived ∆ values must be compared to other global estimates, to ensure these values are not underestimating the true epistemic uncertainty. Atkinson and Adams (2013) computed the value of ∆ for active crustal regions based on the spread of the NGA-West1 models for a range of magnitudes and distances. They suggest a distance-dependent relation for ∆, which ranges from 0.23 at short distances to 0.7 ln units at distances > 300 km. This higher value is controlled by data from M5.0 and M6.0, which were not well-constrained in the NGA-West1 project. Al Atik and Youngs (2014) evaluated the model-to-model variability within the NGA-West2 models and showed that it increases with period and distance. For strike slip faults, which dominate the seismicity in the Dead Sea region, the variability at short periods varies between 0.1 and 0.2 at distances < 100 km, increasing up to 0.4 for distances > 100 km. Weatherhill et al. (2020) show that short period ground motion epistemic uncertainty factors for several site specific studies fall within the range of 0.4-0.5 natural log units. Kowsari et al. (2020) analyzed events from Iran, and used 14 GMMs, both local and global. They found that the epistemic uncertainty at short periods ranges between 0.2 at short distances and 0.7 natural log units at distances > 100 km with an inverse dependence on magnitude (lower variability for higher magnitudes).
Based on the values obtained by Maiti et al (2021), compared to the studies mentioned above, we select a single value of ∆ = 0.4 natural log units for our analysis of hazard at T = 0.01 s, independent of distance or magnitude, which seems to average the different ranges and dependencies presented above.

Hazard runs with NSHMP-haz
We use the NSHMP-haz (v1) open source code, shared by the United States Geological Survey (USGS) on Github (see data and resources section). The seismogenic source model used for the hazard runs follows the magnitude-cutoff model presented by Maiti and Kamai (2020, Fig. 4(a) in their paper). The model is comprised of twelve finite fault sources, which are all part of the Dead Sea Fault (DSF) system (Fig. 13a). The off-fault seismicity, including the Cyprus Arc, is represented as gridded seismicity, using an adaptive kernel approach. The interaction between the on-fault and off-fault seismicity is defined as a magnitude cut-off, such that large events can only occur on the main faults, while ensuring balanced moment-rates. We run hazard analysis for Israel using three different approaches for GM representation: (1) using a single GMM with a full weight of 1.0, (2) using multiple GMMs with assigned weights, and (3) the backbone approach.
Six alternatives for GMM logic trees are used within our hazard analysis: (1) Campbell and Bozorgnia (2008) used as a single GMM, as currently used in many applications in Israel, including the current hazard map (from hereon CB08), (2) Campbell and Bozorgnia (2014) Fig. 13a, for 2% (left) and 10% (right) POE in 50 years. Squares are for Campbell and Bozorgnia (2008). Crosses are for Campbell and Bozorgnia (2014). Circles are for the Campbell and Bozorgnia (2014) used as a backbone with a standard deviation of 0.4. Upward pointing triangles are for 4 NGA-WEST2 models, equally weighted. Left pointing triangles are for the GMM derived in this study. Stars are the GMM derived in this study, used as a backbone with standard deviation of 0.4 of 0.4 and weight of 0.63 for the median and 0.185 for the low and high branches (KYB22 BB). Figure 13 presents state-wide hazard maps for PGA, 2% probability of exceedance (POE) in 50 years, calculated using these six alternatives. Inspecting Fig. 13, it is clear that hazard is lowest in Fig. 13(a), using Campbell and Bozorgnia (2008), as is currently used in the Israeli building standard. Figure 13(d) and (f) show comparable results of the four NGA-West2 GMMs and the GMM derived in this study, KYB22, used as backbone.
In Fig. 14 we compare hazard results at six locations, for two spectral periods-PGA and 1 s, and two hazard levels-POE of 2% and 10% in 50 years. The six locations, shown on the map in Fig. 13(a), are six large cities in Israel, located at varying distances from the main faults. As can be seen in Fig. 14, for PGA (top row)-the locations close to fault sources, namely Haifa (close to the Carmel Fault), Tiberias and Eilat (close to the DSF), show higher variability in the hazard results, than the locations further away, namely Jerusalem, Tel Aviv and Beer-Sheva. The two GMM combinations of KYB22 BB and NGA_ West2 show very similar hazard results (mostly between CB08 and CB14) both close toand further away from faults, at both POEs. At long periods (T = 1 s, bottom row)-KYB22 is consistently higher than all other model combinations, with the greatest differences, of up to 200%, at locations close to the DSF-Eilat and Tiberias. This overprediction is consistent with what was observed and discussed in Figs. 9, 10 and 11.
Naturally, further inspections of the new GMM in different applications are required, but the results presented above show that the epistemic uncertainty is captured and that the hazard results are comparable and similar to what other alternatives of popular GMMs would yield. Therefore, we suggest using this new GMM in hazard calculations for Israel, as a model within the ground motion logic tree, or as the center model used in the backbone approach, mostly because it has the advantage of being validated and calibrated to local seismicity, while also providing full applicability.

Summary
In this paper, we present a new hybrid GMM for Israel. The GMM is hybrid in that it relies on a suite of nine FAS GMMs, developed using both recorded data and simulations, and the database for the regression is a synthetic one, converted to RSP using random vibration theory. It is also hybrid in that some terms in the GMM are fully regressed while others are adopted from the parent model ASK14, because the empirical data is too scarce to constrain them. The model residuals are computed with respect to the empirical database, and the model is slightly altered so that it would not be rejected by the available data.
The new GMM is then compared to other models and GMM combinations in a statewide hazard analysis. We find that hazard results are comparable and similar to what other GMMs yield, and that epistemic uncertainty is captured.
This newly derived GMM is ready to be used in forward engineering applications. It covers the full response-spectral domain from 0.01 to 10 s, includes nonlinear site-effects as well as finite-fault effects, and is calibrated to local Israeli data. Hence, we strongly suggest that this model be considered as one of the branches in the ground-motion logic tree when evaluating seismic hazards for Israel.