A backbone seismic ground motion model for strike-slip earthquakes in Southwest Iceland and its implications for near- and far-field PSHA

The earthquake hazard and seismic risk in Iceland are highest in the Southwest due to the transform faulting in the South Iceland Seismic Zone (SISZ) and Reykjanes Peninsula Oblique Rift (RPOR) being in close proximity to a large part of the population. Reliable probabilistic seismic hazard assessment (PSHA) is therefore critical in this region which in turn requires two of its key elements: the most appropriate ground motion models (GMMs) and the specification of seismic sources. In this study, we address this by employing a suite of new hybrid Bayesian empirical GMMs and a new physics-based finite-fault system model for the SISZ–RPOR. By ranking the GMMs using the deviance information criterion against the Icelandic strong-motion dataset we propose backbone GMMs (i.e., the most appropriate models), which along with their backbone Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Delta }$$\end{document}-factors not only comprehensively capture the characteristics of Icelandic data but also the epistemic uncertainties of their ground motion prediction. We then simulate suites of synthetic finite-fault earthquake catalogues that are consistent with the physics-based fault system and long-term seismic activity of the region. We carry out a Monte-Carlo PSHA for two representative near-fault (the town of Selfoss) and far-field (the capital city of Reykjavik) sites in Southwest Iceland, and compare the results with those of a classical PSHA. We show that the backbone PSHA is consistent with classical PSHA point estimates but more importantly, the backbone approach reveals the “body” of the hazard i.e., the majority of realistically expected hazard values, which dwarfs any differences in the results from the two approaches. The PGA results show that the body of hazard values varies from 0.07-0.13 g and 0.4-0.6 g for Reykjavik and Selfoss, respectively. The scale factor is distance- and period-dependent that results in a constant body of the backbone hazard values for a far-field site at all periods while at long periods and decreasing annual exceedance rates, the body progressively increases. Moreover, the disaggregation showed that the modal event for a 475 year return period is a Mw ~ 6–6.5 at distance ~ 14–26 km for the far-field site, but at ~ 6 km for the near-fault site. We conclude that the backbone approach and the models used in this study make ideal candidates for improved PSHA for Iceland.


