Site specific probabilistic seismic hazard model for Isfahan, Iran: estimates and uncertainties

We present a seismic source characterization model for the probabilistic seismic hazard assessment (PSHA) of the Isfahan urban area, Iran. We compiled the required datasets including the earthquake catalogue and the geological and seismotectonic structure and faults systems within the study region to delineate and characterize seismic source models. We identified 7 relatively large zones that bound each region with similar seismotectonic characteristics and catalogue completeness periods. These regions were used for calculating the b-value of the Gutenberg–Richter magnitude recurrence relationship and for estimating the maximum magnitude value within each region. The earthquake recurrence parameters were then used to build a spatially varying distributed earthquake rate model using a smoothed kernel. Additionally, based on a fault database developed in this study and on a local expert’s opinion supported by tenable constraints about their slip velocity, a fault-based model is also created. We further performed sets of sensitivity analyses to find stable estimates of the ground motion intensity and to define alternative branches for both the seismogenic source and ground motion prediction models. Site amplification is considered based on a Vs30 map for Isfahan compiled within this study. The alternative source and ground motion prediction models considered in the logic tree are then implemented in the software OpenQuake to generate hazard maps and uniform hazard spectra for return periods of interest. Finally, we provide a detailed comparison of these PSHA outcomes with both those presented in the 2014 Earthquake Model of Middle East (EMME14) and with the national seismic design spectrum to further discuss the discrepancies between hazard estimates from site-specific and regional PSHA studies.


Introduction
Probabilistic seismic hazard analysis, PSHA (Cornell 1968;McGuire 2004) is the standard procedure for estimating the ground shaking hazard at a site (McGuire 2004). The nature of the phenomenon dictates the application of probabilistic approaches to blend information about seismic sources, earthquake magnitudes, seismic wave traveling paths, the ground motion intensity and site response with the information about the historical earthquakes and seismotectonic settings. In other words, PSHA considers "all" potential earthquake scenarios that can affect a given site combined with the distribution of ground motion intensity (commonly defined through ground motion prediction equations, GMPEs) to be expected at the site for all such earthquake scenarios. The objective is the estimation of the likelihood of observing at a site in a given time frame ground motion intensity measures, such as peak ground acceleration (PGA) and spectral acceleration (SA), of different severity. Performing seismic hazard analysis is essential for developing resilient cities in seismic prone areas because it provides the fundamental requisite for the design of new safe buildings and for the assessment of the existing ones.
Iran as one of the most seismically active regions in the Middle East, has a long history of hazard studies starting from early sixties. Since 1962 to date, different versions of the seismic design guidelines and codes have been adopted country-wide. The latest version of the seismic design code (i.e., fourth edition of Standard 2800) indicates the level of seismicity for different regions of the country based on a seismic zonation map (ICSRDB 2014). Besides Standard 2800, 2 regional hazard studies of GSHAP 1 (Giardini 1999) and EMME14 2 project (Şeşetyan et al. 2018) also include Iranian territory. Several studies conducted hazard studies specifically for Isfahan urban area including identification of active faults and performing hazard analysis (Beygi et al. 2016;Safaei et al. 2012;Tabaei et al. 2016;Tajmir Riahi et al. 2014).
In this paper, we describe a probabilistic seismic hazard model developed for Isfahan following the concepts and methodologies set by Field et al. (2003) implemented in the OpenQuake engine (Pagani et al. 2014). Our model is based on the latest information in openly accessible datasets and scientific literature and we supplemented it with local experts judgement. The model consists of a combination of distributed seismicity areas and active faults within the region. The former is underpinned by the direct analysis of earthquake catalogues, while the latter is based on the spatial delineation and activity rates of faults derived by geological and geodetic observations together with meaningful constraints on the faults slip rates. We intend to combine these 2 source models into 1 model to produce a robust estimate of the seismic hazard in the urban area of Isfahan. We also discuss our model development assumptions and the corresponding implementation, including the key aspects of the ground motion models, the characterization of model uncertainties, the logic-tree weights, the sensitivity analysis and results. Finally, we critically compare the outcomes of our model both with those obtained in the regional EMME14 model and with those proposed in the national seismic design code of Iran.

An overview of the study area
Isfahan is located in central Iran with the approximate coordinates of 32°38′41″N and 51°40′03″E. Isfahan has an urban area of about 551 km 2 and an average altitude of 1574 m. Figure 1a shows its location on the map of Iran while Fig. 1b shows the extents of its urban area. The Isfahan municipal is the third most populous city of Iran with a population of more than 2 millions according to the 2016 Census results published by the Statistical Center of Iran. 3 In this section we provide a brief summary of the geological and geotechnical properties of the study area.

Geological characteristics of the study area
The Iranian plateau is located between the Arabian Plate in the southeast, and the Eurasian Plate in the northwest making part of the Alpine-Himalayan orogenic belt (Berberian and King 1981). This belt is seismically active with a long history of large devastating earthquakes. Figure 1a shows an overview of the tectonic map of Iran. Our study area is mostly located in the central metamorphic Sanandaj-Sirjan zone and also partially in the high Zagros belt and Zagros simply folded belt, Urumieh-Dokhtar magmatic assemblage and Central Iran. In a more general geological classification, considering similarities of sedimentary facies and structures, the Sanandaj-Sirjan zone and the Urumieh-Dokhtar magmatic assemblage can be considered as part of the Central Iran zone (Berberian and King 1981). A brief description about the activity and the geological features for these four mega zones is provided in the following paragraphs.  Safaei et al. (2012). According to this map, the study area is mostly located in the North Sanandaj-Sirjan zone, in the South Sanadaj-Sirjan and partially in the high Zagros belt, Zagros simply folded belt, Urumieh-Dokhtar magmatic assemblage and Central Iran. Note: The red circle indicates 200 km distance from Isfahan; b Geotechnical setting of Isfahan showing the boundaries of the soil types, topography, and the location and type of the available geotechnical test data.
Note: Panel b also shows the extents of 15 municipal districts of Isfahan Sanandaj-Sirjan-is a quasi-rigid geological zone which is the result of a polyphase deformation that happened in the greenschist and amphibolite metamorphic facies. This deformation is most strongly tectonized in the Zagros orogeny (Mohajjel and Fergusson 2000;Şengör 1990). Both the southern and the northern parts of the Sanandaj-Sirjan zone are made of metamorphic rocks. While the former is most probably metamorphosed in the Middle to Late Triassic, the latter is deformed during the Late Cretaceous and therefore it contains many intrusive felsic rocks (Eftekharnejad 1981;Ghasemi and Talbot 2006;Safaei 2009). Due to the occurrence of several metamorphic phases in this zone, its properties have become similar to those of stable continental regions (SCRs). Even though in this zone the seismic activity has been lower than that of other more active regions of Iran, reports about the occurrence of several historical events indicate that it can be considered a seismically active zone (Ambraseys and Melville 1982). The low deformations in SCR, however, cannot imply that large earthquakes are not possible. Along these lines, Craig et al. (2016) argued that the tectonically stable continental lithosphere can store elastic strain on long time scales, the release of which may be triggered by rapid, local transient stress changes caused by surface mass redistribution, resulting in the occurrence of intermittent intraplate earthquakes. In fact, based on the historical records, earthquakes in this region are expected to be large albeit with long interarrival times. This consideration further highlights the fact that the failure to record large earthquakes along the identified faults in this region cannot indicate that these faults are inactive. As such, as suggested by Berberian (2005), more local investigations to identify the active faults in such regions can significantly help to understand their tectonic behavior.
High Zagros belt and Zagros simply folded belt-are located between the Main Zagros Thrust and the Main Recent Fault while the Zagros simply folded belt is located between the Main Recent Fault and Mountain Front Fault. Most of the seismicity in Zagros is due to thrust faulting (Hassanpour et al. 2018;Jackson et al. 1995).This is one of the most seismically active fold-and-thrust belts globally. Mapped surface faulting associated with earthquakes is very rare, with most seismicity occurring on blind reverse faults buried beneath a sedimentary cover. Therefore, the distribution of earthquakes provides the most useful data about the location of active faults (Karasözen et al. 2019). The 2017 Sarpol Zahab M7.3 earthquake is one of the largest recent earthquakes occurred in this zone (Feng et al. 2018).
Urumieh-Dokhtar magmatic assemblage-is approximately 100-150 km wide and it is located about 150-200 km from the Zagros zone. It contains various intrusive and extrusive rocks along the active margin of the Iranian plates (Alavi 1994;Berberian and King 1981). There are many active faults, such as Kashan and Dehshir faults (Sect. 3.1), in this zone and evidence suggests that the magma ascended mainly in the faults with various direction (Chavideh and Tabatabaei Manesh 2021). The occurrence of several historical earthquakes and the recording of instrumental earthquakes in this zone indicate that it is an active geologic zone (Ambraseys and Melville 1982).
Central Iran-zone covers a large part of the center of Iran. Among the most important events occurred here are the 1962 M7.2 Buin Zahra, the 1968 M7.1 Dasht-e-Bayaz, the 1978 M7.4 Tabas and the most recent 2003 M6.6 Bam earthquakes. The M7.4 Tabas earthquake occurred centuries after the last major earthquake experienced in the region. This indicates that there is a possibility of a devastating and large earthquake in areas of central Iran even if a devastating earthquake has not been reported (Berberian 1979). In Bam, no historical earthquake had been reported before the 2003 event.

Geotechnical characteristics of Isfahan urban area
The geotechnical setting of Isfahan is investigated in several studies (Abdolahi et al. 2014;Ajalloeian and Jami 2010;Farashi and Ajalloeian 2012;Mohammadi et al. 2020). Herein, together with the past studies, we used the results of 369 geotechnical investigations based on standard penetration tests, SPT (303), stratigraphic logs (13) and Vs 30 measurements (52) to identify the geotechnical setting of the study area. Figure 1b shows the location and type of the available geotechnical test data. These resources show that a major part of Isfahan is located on Quaternary alluvial sediments while small areas in the southern border and north western (i.e., areas closer to the hills) are built on the Jurassic rock units. The bedrock consists of alternating layers of "shale" and "sandstone". The alluvial sediments are composed of 2 distinct parts of "flood plain deposits" and "fluvial deposits" and in depth they could be distinguished by three layers of surface fine-grained, middle-coarse and lower fine-grained sediments. The predominant texture of these sediments is clay and mud. The alluvial sediments start with a small thickness in the south of the Zayandehrud River (see Fig. 1b) and the thickness increases by moving towards the north. The internal friction angle and cohesion of the alluvial sediments range between 30-40° and 0.1-0.2 kg/ cm 2 , respectively, with an average elasticity modulus of 250-850 kg/cm 2 .
In this study we use Vs 30 as a proxy to account for the effects of the local soil conditions on the ground motion intensity at the surface. As such, to estimate Vs 30 in the study area, we used an hybrid approach based on 2 resources: (i) available geotechnical test data, and (ii) a statistical method based on the topographical slope proposed by Wald and Allen (2007). Given that the test data do not cover the entire area, we give priority to the tests, where available, while the slope-based data are only used to fill in the spatial gaps. More details about the procedure used in developing the Vs 30 map are described in the supplementary material of this manuscript. It is worth noting that in the current study we have not investigated the sensitivity of our results to the uncertainty in the site condition, which is expressed by the uncertainty in the Vs 30 map. Such a sensitivity could be investigated in a separate study. The Standard 2800 (ICSRDB 2014) Iranian national design code specifies 4 different soil types based on Vs 30 values, namely Type I: Vs 30 > 750 m/s; Type II: Vs 30 = 375-750 m/s; Type III: Vs 30 = 175-375 m/s; and Type IV: Vs 30 < 175 m/s. Figure 1b shows the spatial distribution of the soil types based on this soil classification. According to this map, soil Type III is predominant (about 68% of the total area of interest) in the northern and central regions, while the southern part mainly consists of soil Type II (corresponding to about 32% of the area).