Introduction
Iceland is the subareal manifestation of a localized oceanic plate uplift due to the Iceland hot spot. Its collocation and interaction with the Mid-Atlantic Ridge, the extensional tectonic plate margin that separates the Eurasian and North American Plates drives the intense volcanic and seismic activity of this most seismically active region in northern Europe. The Mid-Atlantic Ridge crosses Iceland from Southwest to north but due to a ridge-jump caused by the hot spot, the ridge is displaced towards the east, forming two major transform zones, the South Iceland Seismic Zone (SISZ) and the Reykjanes Peninsula Oblique Rift (RPOR) in the Southwest and the Tjörnes Fracture Zone (TFZ) in northwest. They are the regions with the greatest potential for occurrence of large earthquakes in Iceland. The SISZ is collocated with the South Iceland Lowland, a relatively flat and populated agricultural region with all critical infrastructures and lifelines of a modern society such as hydroelectric and geothermal power plants, dams, above-ground pipelines and electrical transmission lines. The RPOR however is a relatively remote high basaltic plain (Thordarson and Hoskuldsson 2002).
In these zones, strong earthquakes have repeatedly taken place on a dense array of separate but parallel near vertical north-south trending right-lateral strike slip faults that in fact is continuous from the SISZ across the Hengill triple junction and to the relatively uninhabited RPOR (Einarsson 1991(Einarsson , 2008(Einarsson , 2014Einarsson et al. 2020;Steigerwald et al. 2020). As a result, the earthquake hazard is the highest in this region, the SISZ in particular, and the investigation of earthquake strong-motion and its effects on manmade environment is of particular interest and vital for seismic risk assessment and its mitigation (Tryggvason et al. 1958;Sigbjörnsson et al. 1995;Solnes et al. 2004). For this purpose, probabilistic seismic hazard analysis (PSHA) is used to quantify the probability of levels of earthquake ground motion parameters being exceeded over a specified time period at any given location (e.g., Cornell 1968;Reiter 1991;McGuire 2004).
In order to produce a reliable PSHA, a careful estimation of its inputs and their uncertainty is of paramount importance. The seismicity parameters and selection of appropriate ground motion models (GMMs) are the key PSHA inputs affecting the hazard levels (e.g., Cramer et al. 1996;Cramer 2001;Lombardi et al. 2005;Sabetta et al. 2005). A rigorous assessment of the uncertainties of these inputs is therefore required. In the scientific literature, such uncertainties are divided into two major categories, aleatory variability and epistemic uncertainty (Toro et al. 1997). The aleatory variability is related to the natural unpredictability of the earthquake process whereas the epistemic uncertainty arises from incomplete knowledge and limited data. In the current practice of PSHA, the different estimates of the median ground motions predicted by empirical GMMs are attributed to epistemic uncertainty. In this regard, classical PSHA deals with the epistemic uncertainty in ground motion estimates using a logic tree framework (Kulkarni et al. 1984). However, the logic tree results are highly dependent on how the branches are weighted and populated, which is subjective to expert judgment (Bommer 2012;Delavaud et al. 2012). In addition, the use of GMMs in logic trees cannot describe properly the distribution of ground motions and as a result, the mean hazard curves do not contain any specific information about the lower or higher ground motion estimates (SSHAC 1997;Kotha et al. 2018).
To overcome this problem, a backbone approach to PSHA has been introduced to capture better the center, body, and range of expected ground motions (Atkinson and Adams 2013). In this approach, the median GMM (i.e., the backbone model) represents the center and can be scaled up and down with a factor to capture the body of realistically expected ground motion values i.e., the majority of the distribution of ground motions (while the range would represent nearly the entire distribution). This approach is fast becoming the standard-practice in PSHA as it has been applied in the US national seismic hazard maps (Petersen et al. 2008), for the PSHA underlying the Canadian National Building Code (Atkinson and Adams 2013), for the NGA-West2 models (Al Atik and Youngs 2014), for the European and Middle eastern data (Douglas 2018a, b) and for the Iranian strongmotions (Kowsari et al. 2020a).
However, the application of the recently proposed backbone approach in PSHA has been nonexistent in Iceland. Therefore, in this study, we propose the first backbone seismic ground motion model for strike-slip earthquakes in Southwest Iceland and show examples of its PSHA results in the near-and far-field regions, respectively. In doing so, we take advantage the new 3D finite-fault system model proposed by Bayat et al. (2022) for the bookshelf strike-slip fault system in Southwest Iceland. On its basis, we carry out multiple Monte Carlo simulations of finite-fault earthquake catalogues that are consistent not only with the overall rate of seismicity in the region, but for individual subzones of the system as well, in accordance with the fault system model. The Monte Carlo simulation approach is consistent with previous local PSHA studies in Iceland where this approach had been used long before (Sigbjörnsson et al. 1995) this approach gained popularity (Atkinson 2012) (albeit only for point sources, and empirically estimated activity rates). We also take advantage of the new Bayesian GMM (Kowsari et al. 2020b), establish the first backbone GMM for Iceland thus reducing the epistemic uncertainty in GMM predictions of Icelandic strong-motions and then investigate the effects on PSHA using the backbone approach vs. the classical approach for two key population centers in Southwest Iceland i.e., one in the near-fault region (Selfoss) and the other in the far-field region (Reykjavik). Specifically, we apply a recently developed data-driven GMM-ranking method on several candidate GMMs developed from local, regional, and worldwide datasets (Kowsari et al. 2019b). The datadriven ranking methods help to rank the GMMs in an objective way and thus reduces the epistemic uncertainty associated with the selected GMMs. Thus, the GMMs are ranked and the best GMM is chosen as the backbone model. Then, a scale factor that can be added and subtracted from the backbone GMM to cover the body of ground motions is inferred from different GMMs. Finally, we carry out a backbone PSHA at the two sites in Southwest Iceland using the simulated finite-fault earthquake catalogues for the entire SISZ-RPOR which is based on a new, physics-based 3D fault system model of Southwest Iceland that explains the observed Icelandic earthquake catalogues. We discuss the differences in the results with those of the classical PSHA in terms of the potential implications that the use of the backbone approach is expected to have on a reappraisal of seismic hazard in Iceland.

Study region and strong-motion data
The study region is Southwest Iceland where the transform motion of the eastward jump of the Mid-Atlantic Ridge is taken up by a complex "bookshelf" faulting system i.e., an array of near vertical dextral strike-slip faults across the SISZ and the RPOR (Einarsson et al. 1981(Einarsson et al. , 2020Einarsson 1991Einarsson , 2008Einarsson , 2014Steigerwald et al. 2020). This bookshelf system characterizes the occurrence of strong earthquakes in the SISZ and has recently been shown to be continuous across the Hengill Triple Junction and along the RPOR (Steigerwald et al. 2020). There are five volcanic systems in the RPOR in which normal faulting seismicity occurs episodically during (very rare) intense rifting episodes 1 3 between which transcurrent motion dominates (Saemundsson et al. 2020). Moreover, the corresponding maximum earthquake magnitudes associated with rifting episodes are believed to be approximately M w 4.5-5, and as such they are considered to have a minimum contribution to the seismic hazard. In this region, the hypocentral depth of earthquakes gradually increases from west to east. It is generally at 1-5 km in the western part of the RPOR and down to 12-15 km depth in the easternmost part of the SISZ (Stefánsson and Halldórsson 1988;Stefansson et al. 1993;Panzera et al. 2016;Steigerwald et al. 2020). For the strike-slip faults in this region, the maximum magnitude is estimated to vary from ∼ M w 5.5 in the westernmost part of the RPOR to ∼ M w 7-7.2 in the easternmost part of the SISZ. Thus, the seismogenic potential increases towards the east as confirmed by the historical catalogue and observational reports in annals (Tryggvason et al. 1958;Halldórsson et al. 1984;Stefansson et al. 1993;Ambraseys and Sigbjörnsson 2000).
In this study, we used the earthquake strong-motion data recorded on 30 stations of the Icelandic strong-motion network from six strike-slip mainshock earthquakes that occurred in the SISZ between 1987 and 2008 covering magnitude range of 5.1-6.5, with the distance range being 1-80 km (Sigbjörnsson et al. 2014). In addition, we used the extreme near-fault array recordings of ICEARRAY I during the 2008 Ölfus earthquake Halldorsson and Sigbjörnsson 2009). In Iceland, the recording stations are classified according to the Eurocode 8 classification scheme into the rock class ( V S30 > 800 m/s) and the stiff soil class (360 m/s < V S30 < 800 m/s) (Sigbjörnsson et al. 2014). Figure 1 shows the distribution of the earthquake magnitude versus the Joyner-Boore distance, R JB (i.e., the closest horizontal distance to the vertical surface projection of the fault surface) of the dataset for the two classes of rock and stiff soil, respectively. We note that in this study we focus explicitly on sites classified as rock, as they constitute the majority of the dataset and are associated with greater prediction confidence than those on stiff soil (Kowsari et al. 2020b;Rahpeyma et al. 2022).

Fig. 1
The distribution of the earthquake magnitude versus distance of the Icelandic strongmotion dataset for the two site classes of rock and stiff soil, respectively. The vertical and horizontal bars show the number of observations with respect to distance and magnitude, respectively 1 3

Empirical ground motion models
For the simple characterization of the propagation of seismic ground motion peak parameters with increasing source-distance, we consider several empirical GMMs in three categories: local models calibrated to the earthquake strong-motion data of South Iceland, regional models corresponding to European and Middle Eastern data (primarily those that had been proposed in the previous hazard studies as being suitable models for the seismotectonic environment of Iceland) and finally, as an example, NGA-West2 models developed from worldwide data of shallow crustal interplate earthquakes. Other GMMs that fall into these categories but fail the criteria proposed by Cotton et al. (2006) and Bommer et al. (2010) for the minimum requirement of the functional form for GMMs to be applied in PSHA are omitted. Specifically, we retain GMMs with functional forms that are able to capture the prevalent non-self-similar scaling of ground motion at large amplitudes, in particular in the near-fault region (see e.g., Abrahamson and Shedlock 1997;Halldorsson and Papageorgiou 2005;Abrahamson et al. 2008;Douglas 2018c). The final GMMs are the following: The ones calibrated in Kowsari et al. (2020b) to the Icelandic strong-motion dataset using Bayesian inference with informed priors, named as Kea20-1 to Kea20-6. Each of them has a different functional form, which satisfies the abovementioned criteria, and together they facilitate the estimation of epistemic uncertainty. Then, the regional GMMs such as Ambraseys et al. (2005) Table 1 lists the key metadata of the GMMs used in this study, their magnitude and distance ranges, the spectral period band for pseudo-spectral acceleration (PSA), as well as the inclusion of peak ground acceleration (PGA) and/or peak ground velocity (PGV) and region of origin.
Some of the GMMs such as AB10 and Zh06 were proposed in the ESHM13 as suitable GMMs for PSHA in active oceanic regions. The Am05 that was developed based on European and Middle Eastern dataset, including the Icelandic data available at the time, has been applied in several PSHA studies in Iceland. While it had already been shown that the NGA models cannot fit Icelandic ground motions (e.g., Ornthammarath  Figures 2 and 3 show the attenuation of the selected GMMs versus distance for PGA and PSA at T = 0.3 , 1.0 and 2.0 s for M w 5.5 and 6.5 on rock, respectively. The GMM predictions for rock use V S30 of 800 m/s (see Sigbjörnsson et al. 2014) and the different ground motion intensity measures of the GMMs have been converted to the rotation invariant measure (Rupakhety and Sigbjörnsson 2013) by using the numerical comparison of intensity measures shown in Kowsari et al. (2020b). The color scheme shows the different groups of GMMs: Yellow-to-orange colors represent the local Kea20 1-6 models, blue colors represent NGA-West2, and green colors represent regional GMMs and the purple is Zh06. The red circles are the observed strong-motion data in the two magnitude bins, dominated by recordings of the M w 6.5 and M w 6.4 earthquakes in 2000 and the M w 6.3 2008 earthquakes. From this figure, it is qualitatively clear that the Kea20 models fit the recorded data well in the magnitude and distance range where data is available. The regional and NGA-West2 GMMs appear to exhibit a strong bias against the Icelandic strong-motion data primarily due to the failure of capturing the high near-fault motions and rapid attenuation of strong-motion (Kowsari et al. 2019a, Fig. 2 The attenuation of the candidate GMMs versus distance for PGA and PSA at T = 0.3 , 1.0 and 2.0 s for M w 5.5 on rock. The red circles are the observed strong-motion data 2020b). As a result, they underestimate the Icelandic strong-motions in the near-fault region and overestimate it in the far-field region.

Ranking of ground motion models
Different ground motions predictions by different GMMs (see e.g., Figs. 2 and 3) have the largest contributions to the overall epistemic uncertainty, in particular at long return periods ). Therefore, selection of appropriate GMMs to reduce the epistemic uncertainty is vital for application in PSHA and this can be achieved by using data-driven GMM-ranking methods. The data-driven methods reduce subjectivity in the selection process and thus have an advantage over the classical residual analysis method and other statistical methods such as chi-square test, Kolmogorov-Smirnov statistic, Pearson correlation and variance reduction (Scherbaum et al. 2004).
The log-likelihood (LLH) method proposed by Scherbaum et al. (2009) and Euclidean distance ranking (EDR) method proposed by Kale and Akkar (2013) are at present the most used data-driven ranking methods. In some case, however, their performances have exhibited shortcomings. Namely, the LLH method prefers the predictive model with larger sigma when Fig. 3 The attenuation of the candidate GMMs versus distance for PGA and PSA at T=0.3, 1.0 and 2.0 s for M w 6.5 on rock. The red circles are the observed strong-motion data the observed data congregate away from the median estimations of the two GMMs (Kale and Akkar 2013). On the other hand, the EDR method favors a smaller sigma when two predictions give the same mean, regardless of what the true uncertainty is Mak et al. (2014). Thus, in those cases, sigma of the GMM may lead to ignore the model bias (represented by deviation of the median). Therefore, the selected GMMs for performing either non-ergodic or partially non-ergodic PSHA must accurately represent the aleatory variability of the ground motions in the region under study. Addressing these issues, Kowsari et al. (2019b) proposed a datadriven method using the Deviance Information Criterion (DIC) for selection of the most suitable GMM for application in PSHA. The main advantage of the DIC is to introduce the posterior sigma as a key unknown measure in order to rank the models objectively. The posterior distribution of sigma is then obtained based on the observed ground motions and shows the deviance of predicted values from the observed ground motions that is representative of the aleatory variability in the region under study. In this way, the DIC ranks models more favorably when they are associated with smaller bias and the corresponding posterior sigma is close to the aleatory variability of the ground motions in the region under study, for the given dataset. In the context of the Bayesian statistical framework, the posterior sigma is conditioned on the observed ground motions obtained from the region under study. Since the posterior sigma considers the misfit of the GMM predictions to the observed ground motions, it partially removes the ergodic assumption.
In Bayesian statistics, the posterior distribution describes updated information about unknown parameter ( ) and can be obtained by multiplying a prior distribution to a likelihood function as follows: where p( ) is the prior distribution and p(y| ) is the likelihood function that is a probability distribution that expresses the contained information in the observed data. Let us assume that sigma is unknown ( = ) and that the posterior mean of sigma will be based on the ground motion observations of the region where the GMM model selection is to be carried out for PSHA purposes. In this way, we assume a scaled inverse chi-squared distribution for the prior distribution of sigma (Kowsari et al. 2019b): where is the number of chi-squared degrees of freedom, and s 2 is the scaling parameter. On the other hand, for the likelihood function we assume that the logarithm of the ground motion parameter follows a normal distribution: where N is the number of observations, y is the natural logarithm of the observed ground motion parameter, ( ) is the mean value predicted by the GMM, and is the standard deviation of the GMM. Therefore, the posterior distribution is given by (Kowsari et al. 2019b): The conditional posterior distribution of the sigma is again a scaled inverse chi-squared distribution with the degrees of freedom, N + , and the scaling parameter, . As a result, the DIC of the normal model is given by (for more details see (Kowsari et al. 2019b)): where L is the number of samples, are drawn samples from the posterior distribution using a Markov Chain Monte Carlo method, and −2 is the posterior mean value of sigma squared.

A new backbone ground motion model for Southwest Iceland
A candidate GMM that most appropriately describes Icelandic strong-motions needs to be identified and taken to represent the backbone GMM for application in PSHA for Iceland. Thus, to evaluate the backbone GMM's objectively, we use the DIC data-driven selection method presented above that uses a Bayesian statistical method for inferring the most suitable GMM. Figure 4 shows schematically the results of the DIC ranking against the strongmotion data, presented for PGA and PSA at different oscillator periods of T = 0.3, 1.0, and 2.0 s for the two site classes (i.e., rock vs. stiff soil). We indicate the results of the ranking by the arrangement of these 15 GMMs from 1, the best ranked, to 15, the least ranked, from left to right. In addition, for clarity the color scheme in the ranking results follows that assigned to the GMMs in Fig. 2 with yellow-toorange colors representing the local Kea20 1-6 models, blue colors representing NGA-West2, and green colors representing regional GMMs. Not surprisingly, the results show that the Kea20 GMMs are ranked best almost at all periods (except for PSA at T=2.0 s at stiff soil). The ranking also indicates that NGA-West2 models are much less favorable Ranking of the selected GMMs against the Icelandic strong-motions using the DIC method for PGA and PSA at different periods of T = 0.3, 1.0, and 2.0 s for rock (top) and stiff soil (bottom) site classes. The color scheme shows the different groups of GMMs: Yellow-to-orange colors represent Kea20, blue colors represent NGA-WEST2, and green colors represent regional GMMs and the purple is Zh06 GMMs to describe Icelandic ground motions. More interestingly however, in most cases the regional GMMs that had been applied in the previous seismic hazard studies in Iceland were shown to represent the least favored models. This is consistent with recent results showing that they are strongly biased against the Icelandic dataset and therefore, the results of previous seismic hazard studies in Iceland using these GMMs are cast in serious doubt (Kowsari et al. 2019a). We note that the limitations of the Icelandic strong-motion dataset, in particular the lack of larger magnitude data above M w 6.5, may appear to allow questioning these results. This is not the case however, as (1) the most prevalent characteristics of Icelandic earthquake strong-motion are the high near-fault amplitudes and the rapid attenuation with distance. None of the other candidate GMMs are able to accommodate this in a satisfactory manner and only the locally calibrated GMMs (i.e., Kea20) therefore captures this characteristic. Hence, it is the backbone model. (2) The Kea20 models all have been calibrated using Bayesian inference with informative priors for those GMM-parameters that were essentially unconstrained due to the above limitation. They are thus essentially hybrid models, constrained by the available Icelandic strong-motion data, and the universally observed GMM behavior at larger magnitudes i.e., the near-fault saturation of peak amplitudes, also consistent with physics-based finite-fault modeling of earthquake rupture and strongmotion simulations (see Kowsari et al. (2020b) for details and references).

Epistemic uncertainty of the ground motion models
Having established that one of the empirical Bayesian GMMs of Kea20 is the backbone model at the periods considered in this study (except for PSA at T=2.0 s at stiff soil), the corresponding scale factor Δ needs to be inferred for each oscillator period. For this purpose, we assume the scale factor can be the model-to-model variability that is the standard deviation in the logarithm of the predicted median ground motion, which is proposed by Al Atik and Youngs (2014) as: where log(PGA) i is the logarithm of the median ground motion predicted by the ith GMM, and w i is the probability weight assigned for the ith model. Although a limited Icelandic strong-motion data does not allow us to capture the total epistemic uncertainty, the modelto-model variability can provide useful information about the strong-motion estimates at different magnitudes and distances.
Therefore, we assigned different weights to each model to evaluate the standard deviation of the predicted median for the recalibrated GMMs at five different magnitudes between M w 5-7. We assign a weight of 0.2 for the first ranked model, 0.1 for the second and third ranked GMMs and 0.05 for the rest. We note that the results shown are not overly dependent on these weight values, because the top three ranked models (Kea20 models) have similar predictions and thus ranking scores and are as a result assigned most of the weights. Figure 5 shows the model-to-model variability versus distance at different earthquake magnitudes.
From Fig. 5 (top row), the model-to-model variability grows with distance that is in accordance with the analysis of the NGA dataset of ground motions of shallow crustal events in active tectonic regions (Atkinson and Adams 2013). It generally has minima in the distance range of 5-20 km and increases considerably in the near-fault ( R < 5 km) and far-field ( R > 40 km) regions. Moreover, the variability is highest for smaller events as their relatively much smaller source dimensions allow them to have much greater variability in hypocentral depths, resulting in greater variability in the near-fault motions observed on the surface (Sonnemann et al. 2020). The model-to-model variability is larger at short periods (i.e., PGA and PSA at T=0.3s) than the longer periods (PSA at T=1.0 and 2.0 s). For PSA at T=1.0s, the small standard deviation indicates that the variability between the predicted ground motions from selected GMMs is relatively low. Therefore, the model-to-model variability, can be considered as the scale factor of the backbone GMM. The scale factor can be thought of being representative of the epistemic uncertainty that the GMMs contribute to the total epistemic uncertainty, the latter being unattainable due to the limited dataset that is available. The backbone GMM along with its corresponding Δ are shown in Fig. 5 (bottom row) versus distance for four oscillator periods, on rock and stiff soil, respectively. The backbone model predictions for PGA and PSA at T=0.3, 1.0 and 2.0 s are thus shown as solid black lines and the backbone ±Δ by a gray shaded region around the mean. The red circles are the observed ground motions from the three largest earthquakes (the M w 6.5, 6.4 earthquakes of 2000 and M w 6.3 in 2008) in the SISZ that account for the majority of the data. The corresponding weighted moment magnitude of the data shown is M w 6.4. We see that the backbone GMM±Δ not only appropriately captures the Icelandic strong-motion over the available magnitude and distance ranges but also, the other top-ranked GMMs fall within the ±Δ . The results show that the scale factors slightly increase with distance which is consistent with the results presented in Kowsari et al. (2019a) where the variability among the median ground motion estimates of the Icelandic GMMs was shown.

The SISZ-RPOR earthquake source model
Having now established the backbone GMM and before carrying out any PSHA, we need to define and categorize the seismic sources in the best and most reliable, yet effective, way possible based on all available information. Figure 6 shows the SISZ-RPOR along with the short-term revised instrumental earthquake catalogue compiled by Panzera et al. (2016) for the very active period of 1991-2013 (black dots) and the recently revised historical earthquake catalogue of significant earthquakes since 1904, the ICEL-NMAR (circles) (Jónasson et al. 2021). In addition, mapped and inferred surface fault traces of significant strike-slip faults in the RPOR are also shown based on surface fault mapping, microearthquake relocations and from historical accounts (Roth 2004;Clifton and Kattenhorn 2006;Hjaltadóttir 2009;Einarsson 2014;Steigerwald et al. 2018;Einarsson et al. 2020). The system of north-south trending strike-slip faults in the SISZ-RPOR dominates its strain release during long periods devoid of intense rifting episodes (Saemundsson et al. 2020). For PSHA purposes therefore the relevant source zone in Southwest Iceland is the transform fault system of SISZ and RPOR. Notably, the maximum earthquake magnitudes vary in a systematic manner along the zone and is a manifestation of the systematic tectonic variations across the zone. As a result, this seismic region should not be considered to be a single source zone, but a more realistic zonation as shown in Fig. 6 is more appropriate. Nevertheless, the very scarce earthquake catalogue for some of these subzones all but precludes a reliable estimation of zone seismic activity rates.
Circumventing this issue by approaching it from a physics-based standpoint, Bayat et al. (2022) proposed a new physics-based finite-fault system model for the transform fault systems of the SISZ and the RPOR. Constrained by known historical and mapped fault locations, historical and instrumental seismicity, geophysical information e.g. on maximum seismogenic depth in the region, they modelled the SISZ-RPOR transform zone as an array of finite-size "bookshelf" strike-slip faults of systematically varying seismic potential and accounted for variability in interfault distances guided by fault Fig. 6 The South Iceland Seismic Zone (SISZ) and the Reykjanes Peninsula Oblique Rift (RPOR) along with seismicity from the earthquake catalogue compiled by Panzera et al. (2016). The vertical lines in the SISZ are the suggested source fault locations of historical earthquakes based on fault mapping and/or investigation of historical records SISZ (Roth 2004;Einarsson 2015). In the RPOR, the lines show the strike-slip faults based on Clifton and Kattenhorn (2006), Steigerwald et al. (2020) and Einarsson et al. (2020). The circles are the ICEL-NMAR catalogue (Jónasson et al. 2021) and the small black dots are the short-term catalogue of Panzera et al. (2016) mapping by field surveys and relative relocations of seismic swarms (e.g., Hjaltadóttir 2009;Einarsson et al. 2020;Steigerwald et al. 2020). The activity of the fault system is calibrated to the velocity of plate tectonic extension in the region as reported by long-term deformation studies (Sigmundsson et al. 1995;Hreinsdóttir et al. 2001;Sigmundsson 2006;Decriem et al. 2010;Einarsson et al. 2020). The resulting model is completely specified in terms of suites of 3D finite-fault strike-slip faults in the zone along with their corresponding annual slip rates, thus fully capturing the salient tectonic characteristics of the transcurrent faults in the region that dominate the seismic hazard.
Carrying out a full physics-based PSHA from earthquake rupture and ground motion simulations is outside of the scope of this study. Instead, we utilize the characteristics of the fault system models to develop area seismic zones with seismicity rates that are constrained by the physics-based fault system in each zone . In that study, new estimates of the seismicity rates of each zone (Fig. 6) that are not dependent on the existing catalogue were provided. For this purpose, first the distinct faults of the fault system model that fall inside each zone were identified. On the basis of their total annual slip and moment rates in each zone, the corresponding magnitude-frequency relationship was estimated using an iterative version of the method proposed by Brune (1968) that estimates the fault slip rate from a magnitude-frequency relationship. This entails providing initial values of a and b on the basis of past PSHA studies (that were based on statistical analysis of the catalogue), estimating the corresponding total slip rate in the zone, and then in an iterative manner allowing them to vary over realistic ranges until the magnitude-frequency relationship for all zones reproduces the slip rates of the physical model. Therefore, from the complete description of the physical model, they developed simple and equivalent Gutenberg-Richter models for subregions of the transform zones (something that the limited earthquake catalogues never could produce) which in turn allows a more accurate approach to standard-practice PSHA. Table 2 presents the resulting Gutenberg-Richter parameters of the magnitude-frequency relationship along with the maximum earthquake magnitudes from the physical model for each subregion. These results provide us with a much greater confidence in the zone seismicity rates over those obtained by analyzing the incomplete catalogue.
We took advantage of the availability of these models and simulated multiple finitefault earthquake catalogues for the entire bookshelf system that are both compatible with the earthquake faulting and the long-term seismicity in the region. One realization of the synthetic 3D fault system is shown in Fig. 7. It models 3000 years of strike-slip fault activity in the region limited with a minimum magnitude of M w 5. The magnitudefrequency distribution of these synthetic earthquakes in each zone is consistent with the activity rate derived from the 3D fault system model for each zone that is presented in Table 2. This figure shows one realization, the vertical surface projection of these randomly distributed fault planes with colors indicating magnitude bins (minimum magnitude of 5), clearly reflecting the expected increase in maximum possible magnitudes from west to east. M w,max~ 5.50~ 6.00~ 6.50~ 6.50~ 6.70~ 7.00

Near-fault and far-field PSHA from backbone vs. conventional models
To illustrate the differences between PSHA results using the backbone approach vs. the classical approach, we carry out a site-specific PSHA on rock for two important population centers, one in the near-fault region (i.e., the town of Selfoss) and one in the far-field (i.e., the capital city of Reykjavik). The former is located in the center of zone 5 while the latter is located outside of the seismic zone region but closest to zone 3. For the classical PSHA, we used the new set of Icelandic GMMs (i.e., Kea20) in a logic tree framework with equal weights. In the backbone approach we applied the first-ranked models shown in Fig. 3 as the backbone GMMs. Thus, the logic tree branches are populated with different backbone GMM and the corresponding ±Δ at different periods. The weights of 0.17, 0.66 and 0.17 are respectively assigned for the lower, central and upper values of GMM estimates based on the discrete approximation of the normal distribution when Gaussian quadrature is applied (Miller and Rice 1983).
In this study, we carry out PSHA using a Monte Carlo approach (Kowsari et al. 2021). For this purpose, we use the multiple finite-fault earthquake catalogues simulated for the entire bookshelf system . Then, the horizontal distance to the vertical surface projection of the fault, R JB , known as Joyner-Boore distance is calculated for each discrete earthquake line source (i.e., simulated fault) that affects the chosen site. Moreover, we generate the normalized residual (ε) that represents a measure of the goodness-of-fit of GMMs which is generally assumed to be normal with a mean zero and a standard deviation (σ). Then, we estimate the ground motion intensity measure of interest using the selected GMMs for each generated set of magnitude, distance and epsilon. This process is repeated for a specified number of subcatalogues and then the worst-case ground motion from each of these subcatalogues is selected and sorted. Finally, the ground shaking values for earthquakes corresponding to different hazard levels can be determined from the sorted list.
The PSHA results from the classic and backbone approaches for the two sites are shown in Fig. 8 in terms of hazard curves at four different periods. The hazard curves for Reykjavik are solid while those for Selfoss are dotted. The color-coding corresponds to which Fig. 7 One realization of a hypothetical 3D fault system model for each zone consistent in activity with the magnitude frequency distribution of the zone, as derived from the 3D fault system model realizations. This synthetic finite-fault catalogue results in a collective zone-average fault slip-rates that is consistent with that postulated by the 3D fault system model. The locations of the faults are allowed to vary randomly, and the colors indicate the distribution of magnitudes approach is employed. Moreover, the black dashed lines show the 15% and 85% fractile hazard curves from the classical approach. The horizontal blue dashed lines show two hazard levels that are important for seismic design purposes i.e., earthquakes corresponding to 10%-in-50 years and 2%-in-50 years probability of exceedance that has a return period of 475 years and 2475 years, respectively. The results show that the mean hazard curves using the backbone and classic approach are similar for both sites, which was to be expected. However, the gray shaded areas around the backbone (i.e., "central") hazard curve show us the "body" of the results i.e., the variation of the hazard values associated with the mean ±Δ effects when applied to the backbone model. In other words, the body shows us the range of realistically expected hazard values for the site. That inevitably leads the distribution of backbone hazard estimates for any given annual exceedance rate to be symmetric around the central backbone hazard curve, and at all oscillator periods. We observe that while the body of the backbone hazard values for a far-field site is relatively constant for Fig. 8 The seismic hazard curve for Reykjavik and Selfoss at PGA and PSA at T = 0.3, 1.0 and 2.0s. The horizontal blue dashed lines show the annual exceedance rate of 0.0021 and 0.0004 that is corresponding to an earthquake with a return period of 475 years and 2475 years, respectively. The gray shaded area shows the uncertainty in the backbone PSHA representing the lower and upper predictions. The black dashed lines show the 15% and 85% fractile hazard curves for the classical approach a given annual exceedance rate at all periods, at long periods for a near-fault site it is seen to progressively increase with decreasing annual exceedance rate. This is a direct result of a much larger Δ at small distances at long periods than observed at short periods (see Fig. 5), emphasizing that the scale factor is distance-and period-dependent and indicating that the contribution to the hazard at longer periods would be dominated by very close-by earthquakes. Then, at short periods (PGA and PSA at T = 0.3s) the body of hazard values is observed to be wider at far-field sites than at near-fault sites and at all annual exceedance rates. As shown in Fig. 5, the Δ values available for the far-field site are larger than for the near-fault site. The reverse is observed at long periods due to the same reason as mentioned above. Furthermore, the 15% and 85% fractile hazard curves represent the extent of the epistemic uncertainty associated with the classical approach. For the classical PSHA, the logic tree is populated with six Bayesian GMMs with different functional forms (i.e., Kea20 GMMs) that all were calibrated to the Icelandic strong-motions. Thus, such uncertainty is narrower than the uncertainty in the backbone approach, particularly for the farfield site the hazard of which is dominated by events at ~ 15-30 km distances where the Kea20 models predict very similar values, as this is the distance range in which majority of the available Icelandic data lies. The converse is true for the near-fault site where we see the variability of the ground motion predictions increasing due to the paucity of near-fault data (see Fig. 5).

Disaggregation of seismic hazard
Providing further insight into the backbone hazard results we present the disaggregation of the hazard in Fig. 9 which also allows us to identify the individual earthquake scenarios that contribute most to the hazard at these two sites, respectively. The results are shown for the annual ground motion probability of exceedance of 0.0021 (1/475 years). The disaggregation of the seismic hazard suggests the earthquake that has the dominant effects is expected to be 6.0 ≤ M w < 6.5 at a distance 14 ≤ R JB < 26 km in Reykjavik at all periods, with the contribution of larger magnitude and larger distance events becoming more pronounced with lower oscillator periods, consistent with previous disaggregation efforts ( Solnes et al. 2004), and 6.0 ≤ M w < 6.5 at a distance 0 ≤ R JB < 6 km in Selfoss, respectively.
In fact, the easternmost source zones that are capable of the largest earthquakes in Southwest Iceland appear to be sufficiently far away from the capital region to have relatively low hazard contribution. This is not the case for the town of Selfoss which is located inside the SISZ itself and in the extreme near-fault region of the causative faults. As a result, as Fig. 9 shows, earthquakes of M w ≥ 5.5 at a distances 0 ≤ R JB < 6 km dominate the hazard at all oscillator periods and again, the contribution of larger and more distant earthquakes is seen to increase progressively with lower period, signaling the increased contribution of E-SISZ to the hazard for Selfoss as larger earthquakes radiate longer period seismic waves to a greater extent than smaller events.
These results are consistent with the geometry of the earthquake source zone with respect to the sites in question and the systematic variation in maximum magnitudes along the source zone. Namely, the three seismic zones that lie in greatest proximity to Reykjavik are the C-RPOR, E-RPOR and HTJ, having maximum magnitudes set at M w 6, 6.5 and 6.5, respectively, at distances of R JB 14-26 km, with larger earthquakes taking place beyond ∼ R JB 32 km from zones further to the East. The very rapid decay of ground motion amplitudes with distance that is a well known and unique characteristic of Icelandic ground Fig. 9 2D representation of disaggregation of the seismic hazard for Reykjavik (left column) and Selfoss (right column) at PGA and PSA at T = 0.3, 1 and 2 s for associated a return period of 475 years motions, and the non-self-similar ground motion scaling of large earthquake near-fault ground motions (i.e., ground motion saturation), leads to the effects of larger earthquakes not being as far-reaching as might be expected (Kowsari et al. 2020b(Kowsari et al. , 2021. The PSHA results of this study allow us to put previous hazard results for the SISZ and the Reykjavik area into context. Solnes et al. (2004) presented a national hazard map for Iceland in which the 10% probability of exceedance of PGA in 50 years was 40% g or larger for Selfoss and in the range of 7-10% g for the capital region of Reykjavik (Solnes et al. 2002). The national annex for the earthquake resistant design of structures according to Eurocode 8 that had been in effect since 2003 showed similar values, until an update in 2010 which saw the PGA increasing to 50% g for the South Iceland Seismic Zone and to 10-20% g for the capital region (Standards Council of Iceland / Staðlaráð Íslands (SI) and Halldorsson 2010). The latter range of values is approximated from the mean hazard values at locations in Reykjavik that are furthest away and closest to the seismic zones i.e., reflecting the variability over the capital region. However, the variability around the mean values was not reported in the previous studies. In that context, we note that backbone hazard results in Fig. 8 show essentially three point-estimates at any given annual exceedance rate i.e., three hazard curves for each site and oscillator period. Namely, those corresponding to the backbone GMM mean and the GMM mean±Δ values, representing the center hazard values, and upper and lower estimates of the body of realistic hazard values. While these three estimates give a good idea of the range of hazard values, Fig. 5 clearly shows that the Δ-factor itself is associated with varying degrees of uncertainty.

Uncertainty analysis
For completeness therefore, we carry out uncertainty analysis to investigate further the effect of the variation in the scale factor with distance on the distribution of site-specific PSHA results. With this we hope to improve our understanding of the effects of relevant variables by considering some distributional form over a sample space of their possible values (Limbourg and De Rocquigny 2010). We therefore model the scale factor uncertainty by treating it as a random Gaussian variable with a period-and distance-dependent mean and standard deviation (see Fig. 5) and sample its parametric space through Monte Carlo simulations. In this way, we generate Δ-modulated backbone model predictions for each site, resulting in a distribution of hazard curves that provides us with a complete picture of the center, body and range of realistically expected hazard values. Figure 10 shows the results of the analysis, in the form of PSHA histograms for PGA and PSA at T = 0.3, 1.0 and 2.0s for ground motion exceedance associated with a return period of 475 years for Reykjavik and Selfoss, respectively. The blue and red smooth density functions that have been fit to the histograms are associated with the -Δ (i.e., lower) and + Δ (i.e., upper) hazard estimates, respectively, and the histogram's mean and standard deviations are given accordingly. The peaks of the histograms in each plot in Fig. 10 represent the two outer dots in each of the plot of Fig. 8, but here their dispersion reveals the variability of each of the point on the ±Δ hazard curves (for this return period). The results show that the observed variability of Δ around its mean is not reflected in wide histograms of the associated hazard estimates. In fact, within each oscillator period, the coefficient of variation (i.e., extent of dispersion) of the ±Δ hazard estimates is very low, and on the average lowest for the near-fault site. As a result, we are confident in the backbone hazard estimates capturing the center and body of the hazard and representing well the extent of epistemic uncertainty on the hazard.

Conclusion
The backbone approach to PSHA is fast becoming the de-facto standard practice in modern PSHA because it provides a more realistic and robust treatment of the contribution of GMM's epistemic uncertainty to PSHA results. In this study, several GMMs that have been proposed or used for PSHA in Iceland in addition to several frequently used worldwide ones and all of which satisfy the conditions on the GMM-functional form required of GMMs in PSHA have been considered as candidates for the backbone GMM for Iceland, specifically Southwest Iceland from which the majority of Icelandic earthquake strongmotion recordings exist. By applying the DIC data-driven method to rank the GMMs against the Icelandic strong-motion dataset of SISZ mainshocks, the best GMMs i.e., the backbone models have been revealed for PGA and PSA at different oscillator periods. The results show that the backbone model for Iceland is one of the recently proposed hybrid Bayesian empirical GMMs (Kowsari et al. 2020b). The backbone models at each period are associated with a scale factor, Δ , that is the model-to-model variability quantified from the selected GMMs. Thus, Δ represents the epistemic uncertainty associated with the GMMs and scaling the backbone GMMs±Δ not only captures the top-ranked GMMs but also comprehensively captures the characteristics of Icelandic strong-motions over different magnitude and distance ranges, with Δ taking the highest values in the near-and far-field regions i.e., showing minimum values at the distance ranges where most data is available.
We then carried out site-specific Monte Carlo PSHA for two representative far-field and near-fault sites using the backbone approach and compare the results with the classical approach using the logic-tree approach. The seismic source specification for the PSHA in this study is based on a simulation of a catalogue of finite-fault earthquake sources that are consistent with a physics-based model of the bookshelf strike-slip fault system of SISZ-RPOR that captures the systemic variation of the seismogenic potential of the transform zone and the seismic activity of which in fact explains the historical catalogues. The mean hazard curves from both approaches are similar but the advantage of the backbone approach is that now the mean and the ±Δ hazard curves essentially provide us with the "body" of realistically expected hazard values. This gives us a better insight into not only the hazard values themselves, but the extent of our confidence in those values as suggested by the varying width of the body at near-vs. far-field sites at different periods and different hazard levels. Previous hazard values for these near-and far-fault sites are found to be confined within the body of the backbone hazard values but now the extent of the body can be used to put the differences between past hazard estimates into proper context.
For completeness, we carried out the disaggregation of hazard in terms of magnitude and distance contributions. The modal event for the return period of 475 years for the farfield site of Reykjavik is a M w ~ 6-6.5 at R JB ~ 14-26 km, while that for the near-fault site of Selfoss is in the same magnitude range but at a distance less than R JB ~ 6 km. Other lower contributions are seen to extend over a greater distance range and towards larger magnitudes for the far-field site, while extending to both lower and higher magnitudes for the near-fault site, but always within R JB ~ 20 km. We have shown that the backbone models of this study along with the use of the effects of their Δ factors on the hazard estimates provide a clearer view of the hazard and its epistemic variation, as compared to the classical approach. We therefore conclude that the backbone approach is the ideal approach for future PSHA for Iceland, and in particular when used in conjunction the best ranked GMMs for Iceland, and when based on finite-fault earthquake catalogues that are consistent with the new physics-based fault system model of bookshelf faulting in Southwest Iceland.