Active faults
To provide an input for active fault-based PSHA, we created a new dataset of active faults considering a region within 200 km of Isfahan. However, the focus in compiling our database has been on a region of about 100 km around the site. Analysis of satellite data (Landsat ETM + and Aster) and images of geographic browser (Google Earth and Sasplanet) followed by field studies enabled the recognition of many previously unmapped faults. As a result, we generated and further processed a raw database including the characteristics of 125 faults. It should be noted that in absence of extensive trenching and excavation data, we could not precisely determine all the features of these faults. However, based on remote sensing, aerial geophysical data processing, field surveys in previous studies (Beygi et al. 2016;Safaei 2009;Safaei et al. 2012) as well as experts judgments, we parameterized these faults to the extents possible.
Given the dextral transpression tectonics of the region, we expect to develop and extend the faults with different mechanisms such as riedel shears faults, antithetic faults, synthetic faults and tension faults, which are well identified in the region. Other structures expected in transpression environments are the positive flower structure or upward-diverging faults (Fossen 2010). In this structure, when a main basement fault reaches near the ground, it turns into several parallel faults with the same mechanism. As a result, several parallel fault traces can be seen on the surface, while at depth, they merge into single main fault. In our study area, a flower structure is well observed and a large number of faults with northwest-southeast strike direction with this structure is identified. Based on these premises, we grouped the faults into 3 classes that reflect their level of importance and activity. In this process, only 1 representative fault is selected for a flower structure and the faults with short lengths or those with no associated earthquakes are classified as belonging to Classes 2 or 3. Note that in our PSHA study, we used only the most important class 1 faults, still, we emphasize that all major faults capable of generating large earthquakes in the region are kept in our final fault database. Figure 2 shows the extent of the study area. Figure 2b shows the traces of the 23 faults classified Class 1 while Table 1 lists more information about each fault. This Figure Table   2C of Wells and Coppersmith (1994) ǂ ǂ Strike is the average value estimated along the fault length while the dip and rake are adapted based on the global fault database proposed by Kaklamanos et al. (2011)  describing the geometry, kinematics, slip rate, and the associated epistemic uncertainties. Note that, from a geological point of view, the most reliable parameters of this database are the fault trace and mechanisms while the values of other parameters such as slip rate, dip and rake are inferred and constrained based on our best estimates and global data. The slip rates were first estimated based on expert judgment of the geodetic and seismic data, as well as considerations of geomorphic expression of similar but better studied faults in the active regions of Iran (Alipoor et al. 2012;Jamali and Hessami 2011;Khodaverdian et al. 2015;Walpersdorf et al. 2014) and for stable continental regions worldwide (Allen 2020;Calais et al. 2016). This provided a preliminary understanding of the relative movement of the faults with respect to each other. In a second step, the central values of the slip rates (i.e., SRM) were further calibrated based on the observed earthquake rates and the magnitude frequency distribution obtained for the seismogenic super zones of the study area (see Sect. 4.1). In addition, to account for the uncertainty on slip rates, an upper and lower bound value, SRL and SRU, were also introduced. We inferred SRU/SRL values based on 2 assumptions: (i) the regional variability of the slip rate estimates observed in the Iranian active faults database such as the ones given in EMME14 (Danciu et al. 2017); and, (ii) the error range of the mean slip rates provided by the local geologist. On these premises, for super zones (such as Zagros) with better known history of seismic activities where earthquake data for M w > 6.0 was available, SRU/SRL ratios as low as 2.0 was adopted. For less known regions or where no M w < 6.0 was reported, SRU/SRL ratios of 3.0 and 4.0 were considered. We emphasize that, for a more precise modeling of the faults, an extensive study to better estimate the slip velocity measurements of the new faults introduced within this study is strongly recommended in future updates of this model.
For each fault, a distribution of maximum magnitude (M max ) is calculated based on its identified features. M max estimation along a fault segment has been a topic of discussion for many years (Mignan et al. 2015;Zöller et al. 2013). Herein, we consider multiple approaches to estimate the values of M max : (i) based on the maximum observed magnitude along the fault ( M obs max ), (ii) as a function of the seismic moment, M M 0 max (Hanks and Kanamori 1979); and (iii) based on faults capacity using the empirical relationships, M cap max (Leonard 2010;Wells and Coppersmith 1994). These relationships are functions of the faults physical parameters, such as surface fault length, subsurface fault length, down-dip fault width, fault area and average displacement per event. Table 1 lists a summary of the calculations based on different methods. Based on methods (ii) and (iii) and by variation of the fault parameters (e.g., dip ± 15° and the corresponding width) as well as the distribution of M max reported in the above studies, we sampled M max for M

Earthquake catalogue
A PSHA study commonly starts with assembling a complete (in space and time) earthquake catalogue of events characterized by a homogeneous magnitude definition. The catalogue, which is required for proper definition of the past earthquake occurrences (and for forecasting rates and location of future ones), is the basis for the development of area seismic source recurrence characteristics and for performing smoothed seismicity analyses. Earthquake catalogues typically consist of a historical part, which summarizes historical reports on observations of earthquakes occurred in the distant past, and an instrumental part, which relates to the recent seismicity detected by seismic networks. In the following, we provide a description of the catalogue resources and of the earthquake magnitude harmonization procedure. Table 2 lists the sources of the catalogues as well as the covered period and magnitude ranges. Figure 3 shows the spatial distribution of the events and their estimated moment magnitude (M w ) values.  Ambraseys and Melville (1982) collected information about historical earthquakes in Iran, originally based on chronicles of historical manuscripts. Herein, we identified 17 historical events that are reported in terms of macroseismic intensity (MMI) within the study region. The instrumental catalogue obtained from 2 sources, ISC and IIEES, spans an area between 29.0-36.0°N and 48.0-55.0°E. ISC is a more complete database than IIEES, it provides more precise coordinates of the epicenters but is only available from 2006 to date. As such, we compiled the instrumental catalogue using the IIEES database in the interval of 1909-2005 and using the ISC database for the interval of 2006-2020. The magnitude values in the entire catalogue were then harmonized to moment magnitude using the conversion relationships in Table 3. For the historical events, we considered a 2-step process that converts the MMI to surface magnitude (M s ) through the conversion equations proposed by Papadopoulos et al. (2014), and from M s to moment magnitude (M w ) using the conversion relations of Scordilis (2006). It should be noted that this 2-step approach may introduce additional uncertainty to the model. In fact, alternatively one could directly convert MMI to M w (Álvarez Rubio et al. 2011;Atkinson et al. 2014a). IIEES reports magnitudes in terms of body wave magnitude, m b , surface magnitude, M s , and local magnitude, M l , while ISC is based on Nuttli magnitude (M N ). All different magnitude types are converted to M w using the equation of Scordilis (2006) for m b and M s , that of Zare et al. (2014) for M l , and that of Sonley and Atkinson (2005) for M N .
The instrumental catalogue includes all events and not only mainshocks. Assuming a Poissonian recurrence of the earthquakes in terms of location and magnitude, typically only mainshocks are used to assess the recurrence parameters in a time-independent seismic hazard assessment (Frankel 1995). Therefore, conforming to this assumption, herein we removed the clustered events from the instrumental catalogue using the Gardner and Knopoff (1974) declustering algorithm based on the standard parameter values of that method. The declustered instrumental catalogue is left with 3,439 out of the 11,416 events spanning the magnitude range of M w = 2.0-7.3 contained in the original catalog. The choice of different methods, such as the declustering algorithm or the homogenization of the catalogue, may affect the final PSHA results (Beauval et al. 2020). However, the impact is expected not to be critical in this regional study and, thus, no sensitivity analysis on the matter was done here. We assume that the adopted approaches represent a central path. Albeit, one could compile multiple catalogues from different resources to represent the inherent variability of the regional seismicity, this is not doable nor practical as no alternative seismic networks are available in the region. Thus, the earthquake catalogue in this investigation represents the most updated and relevant for the region of interest.  Sonley and Atkinson (2005)

Seismic source characterization
As alluded to earlier, we generated a source model consisting of a combination of 2 components based on distributed seismicity and identified active faults. The former stands on the analysis of the harmonized earthquake catalogue, while the latter uses the database of active faults and their slip rates compiled in the current study. The distributed seismicity model is well calibrated to account for spatial variability of the earthquake process without association with the fault systems, whereas the active faults constrain the long recurrence interval earthquakes of localized seismicity around the known major tectonic structures. As such, we intend to compensate their limitations through a reasonable combination of both.
In the following sub-sections, we briefly discuss the main characteristics of the proposed source model.

Distributed seismicity model
Seismic source zones (SSZs) are delineated on the basis of geological, geophysical, and/ or seismological similarities (Reiter 1992). SSZs are basically geographic polygons of assumed uniform seismicity. In addition, the maximum magnitude, depth, style of faulting, strike, dip and rake are assumed spatially invariable within a zone. Herein, we follow the prescriptions of Vilanova et al. (2014) and Danciu et al. (2016) for the definition of the zones. Our zonation is mainly based on (i) the analysis of the earthquake catalogue in terms of completeness intervals, seismic activity rates and seismogenic depths; and, (ii) the surface geometry of the faults and the geological settings as well as the style of faulting. Based on this concept, we discretized the study area into 7 independent source zones, Z1-Z7. In this zonation process, we have assured the compatibility with the EMME14 seismogenic source model for this region (Danciu et al. 2017;Şeşetyan et al. 2018) and the geological zonation studies of Iran (Arian 2015; Nogole-Sadat and Almasian 1993; Oveisi et al. 2019). Figure 3 shows the extent of the zones and the spatial distribution of the events within each zone. Z1 and Z3 lie within a tectonic environment that is adjacent to major plate boundaries. However, they have not generally undergone any recent large tectonic deformations or large earthquakes. It should be noted that characterizing these regions only based on the observations and background seismicity may result in an underestimation of the seismic hazard. We compensate this limitation by thoughtful considerations on the choice of M max and by relying on the fault source model for generation of long recurrence interval large magnitude events. Z3 is basically located in Sanandaj-Sirjan zone (see Fig. 2b) and, as described earlier, despite the few events in this region with no events of M w ≥ 5.0, it cannot be considered inactive.
Within the framework of EMME14, the local geologist and experts have proposed Z1 and Z3 to be classified as stable shallow crust (i.e., SCR) similar to the Central and Eastern North America, CENA Johnston et al. 1994) or Australia (Crone et al. 1997). This region consists of relatively old crust, with very low deformation rate, resulting in episodic earthquakes of low magnitude and very low stress drop. Furthermore, this domain is also characterized by slow anelastic attenuation, as seen from the observed ground motions from the Iranian strong motion database. In addition, the number of active faults identified in Z3 and the occurrence of small earthquakes caused by the movement of these faults (Safaei 2013) imply that the region might be considered also as an active shallow crustal (ASC), despite the few historical large events. However, it should be noted that this assumption may not be shared by other local studies in Iran. To address such a difference in interpretation, we considered Z1 and Z3 as hybrid (ASC or SCR) regions in our logic tree. The impact will be reflected in the GMPE selection and the relative weights, which will be discussed in the following sections. Although Z1 is outside of the 200 km distance from the center of Isfahan and the events in this zone have a negligible impact on the urban area of Isfahan, the contribution to the hazard from Z1 is also accounted for here for the sake of completeness. The earthquake activity in Zagros (i.e., Zone Z2) is significantly larger with 62 events with M w ≥ 5.0. Zones Z4, Z5 and Z6 are delineated in order to lie within Urumieh-Dokhtar, Alborz, and Central-Iran geological regions, respectively. Note that Z2, Z4, Z5 and Z6 are classified as ASC regions.
In characterizing zones Z1-Z6, we have considered all the events with hypocentral depth shallower than 50 km. Figure 3b, which shows the spatial distribution of events depths, shows that most of the deep events (with hypocentral depth ≥ 50 km) are concentrated beneath the Zagros fold. According to the geological studies, this region consists of mainly carbonate and marly facies, i.e., early Cretaceous (Nogole-Sadat and Almasian 1993). As such, zone Z7, which lies partially within Z2, is also considered to represent the deep seismicity region of Zagros.
For each 1 of the 7 zones we derived a doubly truncated Gutenberg-Richter (GR) relation from the declustered catalogue following a completeness analysis based on the method of Stepp (1972). For each zone, the a-value (a GR ) and b-value (b GR ) of the GR relation were then calculated using the maximum likelihood estimation method that considers unequal completeness intervals for different magnitude ranges using the Weichert (1980) algorithm. To achieve a good fit with the data, in the completeness analysis of the catalog we have considered multiple magnitude thresholds (i.e., M GR lower ) and magnitude bins. These parameters are iteratively adjusted to provide the best solution of the GR relation and, as such, they vary between zones.
Note that when performing the hazard calculations, a common minimum magnitude threshold of M min = 4.0 has been considered for all zones, assuming that this is the lowest magnitude of an event that is capable of causing damages to the structures in the urban area of Isfahan. With a seismic risk analysis perspective, Bommer and Crowley (2017) defined M min as the lower value of the earthquake magnitudes that, if reduced, does not result in higher risk estimates even though it may result in higher seismic hazard estimates.
The estimation of the maximum magnitude in PSHA is a complex, and often controversial, aspect that should be guided by information from geology and the seismotectonic of a seismic source. Herein, for definition of M max for each zone, we considered 3 approaches including: (i) the maximum observed event within the zone (M obs ) plus 0.5 unit ( M obs max ) ; (ii) the non-parametric Gaussian method of Kijko (2004) Table 1) based on empirical relationships. Based on these 3 estimates, the maximum of the 3 values is selected as the M max of the zone. Regardless, no M max values lower than 6.5 were selected for any zone (Wheeler 2011). However, to account for the uncertainty on this value, we considered different models with M max ± 0.3 in the logic tree adopted in the PSHA. Further discussions on the structure of the logic tree are provided in Sect. 6. Figure 4 shows the GR fit for each zone in comparison with the observed data. Table 4 provides a summary of the results of the zonation analysis and also compares the corresponding GR parameters and M max values that were adopted in EMME14. Note that the delineation of the zones in this study is different than that in EMME14 and, thus, we could only approximately extract the GR parameters from EMME14 using the definition of  We analyzed the depth of the events from the catalogue within each zone to characterize the hypocentral depth distribution in hazard calculations. Figure 5 shows the normalized histogram (i.e., frequencies of being in the interval) of the event depth for each zone. We considered 5 depth intervals for zones Z1-Z6 that range between 0 and 50 km and 4 depth intervals from 50 to 140 km for Z7. These frequencies are then utilized in the source modeling as probabilities that the depth of the next earthquake is in the corresponding ranges. It should be noted that before making a statistic of the data, we have eliminated the depth values of 0, 5, 10 and 33 km from the database, assuming that when no information was available, such round numbers were manually inserted in the database. With reference to Note that we found more than 160 events deeper than 50 km in our catalogue most of which occurred in the Zagros region (see Fig. 3b and the area in Z7). This is in contrast with the findings of several past studies (Kalaneh and Agh-Atabai 2016; Karasozen et al. 2018;Talebian and Jackson 2004). This may imply that there may be a systematic error in the reported hypocentral depth of deep events of our catalogue. This topic, nonetheless, is debatable as several geological evidences are reported in previous studies (McQuarrie 2004;Sepehr and Cosgrove 2004) addressing the possibility of deep seismicity in the region, due to the differences in time and the manner of collision of plates in different parts of the Arabian and Iranian plates. Such phenomenon has clearly caused an increase in the thickness of the Iranian crust in the central parts of Zagros and thus provides an evidence for deep seismicity in zone Z7. In the delineation of Z7 we used satellite images to recognize the borders of the outcrops, which are mainly carbonate and marly facies of early Cretaceous, and to distinguish between the deep and shallow parts of the region. Finally, as stated earlier, based on the data in our possession, we decided to keep deep zones in our study, although more investigation on this subject is needed. We emphasize that, given the uncertainty of such an assumption, it would be worth performing a sensitivity analysis on the depth of events in zone Z7 if the contribution to the hazard estimates from this zone were significant. However, given the large distance of Z7 from our study area, the hazard contribution from this zone are indeed limited (See Figure S20 in the supplementary document) and, thus, further investigations were not carried out.
To account for the spatial variability of the seismicity, as described in Poggi et al. (2020), we use a smoothed seismicity approach (Frankel 1995) to redistribute the annual occurrence rates within each polygon (i.e., Z1-Z7) considering the spatial density of earthquakes using an isotropic seismicity smoothing kernel. This approach enables the use of zone-specific b GR and seismic modeling parameters such as rupture mechanism, depth distribution and M max . We considered a set of spatially distributed grid points with 0.1° × 0.1° spacing covering the entire area.
In addition, to account for the modeling uncertainty, we used three different correlation distance values, C d , (Frankel 1995) equal to 10, 40 and 70 km. As an illustrative example, Fig. 6 shows the cumulative annual rates (M w > 0) within the zones using C d = 10 and 70 km. It should be noted that for each grid point of the distributed seismicity model we have adopted the dip and strike angle from the closest fault while the depth aleatory uncertainty is defined according to the depth distribution of each zone, as shown in Fig. 5.

Faults source modeling
In this section we describe the fault source modeling approach including the style of faulting adopted and magnitude-frequency distribution (MFD) derivation. We herein follow the modeling approach referred to as "simple fault" in the OpenQuake engine (Pagani et al. 2014), in which faults are represented by the surface projection of the fault traces with a constant dip angle. In absence of detailed fault information, this is an easy way to account for complexity of the fault geometry with no considerations about the changes in values of the dip angle along the fault strike. We assume that each fault is bounded within the seismogenic thickness, which we herein considered in the range of 5-30 km. The spatial extent of the ruptures along the fault is defined based on the magnitude scaling relation of Wells and Coppersmith (1994) as well as on the rupture aspect ratio.
To characterize the seismic activity of the faults, we use the truncated exponential model of Youngs and Coppersmith (1986) in agreement with the distributed seismicity model. Basically, the earthquake occurrences are defined based on the slip rates and the fault area. Note that the b GR value for each fault is derived from the statistical estimates obtained from the entire zonation study (see Table 4). In the MFD calculations for all faults we assumed a crust shear modulus of μ = 3.0 × 10 11 dyne/cm 2 (Dziewonski and Anderson 1981). In our modeling we have assumed that the slip rate is constant across the fault area and that there is no creeping. It should be noted that we investigated the sensitivity on the uncertainty of the faults slip rates (via SRL, SRM and SRU listed in Table 1) in our logic tree. The sensitivity of the PSHA results to geometric parameters of the faults, such as fault depth, dip and rake angles, that are less critical in a regional study was not carried on here. Instead, to keep the model manageable, we have considered the best estimates of these geometric parameters. However, we would strongly recommend a more detailed sensitivity analyses on the values of such parameters in a site-specific hazard analysis for important structures in this region.

Final seismogenic source model
As shown in previous studies (Hofmann 1996;Wesnousky et al. 1983) and implemented in other urban (Ebrahimian et al. 2019) and regional (Danciu et al. 2017;Poggi et al. 2020) hazard analysis studies, modeling a unique MFD relationship for the faults for the entire magnitude ranges (small, moderate and large) is not recommended. Such MFD modeling may result in an overestimation of the annual rates of small to moderate earthquakes compared to those observed. On the other hand, using the MFD modeling only based on distributed seismicity may under predict the annual rates of larger events. Moreover, given that in the current study the location of the major faults is relatively well constrained, it is more realistic to assume that the larger events occur along the identified faults. A common practice to alleviate this complexity is to use fault modeling only for the larger magnitude earthquakes and to use a distributed seismicity model for small to moderate events. Such a modeling approach will not ignore the possible occurrence of events off the faults.
Herein, we adopted fault source modeling for magnitudes larger than a threshold value, M w > M thre . In addition, for each fault we considered an area within the fault surface projection (in-buffer) to distinguish it from the area outside (off-buffer). A map showing the inbuffer and off-buffer grid points is presented in Figure S3 of the supplementary document. Parallel to the fault source model, we considered a distributed seismicity model split into 2 sub-models containing the in-buffer and off-buffer grid points. The in-buffer grid points of the distributed seismicity model cover the magnitude range of M min ≤ M w ≤ M thre , while the off-buffer grid points are defined for the magnitude range of M min ≤ M w ≤ M max , where M max is the maximum magnitude of the zone where the grid point lies in (see Table 4). Note that to investigate the sensitivity of the results to the choice of M thre , we considered the three values of M thre = 6.0 and 6.5 in our logic tree. More details about the logic tree and the weights are provided in the following sections.

Ground motion prediction equations and logic tree
GMPEs are one of the main factors that influence the earthquake hazard and risk assessment results (Crowley et al. 2005). In parallel with the improvements in strong motion networks and the availability of ground-motion data, the number of GMPEs has significantly increased (Douglas and Edwards 2016). The large database of available GMPEs makes the procedure of GMPE selection a scientific challenge. In this process, several aspects of the GMPEs including the tectonic regime, the site effects modeling, regional or global database utilized for the development as well as distance and magnitude applicability of the model should be accounted for to avoid arbitrary choices (Danciu et al. 2016).
Before proceeding with the ground motion selection of this study, three seismotectonic regions-active shallow crust, stable shallow crust and deep seismicity-are considered for the study area (i.e., within 200 km around the Isfahan city). Noteworthy is that this classification is subjective and aims at capturing the ground motion variability due to source and wave traveling path. From the seismic source point of view, the earthquakes sources and mechanisms are generally located in the shallow crust and more distant deep earthquakes are located beneath Zagros region (see Fig. 3).
In the recent years, 2 frameworks have been commonly used to handle the ground motion epistemic uncertainties through logic tree with multiple independent models and logic trees built upon a single model (so called backbone approach). For this reason, several methods have been used in the literature to select and evaluate the performance of GMPEs in PSHA using observed ground motions from the region of interest. These methods include residual analysis methods (Scasserra et al. 2009), likelihood, LH, and log-likelihood, LLH (Mousavi et al. 2012) methods, and Euclidean distance based ranking, EDR, methods (Kale and Akkar 2013), among others.
According to Danciu et al. (2016), a regional hazard assessment requires the use of multiple GMPEs that can both describe the local characteristics of the ground motion database of the region and at the same time they should account for the completeness of the dataset at a large scale. Multiple selection models are subject to biases of experts elicitation (Delavaud et al. 2012). This is because different models are often derived from the same strong motion datasets and thus the independency of the models is commonly not achieved. On the other hand, the backbone models are also built upon subjective choice of a single model, that is then scaled up and down to describe the epistemic range of the median ground motion (Atkinson et al. 2014b;Douglas 2018;Weatherill et al. 2020). In regions of low seismicity (such as a large area surrounding Isfahan) with limited strong motion recordings, it is obviously difficult to follow the state-of-practice selection and testing of the GMPEs, and therefore herein we are left with an expert elicitation procedure to define a multiple GMPEs logic tree as follows.
Firstly, a literature overview was used to identify the GMPEs used in various projects and models. The starting point was to revisit the EMME14 GMPEs selected for shallow crust earthquakes: i.e. ASB14 (Akkar et al. 2014), CY08 (Chiou and Youngs 2008), AC10 (Akkar and Cagnan 2010) and Zetal06 (Zhao et al. 2006). Since the completion of EMME14, many of these GMPEs are superseded by more recent publications, and we aimed at selecting recent models due to their updated datasets and functional forms. Secondly, a literature review indicated recent studies that provide insights of the ground motion performance with global, regional or local strong motion data, for similar tectonic context. According to Fallah Tafti et al. (2017) and based on the EDR method, CB14 (Campbell and Bozorgnia 2014), BSA14 (Boore et al. 2014), CY14 (Chiou and Youngs 2014), ASK14 (Abrahamson et al. 2014) and Zetal06 (Zhao et al. 2006) were ranked as the most appropriate GMPEs for the seismicity in Iran. Zafarani and Mousavi (2014) investigated the best fitting GMPEs for Iran based on two methods of LH and LLH. Their study supports the use of AS08 (Abrahamson and Silva 2008), CY08 and Getal09 (Ghasemi et al. 2009) based on the former method; whilst Getal09, AS08 and CY08 were ranked as the best GMPEs according to the latter one. More recently, Farajpour et al. (2020) concluded that the NGA-West GMPEs of BSA14 and ASK14 perform better for Iran. This study also showed that, among the local GMPEs, those proposed by Darzi et al. (2019) and Farajpour et al. (2019) fit well with the recorded data at short periods, e.g., between PGA and Sa(0.2 s), while for longer periods the models of Zafarani et al. (2018), Sedaghati and Pezeshk (2017) and Ketal15 (Kale et al. 2015) are preferable.
Following these findings and suggestions while also seeking to some extent a consensus with the aforementioned local studies for Iran, we selected the following four GMPEs: Ketal15-to describe the Iranian strong motion datasets; ASB14-as the regional model that includes a small number of the Iranian data, but it captures the strong motion of similar tectonic regimes in Turkey and Southern Europe; CY14-which is probably the most robust GMPE of the NGA West database; and, CF15- (Cauzzi et al. 2015) to capture the global variability of the median ground motion, due to large magnitude events. Given the relatively fewer GMPEs available in the literature for SCR, we have a limited number of choices. The evaluations performed for CENA (Goulet et al. 2018) (2007) as top performers. However, the first 2 models are derived only for a unique Vs 30 value of about 2500 m/s, a soil condition that is not compatible with the requirements of this analysis (i.e., the use of a range of Vs 30 values). Thus, for SCR we have selected 2 models of AB06 (Atkinson and Boore 2007) and AB11 (Atkinson and Boore 2011) and we use the ASC models to complement these 2 GMPEs, in line with Delavaud et al. (2012). This combination of models is referred to as a "hybrid logic tree" where we use multiple GMPEs to represent the epistemic uncertainties of median ground motion for SCR. Note that given that considering four ASC and 2 SCR GMPEs with equal weights in the logic tree for the hybrid zones implies more weight for ASC than SCR, i.e., about 66% versus 34%, respectively. Finally, we used the 2 GMPEs of Mont17 (Montalva et al. 2017) and BChyd15 (Gregor and Addo 2015) for the deep seismicity zone (i.e., Z7). Figure 7 compares the median estimates of the 4 GMPEs selected for ASC and SCR in trellis plots. This Figure shows the spectra in a period range of 0.01-4.0 s for 16 different relevant scenarios to our study area for Mw of 4.5, 5.5, 6.5 and 7.5 and R jb distance of 10, 20, 50 and 100 km. This figure is meant to identify the center and range of median ground motion values for different magnitude and Joyner-Boore distance pairs (i.e., M w and R jb ) Fig. 7 Trellis plots for the median spectral quantities of the six GMPEs of Ketal15, ASB14, CY14, CF15 (for ASC), AB06 and AB11 (for SCR) for different scenarios corresponding to M w = 4.5, 5.5, 6.5 and 7.5 with a hypocentral depth of 10 km. The response spectra are computed for R jb = 10, 20, 50 and 100 km for a hypothetical rock site with a Vs 30 = 760 m/s subjected to a rupture on a vertical fault with strike slip mechanism generated using the selected GMPEs for a Vs 30 = 760 m/s. Each plot represents an earthquake scenario described by a vertical fault, strike slip fault mechanism, with a hypocentral depth of 10 km. Note that slight inconsistencies among curves may arise from the different distance metrics used by each GMPE and also by the different definition of the spectral quantiles plotted (i.e., geometric mean of two horizontal components, median horizontal, average horizontal, median of the geometric mean of the rotated records rotated over the usable range of oscillator periods, i.e., GMRotI50).
With reference to Fig. 7, the median values of CF15 and AB11 are in the upper range for the M w -R jb scenarios most relevant for this PSHA, i.e., M5.5 and M6.5 at R jb at about 20-50 km from the site. The CF15 is based on a global dataset (see Table 1 of that paper) with large magnitude events (i.e., M w = 8.0-M w = 8.5) and this characteristic may partially be responsible for large amplitude ground motion estimates. The median estimates of AB11 indicate a slow attenuation with distance, more obvious for long distances R jb = 100 km. The lower bound of the range is provided instead by the regional model of ASB14 and Ketal15, whereas the CY14 and AB06 are representing the central models.
A comparion of the aleatory variability described by the total standard deviation of the selected GMPEs is included in the supplementary material (see Figure S7 and Figure S8) for a range of Vs 30 as a function of spectral periods. Noticeably, at short periods (T < 0.1 s) there is an increasing trend of the sigma values, with the highest value for CY14. The highest value of sigma is at 0.15 s for Ketal15 and above the T > 0.2 s the sigma reported by CF15 increases. Overall, the total sigma of the ASC models increases with spectral periods whereas the total sigma of the SCR models is constant and relatively lower. Regarding the Vs 30 , obviously the softer the soil (i.e., low Vs 30 values) the higher the median ground motion due to site amplification. The higher SA median values are overall given by the CF14, AB06 and AB11 relationships. CY14 yields larger values for Vs 30 around 200 m/s while for the remaining Vs 30 values, this model is in the range of Ketal15 and ASB14.
Following the main findings of the sensitivity analyses of other studies (Danciu et al. 2016;Delavaud et al. 2012;Sabetta et al. 2005) we assumed equal weights for all GMPEs, including those for the deep seismicity zone. These studies showed that there is little to gain from investing a significant effort on fine-tuning the weighting schemes because only large differences on the assigned weights have an impact on the hazard results. Thus, by assigning equal weights to these GMPEs, we simply acknowledge that these models are built upon recent strong motion datasets, their functional form is appropriate for describing the source, path and site relevant for our investigation, and finally that there is no bias towards a single ground motion model. A summary of the GMPE logic tree is listed in Table 5. As alluded to earlier, for the hybrid zones of Z1 and Z3 to account for the

Handling uncertainties in the source model: logic tree
To account for uncertainties in the most important parameters of the source model and the ground motion intensities for all the different earthquake scenarios, we created a logic tree structure. A logic tree is designed to handle the epistemic uncertainties of the model components represented as branching levels. The range of possible values of the parameters are described by a 3-point distribution defined by the upper, center and lower values. The associated weights follow this 3-point discretization, where the middle, or central, branch is assigned the most probable value of the parameter, which is reflected in the highest weight of 50%. The assigned weights of the upper and lower values are equally partitioned to 25%. If weights different than 25%, 50% and 25% are utilized for the lower, center and upper values, a justification is discussed in the manuscript but, otherwise, these 3 values are adopted. The logic tree for the source model parameters (Fig. 8) includes 3 earthquake sources: faults, distributed (in-buffer) sources and distributed (off-buffer) sources. The logic tree accounts for uncertainty in 5 main parameters: slip rate, b GR , M max , correlation distance values (C d ) for distributed seismicity rate calculations, and M thre . The slip rates are considered in 3 branches with the lower (SRL), mid (SRM) and upper (SRU) values. To account for the uncertainty in b GR fitting, we considered 3 sub-branches using the mean b GR (i.e., bM) and mean-b GR ± 0.05 (i.e., bL and bU). Note that the impact of the variations in the b GR is also reflected in the activity rates, i.e. a new fit using meanb GR ± 0.05 is used to compute the corresponding mean-a GR . In another branching level, we considered 3 alternative values for the maximum magnitude equal to the mean M max (i.e., MM) and mean M max ± 0.3 (i.e., ML and MU) magnitude unit. Three independent branches using C d = 10, 40 and 70 km were considered for the spatial variability of the activity rates in the distributed seismicity model. Herein, we assumed C d = 40 km as a suitable representative value of the spatial distribution of the events and, therefore, we have considered a larger weight of w = 0.5 for this central branch. On the other hand, we assigned an uneven weighting scheme for the lower and upper values of the range, namely w = 0.10 for C d = 10 km and w = 0.4 for C d = 70 km. Finally, we also considered 2 equally weighted branches, MtM and MtU, that account for the uncertainty in the choice of the threshold magnitude to be used in the Faults source model, with the following M thre = 6.0 and 6.5 values.

Results of the PSHA
We performed hazard calculations for multiple intensity measures including PGA and pseudo-spectral accelerations (SA) at several periods. The PSHA calculations are performed using the OpenQuake version 3.10.1. 4 The calculations are performed for the 3174 sites spread over the urban area with a grid spacing of approximately 500 × 500 m. The results of the PSHA are computed accounting for the soil conditions via the Vs 30 discussed in Sect. 2 and shown in Figure S6 in the supplementary material. Here we present some of the selected results in terms of hazard curves, disaggregation analysis, hazard maps and uniform hazard spectra (UHS). Based on the structure of the logic tree that includes both different source model parameters and different GMPEs, a total number of 7776 branches are created and run by OpenQuake. This is the result of 3 (for b GR ) × 3 (for Slip rates) × 3 (for M thre ) × 3 (for M max ) × 3 (for C d ) × 64 (for GMPE combinations). Note that the total number of GMPE branches (i.e., 64) is the result of the 8 (hybrid) × 4 (ASC) × 2 (Deep seismicity) GMPEs shown in Table 5. Figure 9 shows the mean as well as 16-th and 84-th percentiles of the hazard curves in terms of mean annual rate of exceeding PGA, SA(0.2), SA(0.75) and SA(1.0) for a selected site in the urban area of Isfahan. Figure 9 also shows in grey the hazard curves from the 7776 individual branches. Figure 10 shows the results of the disaggregation analysis (Bazzurro and Cornell 1999) for PGA corresponding to 100, 475, 975 and 2475 year return periods (T r ). As expected, at the shorter return periods the contribution from the events at farther distances is fairly large while for the longer return periods the trend reverses. For instance, while for a 100 years  Figure S6 of the supplementary material for T r = 100, 475, 975 and 2475 years covering the urban area of Isfahan. On the maps, the 15 municipal districts of Isfahan are also shown. The PGA values range between 0.14-0.20 g and 0.31-0.43 g, for 475 and 2475 year return periods, respectively. Figure 12 compares the PSHA estimates in terms of Uniform Hazard Spectra (UHS) obtained in this study with the seismic design spectrum proposed in Standard 2800 for soil types II and III. In this Figure the UHS are based on the 100, 475, 975 and 2475 years hazard levels. The grey lines in this Figure show the UHS obtained for individual sites classified as soil Type II or III while the mean of all sites of a given soil type is also presented by lines in color. Note that the T r = 475 years is assumed to be the basis for the design Fig. 9 Seismic hazard curves (mean and percentiles) for PGA, and spectral accelerations at 0.2 s, 0.75 and 1.0 s for a selected site at 51.604°E and 32.652°N in Isfahan urban area with soil conditions characterized as Soil Type III spectrum of Standard 2800 (Kohrangi et al. 2018). According to Standard 2800, the reference 475 years PGA on bedrock for Isfahan is A = 0.25 g (where g is the acceleration of the gravity). This comparison shows that, for both soil types, the design spectrum is significantly higher than all the 3174 475 years UHS computed in this study and, in our opinion, it provides conservative values for all spectral ordinates.

Sensitivity analysis
Among the crucial steps that analysts need to carry out and document for a defensible site-specific PSHA study (Budnitz et al. 1997), a fundamental one is a sensitivity analysis. The objective of this analysis is to quantify whether the range of epistemic uncertainties in the PSHA modeling is properly captured by the logic tree. Herein, we performed a sensitivity analysis to investigate the effect of our modeling assumptions in the definition of the source model and GMPEs branches in the logic tree. Again, these assumptions include, for instance, the values of the parameters such as b GR , M min , M max , C d , threshold magnitude in the fault model, M thre . Figure 13 shows a comparison of the mean and of the 16th and 84th percentiles of UHS (i.e., black lines) as well as the individual UHS obtained from the 7776 individual branches For more details about its definition, see Baker and Cornell (2006) of the logic tree (i.e., the grey lines) for a selected Isfahan site. Figure 13 also shows the UHS obtained from the combination of the EMME14 source model with the GMPE logic tree of this study (i.e., dotted green lines) as well as the combination of the source model of this study with the GMPE logic tree of EMME14 (i.e., dotted red lines). This comparison shows the relative agreement of our source model and GMPE logic tree with that of EMME14 as significant discrepancies in the mean UHS estimates are not observed. A more quantitative comparison on the impact of different values of the modeling parameters on hazard estimates is depicted in Fig. 14. This figure shows for a selected site with soil Type III, the percentage difference among the mean estimates, for example, of SA(0.2 s) for T r = 100, 475, 975 and 2475 years obtained from single branches of the source model and those from the full source model adopted in the logic tree structure of this study.
The b GR ± 0.05 variation has an impact of up to about ± 10% on the hazard for the lower IM levels which reduces to about ± 5% for those at longer return periods. M max ± 0.30 has a maximum impact of 8% on the values of SA(0.2) while the choice of M thre showed higher sensitivity on the results up to ± 12% for the lower bound (i.e., MtM) values and up to ± 8% for the upper bound (i.e., MtU) values. The sensitivity of slip rates branches (i.e., SRU, SRM and SRL) on SA(0.2) introduce a relatively high result variability ranging between about (± 18%,-± 18%). The lowest observed impact on hazard values (as low as ± 5%) is related to the choice of C d . The reason for such low sensitivity is related to the location of Isfahan where there are almost no events nearby. With reference to a comparison between  Fig. 6a, b, the background seismicity in a region of 40 km around the Isfahan urban area is almost invariant to C d for low return periods. Having said that, given that the source model generated here covers a relatively large area, the model could be freely used for regions far from Isfahan urban area and, as Fig. 6 suggests, a large impact of C d for farther sites is expected. Figure 14 also shows the sensitivity of the results to the choice of "hybrid" versus "ASC only" tectonic regime to be considered for zones Z1 and Z3. For the case of SA(0.2) shown in Fig. 14, this impact is lower than 10% but for other IMs, such as PGA, a larger impact is observed (see Figure S16 to Figure S19 in the supplementary material document). Finally, Fig. 14 also shows the impact of different branches of the GMPE logic tree on the hazard estimates. The upper and lower bound GMPEs switch their positions as a function of IM and IM level (Fig. 7). In Fig. 14, the GMPEs starting with "ASC" (e.g., "ASC: CY14"), represent the "ASC only" branches while, the ones starting with "Hybrid" (e.g., "Hybrid: SP16") are related to branches where Z1 and Z3 are assumed to be 34% (i.e., 2 GMPEs out of 6) SCR and 66% (i.e., 4 GMPEs out of 6) ASC. Finally, note that complementary sensitivity analyses for the hazard curves (see Figure S9 to Figure S16) as well as the corresponding tornado plots for other IMs and return periods (see Figure S16 to Figure S19) are provided in the supplementary material document of this manuscript.
While the definition of our logic tree was basically guided by our sensitivity analysis, the final selected logic tree is based on educated decisions made by the authors. Therefore, another group of researchers starting from the same raw datasets may very well end up with a materially different logic tree definition and weights. Nevertheless, with reference to our observations in this sensitivity analysis, we are confident that we have successfully captured the central value of all different variations in the modeling parameters, critical assumptions and weighting procedure considered in the logic tree.

Summary and concluding remarks
We performed a probabilistic seismic hazard analysis (PSHA) for the urban area of Isfahan (central Iran). We carried out a detailed source characterization for PSHA using the most recent catalog of historical and instrumental earthquakes and an original database of seismogenic sources produced by this study. Our seismogenic source comprises 23 individual parameterized active faults identified within about 200 km radius from the urban area. Based on these 2 datasets, we developed a seismogenic source model that includes both (i) a distributed seismicity model; and (ii) a faults and background seismicity model. The study area is subdivided into 7 distinct zones with distinguished geological and seismicity characteristics. We then subjectively classified these 7 zones into the 3 tectonic regions of active shallow crust, stable continental and deep seismicity regions. In the distributed seismicity model to efficiently account for the spatial distribution of the events we have redistributed the activity rates over the area source zones (Danciu et al. 2021;Poggi et al. 2020).
For each seismotectonic region, we selected a group of recently developed GMPEs to quantify the ground motion intensity on rock and soil conditions. The effects of local soil conditions are accounted for through Vs 30, which is as a proxy variable available in all of the GMPEs selected in this study. The Vs 30 values in Isfahan were estimated based on the available geotechnical test data and on the slope versus Vs 30 statistical correlation. Finally, we accounted for the uncertainty in the source modeling assumptions and ground motion intensity characterization via a logic tree approach.
We implemented the seismogenic source model in OpenQuake and performed the calculation for a number of intensity measures including PGA and pseudo-spectral acceleration for oscillator periods ranging from 0.1 to 4.0 s. We then generated hazard maps and uniform hazard spectra (UHS) for multiple intensity levels corresponding to 100, 475, 975 and 2475 years return periods. We compared the seismic hazard results obtained in this study with the results of the regional earthquake model for Middle East (EMME14) and the seismic design spectrum of the Iranian national design code. The results of our model are generally in good agreement with those of EMME14. The differences between the hazard estimates from the 2 models can be due to multiple reasons including the definition of the new fault model developed in this study, dissimilarities in both the seismogenic source delineations and the GMPEs adopted in the two studies.
The comparison of 475 yr UHS of this study with the seismic design spectrum of the Iranian code shows that the latter is significantly higher than the former and, therefore, we speculate that the design code is quite conservative. Although certain level of conservatism for low seismicity regions is understandable and somewhat expected in design codes, this level of conservativism could perhaps be transparently acknowledged so that the end users could better gauge the level of safety achieved by the structures they design using the code. Two alternative ways to do so could involve in low to moderate seismicity regions either prescribing to design for ground motions with return periods longer than 475 years, as is done in ASCE7 (2016) for CENA, or by adopting for design a percentile of the at 475 yr ground motion intensity (e.g., 84th percentile) higher than the median or the mean (Abrahamson and Bommer 2005). One such modifications seem necessary for many sites in Iran especially those located in the Sanandaj-Sirjan region where the so-called stable continental region of Iran is located.
Although to our knowledge this study provides the most comprehensive PSHA done for the city of Isfahan to date, we emphasize that it could be further enriched and updated by considering multiple aspects. First, a thorough study as well as field investigations are essential to better characterize the fault database in terms of their slip rates and their geometry. On the other hand, Isfahan is mainly located on a low seismicity region (i.e., Zone Z3 of this study) and even though using faults as seismic sources in the PSHA model is in our opinion highly recommended, whenever possible, such an approach requires caution and attention. In such quasi-rigid zones, steady displacements are rather low and thus accounting for large events with long interarrival times via fault modeling has been debated in the community. This is indeed an ongoing research topic  but dwelling on such a matter is beyond the practical scopes of this study. Second, in the current study we used Vs 30 as a proxy to account for the effects of local soil conditions on PSHA estimates at the surface. Although what was done here is the traditional approach, if better geotechnical data were to become available more precise analyses could be adopted to account in the PSHA for the amplification factors as well as their uncertainty at any given site (Bazzurro and Cornell 2004;Rodriguez-Marek et al. 2014) in the Isfahan urban area.