A probabilistic deformation-based seismic hazard model for Iran

Probabilistic seismic hazard analysis, as the most prevalent approach to evaluate earthquake hazard, is commonly based on earthquake catalogs. Although previous studies show that the recurrence time of large-magnitude (Mw> 7.0) events in Iran is more than ~ 1000–2000 years, the available instrumentally recorded earthquakes are limited to less than 100 years. With the idea of having another proxy for seismicity rates, we propose a methodology to evaluate activity rates from strain rates in a combination of regional estimates of β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and mmax\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\max }$$\end{document}. Considering the comprehensive deformation model of the Iranian Plateau, we found that the deformation-based occurrence rates are more than the catalog-based occurrence rates in all seismotectonic provinces in the Iranian Plateau. The ratios of the deformation-based to the catalog-based occurrence rate are between 1.34 and 2.35. For the first time, a probabilistic deformation-based seismic hazard model for the Iranian plateau is also developed and Peak Ground Acceleration (PGA) values for 10% and 2% POEs in 50 years are estimated. The highest levels of PGA are found in the Azerbaijan and Alborz seismotectonic provinces, where the highest value of strain rates is located. For all provinces, except the Zagros region, the spatial averages of PGA with 10% POE in 50 years from the deformation-based model are up to ~ 23% higher than those from a traditional probabilistic seismic hazard.


Introduction
Large-magnitude earthquakes have occurred in the Iranian Plateau during the long history of this ancient country. Probabilistic Seismic Hazard Analysis (PSHA) is the most prevalent approach to estimate seismic hazard and ground motion parameters. Seismic source model, including source geometry and earthquake occurrence rates, and Ground Motion Prediction Equations (GMPEs) are two main components of PSHA, which play a major role in the Probability of Exceedance (POE) of ground motion parameters at any given site. In the common practice of PSHA in Iran, a catalog of historical and instrumental earthquakes is used for the identification of seismogenic sources and quantifying the earthquake occurrence rates (Shabani and Mirzaei 2007;Khodaverdian et al. 2016b;Khoshnevis et al. 2017;Zafarani et al. 2020). Due to the limited number of large-magnitude instrumentally recorded earthquakes and lack of historical evidence for all faults, rates of large-magnitude earthquakes are usually extrapolated from rates of small-to-moderate ones in the prevalent PSHA (Ward 2007). However, Beauval et al. (2008) showed that, for an earthquake with a 475 years return period, a minimal 12,000 years observation independent of the seismic activity level is needed to estimate the occurrence rate with a 20% uncertainty. This long required period of observation could be justified by the fact that large-magnitude earthquakes affect the stress state and consequently the size and recurrence time of earthquakes on other adjacent faults, finally causing earthquake rate changes in time (Stein 1999). This fact is also supported by physics-based long-term earthquake simulations in Iran (Khodaverdian et al. 2016a, c), which showed that long-term seismicity rates should be considered for seismic hazard analyses. Paleoseismological studies (Le Dortz et al. 2009Dortz et al. , 2011Foroutan et al. 2012Foroutan et al. , 2014) also indicate that the recurrence time of large-magnitude earthquakes on a given fault in Iran could be more than 1000-2000 years. In conclusion, while the Iranian earthquake catalog is one of the longest in the world, considering the ancient history of the country, it is short relative to the long return period of large-magnitude earthquakes. Occurrence rates of large-magnitude events could be hence underestimated if only incomplete limited available historical evidence, archeological records (Ambraseys and Melville 1982, Berberian 1994, Berberian and Yeats 1999, Berberian and Yeats 2001, and short instrumental catalogs are considered. To address the issue of the lack of observed earthquakes, various types of data (e.g., geological fault slip rates and geodetic data) have been proposed to be used with the accompanying of incomplete earthquake catalogs (Godano and Pingue 2000;Ward 2007;Mazzotti et al. 2011;Danciu et al. 2018). A model that profits from a single type or a combination of deformation data sets can be an alternative way to the estimation of seismic activity rates (Ward 2007). For instance, in the Uniform California Earthquake Rupture Forecast (UCERF3), three deformation models were developed based on geodetic/geologic data sets and used in the evaluation of seismic activities (Petersen et al. 2007;Field et al. 2009Field et al. , 2014. As another example, Mazzotti et al. (2011) carried out the probabilistic seismic hazard assessment for continental western Canada based on geodetic data. They used 179 GPS station velocities to obtain regional strain rates. PGA values agree well with catalog-based results in only two zones and for other zones are 2-5 times larger than those derived from earthquake catalogs.
Several geological (Berberian and Yeats 2001;Bachmanov et al. 2004;Allen et al. 2011;Yaminifard et al. 2012) and geodetic investigations (Nilforoushan et al. 2003;Vernant et al. 2004b;Masson et al. 2005Masson et al. , 2007 in the Iranian Plateau have been carried out at different scales. GPS measurements were also initiated in 1999 and geodetic data is used to quantify the active tectonics of the Iranian Plateau and to evaluate the fault slip or strain rates. The first earlier models have built for the whole Iranian Plateau and provided present-day kinematics of the Iranian Plateau using geodetic data (Nilforoushan et al. 2003;Vernant et al. 2004b;Masson et al. 2005Masson et al. , 2007. On the other hand, some local models have been built for smaller seismically active regions (Tatar et al. 2002;Vernant et al. 2004a;Bayer et al. 2006;Hessami et al. 2006;Masson et al. 2006;Walpersdorf et al. 2006Walpersdorf et al. , 2014Tavakoli et al. 2008;Djamour et al. 2011;Mousavi et al. 2013). A comprehensive summary of these models is presented in Sect. 3. Considering the mentioned studies, there is an opportunity of developing a hazard model from all available deformation data for the Iranian Plateau.
Up to now, several researchers have built seismic hazard models for the Iranian Plateau with different assumptions (Shabani and Mirzaei 2007;Abdollahzadeh et al. 2014;Khodaverdian et al. 2016b;Khoshnevis et al. 2017;Zafarani et al. , 2020. Khodaverdian et al. (2016b) built a catalog-based seismic hazard model for Iran using the smoothed seismicity approach. Khoshnevis et al. (2017) provided a hazard model for the northern part of Iran based on observed seismicity using the same approach. The Earthquake Model of the Middle East (EMME) is a study, in which both geological data and observed seismicity have been used for seismic hazard assessment of the middle east, including the Iranian Plateau . The EMME consists of two independent seismogenic source models, area-and fault-based source models. Available geological information (i.e., fault trace, and fault slip rates) by that time, together with seismological data (i.e., focal mechanism parameters) have been considered in the development of the fault-based source model. In summary, apart from the EMME, which profits from some deformation data, all available hazard models for the whole or a part of the Iranian Plateau are based on observed seismicity.
In the present study, we are aimed to present the first probabilistic deformation-based hazard model for Iran. In the following section, a brief summary of the tectonic setting of the Iranian Plateau is first presented. Next, we summarize available deformation models for Iran and select the most comprehensive one, which is used for activity rate estimation and hazard model development. In section four, a methodology for earthquake occurrence rate estimation from the strain rate is presented, and then estimated rates for the Iranian Plateau are compared with catalog-based seismicity rates. In section five, a deformation-based seismic hazard model, developed based on the estimated earthquake occurrence rates, is presented. The characteristics of seismogenic sources and uncertainty considerations are also summarized. The outputs of our deformation-based seismic hazard model have been presented as hazard maps and hazard curves, and a comparison with other recent seismic hazard models is presented in the last part.

Tectonic setting of the Iranian plateau
The Iranian Plateau is located between the northward moving Arabian Plate and the Eurasian continent (Berberian 1981). The main tectonic regime has resulted from a convergence between the Eurasian and Arabian plates, which occurs at a rate of ~ 31 mm/year at the longitude of 52°E (Vernant et al. 2004b). Approximately rigid aseismic blocks (e.g., Central Iranian Block (CIB), South Caspian Block (SCB), and Lut Desert Block (LDB)) are surrounded by deforming belts (e.g., Zagros, Alborz, Kopeh Dagh, and Talesh) (Nowroozi 1976;Jackson et al. 1995;Khodaverdian et al. 2015).
Several studies categorized the Iranian Plateau into various seismotectonic provinces according to the tectonic, geologic information, and the quantitative observations of seismic events (Stocklin 1968;Takin 1972;Berberian 1976a,b Ansari et al. 2009). The seismotectonic map, proposed by Mirzaei et al. (1998), is one of the commonly used categorizations, in which the Iranian Plateau is divided into five major tectonic zones, including Azerbaijan-Alborz, Kopeh Dagh, Zagros, Central-East Iran, and Makran. We considered the map of the seismotectonic provinces presented by Mirzaei et al. (1998) as a base and separated Azerbaijan-Alborz and Central-East Iran into two parts so that seven seismotectonic provinces (i.e., Azerbaijan, Alborz, Kopeh Dagh, Zagros, Central Iran, Eastern Iran, and Makran) are finally considered for the Iranian Plateau (Fig. 1). Also shown in Fig. 1 is a homogenized earthquake catalog compiled by Shahvar et al. (2013) and updated by Khodaverdian et al. (2016b). All earthquakes in the Iranian plateau are shallow crustal events, with focal depths less than 40 km, except events in the Makran subduction zone, southeast Iran, where deeper earthquakes (in-slab or interface earthquakes) could occur (Mirzaei et al. (1998). The Makran region is excluded from the current study as shallow crustal deformation data could not be directly informative for seismicity estimation. The major faults and seismic characteristics of provinces are summarized in the following.
Alborz seismotectonic province is one of the most active regions in the north of the Iranian Plateau. The Alborz mountain belt is 600 km long and stretches from the Talesh mountains in the west to the Kopeh Dagh mountains in the east along the southern side of the Caspian Sea. Between Central Iran and the SCB, there is oblique convergence, which is  Mirzaei et al. (1998) observed in the form of the left slip along the Astaneh, Firoozkouh, Mosha, and Taleghan faults and perpendicular movement along the Caspian and Alborz faults (Ritz et al. 2003;Bachmanov et al. 2004;Nazari 2006;Nazari et al. 2009). The Khazar and North Alborz faults (Bachmanov et al. 2004) are thrust faults located in the north part of the province and some major active faults (e.g., Pishva, Eshtehard, and Ipak faults) are delineated in the south. The North Tehran fault with a 110 km length, is the most hazardous fault near the capital city of Tehran (Ritz et al. 2012). Several historical earthquakes in 743, 958, 1177, 1665, and 1830 are reported for the area close to the city of Tehran (Ambraseys andMelville 1982, Berberian 1994).
Zagros seismotectonic province is the most active seismic region in the Iranian Plateau, which absorbed a large amount of the total deformation from the Arabian-Eurasia convergence. Zagros fold-and-thrust belt with 250-400 km width is extended from eastern Turkey to the Strait of Hormoz. The northwest part of Zagros is affected by orogen-parallel thrust faults, and the large Main Recent Fault (MRF), showing partitioning of the oblique motion to the north (Berberian 1995). The right-lateral strike-slip faulting along the MRF is approved by tectonic geomorphology and focal mechanism solutions of earthquakes (Berberian 1995;Berberian and Yeats 2001;Talebian and Jackson 2004;Yamini-Fard et al. 2006). The southeast part of Zagros is characterized by lots of shallow thrust faults and the NW-SE strike of the Zagros fold range turns to almost W-E in this part, east of 52°E, and the belt becomes perpendicular to the Arabia-Eurasia shortening direction.
Azerbaijan seismotectonic province at the northwest of Iran is mostly affected by the North Tabriz fault system with a dominant right-lateral faulting mechanism and NE-SW trend (Jackson and McKenzie 1984;Berberian and Yeats 1999;Karakhanian et al. 2004). Tectonics of the Turkey-Iran-Caucasus region are investigated in several studies (e.g., (Copley and Jackson 2006)) that revealed this region is characterized by partitioning of oblique plate convergence. Earthquakes focal mechanism solutions (Jackson 1992) show thrust faulting in the eastern part of the Greater Caucasus, and the right lateral motion accommodated in the Turkish-Iranian Plateau on strike-slip fault systems (e.g., Chalderan and North Tabriz faults). From 1721 to 1786, during 65 years, the North Tabriz fault has ruptured three times, which includes the M s 7.7, 1721, the M s 7.7, 1780, and the M s 6.3, 1786 earthquakes (Ambraseys andMelville 1982, Berberian 1994). More recently, extensive damages and human losses are caused by the M w 6.4 and 6.3, 2012 Ahar-Varzaghan double earthquakes.
The Kopeh Dagh province, forming the boundary between the Iranian Plateau and Eurasian plate. The Kopeh Dagh range is confined by the NW-SE trending Main Kopeh Dagh Fault (MKDF) in the northwest and thrust faults in the southeast. The M w 7.3, 1948 earthquake has been ascribed to the MKDF. In the central part of the Kopeh Dagh range, an array of right-lateral strike-slip faults with NNW-trending (i.e., the Baghan-Quchan fault system) obliquely crosses the range and allows Central Iran to move northward with respect to the stable part of Eurasia (Shabanian et al. 2009). There are also some major faults (i.e., Neyshabur and the Binalud reverse faults), which were responsible for some large historical earthquakes with M w > 7, that struck near the Neyshabur city between 1209 and 1405 (Berberian and Yeats 1999).
The central and eastern seismotectonic provinces of Iran encompass aseismic rigid blocks (e.g., CIB), which have been surrounded by seismically active belts (e.g., Alborz and Zagros) (Nowroozi 1976;Jackson et al. 1995;. LDB is also surrounding by some fault systems (e.g., Doruneh and Dasht-e-Bayaz faults). These faults are seismically active and responsible for the major destructive earthquakes in the region (e.g., the M w 7, 1968, the M w 7.3, 1979, and the M s 7.3, 1978 earthquakes) (Berberian and Yeats 1999;Aghanabati 2004;). Infrequent and large-magnitude earthquakes with recurrence intervals of more than several hundred years are characteristics of Central-Eastern Iran (Berberian and Yeats 1999;Berberian 2014). Although the region experienced infrequent earthquakes during the approximately 2750 years of observation time, according to paleoseismological studies (Le Dortz et al. 2009Dortz et al. , 2011Walker et al. 2010;Foroutan et al. 2012Foroutan et al. , 2014, several large-magnitude earthquakes occurred in the region during the Late Pleistocene and Holocene. These studies have shown that faults without historical seismic evidence in Central-Eastern Iran are responsible for large-magnitude (M w ~ 7) earthquakes during the Holocene and the contribution of recently silent faults to seismic hazard assessment must be considered.

Deformation and strain rate models for Iran
Several deformation models have been developed for the Iranian plateau based on geological and geodetic information and these models could be classified into two categories. The first category is global models that provide a large-scale perspective of the Iranian Plateau's present-day kinematics, plate motions in the Middle East, and gross deformation across various active belts (Nilforoushan et al. 2003;Vernant et al. 2004b;Masson et al. 2005Masson et al. , 2007. The second category is local models that have been built for smaller seismically active regions in the Plateau using block modeling techniques to provide detailed active deformation of different fault systems in Iran (Tavakoli et al. 2008;Djamour et al. 2011;Mousavi et al. 2013;Walpersdorf et al. 2014). However, some drawbacks of the block modeling technique are: (1) some blocks are poorly represented by the rigid block model due to the existence of small faults in the region, (2) an array of active strike-slip faults are simplified and modeled by only one boundary fault, and (3) the GPS stations are still too sparse, particularly around the secondary faults, and slip rate for few secondary faults is derived, while lots of minor faults in the Iranian plateau need to be considered in a proper strain modeling (Khodaverdian et al. 2015). Khodaverdian et al. (2015) presented the first comprehensive deformation model, in which the long-term crustal flow of the Iranian Plateau is computed by using various data sets, including fault traces, geologic fault offset rates, GPS velocities, principal stress directions, and velocity boundary conditions. As input, the long-term geological offset rates for 33 of 171 fault traces were collected based on the available data on the relative displacement of geologic features. Also, geodetic velocities of 239 GPS benchmarks, together with data from the world stress map were considered. The main outcome of the model is strain rates and fault slip rates, which are the most up-to-date evaluation at the regional scale. The estimates of fault slip rates have been also validated with slip rates resulted from analyzing geodetic velocities, or geological/paleoseismological studies, which have not been used in the model calibration.
When solely GPS data are used in developing a deformation model, an important step is required to limit the perturbing impacts of transient and surficial processes in the earth. Bird and Carafa (2016) and Carafa and Bird (2016) demonstrate that these effects of transient and surficial signals are significantly minimized when the augmented covariance matrix is utilized in the objective function for a deformation model, primarily driven by GPS data. They also showed that the successful application of their method should be investigated by comparing the GPS-based deformation model with a comprehensive model, in which several datasets (e.g., plate-tectonic velocity boundary conditions, geological fault traces and their slip rates, and the most compressive horizontal principal stress directions) are used. Since the Khodaverdian et al. (2015)'s model is developed based on several different long-term input datasets, and since their estimated rates are consistent with independent long-term geological rates, their output (i.e., strain rates) provides a picture of long-term deformation in the Iranian Plateau and is the viable candidate for further calculation of moment and seismic activity rates (see Sect. 4.2 for more details).

Seismic activity rate from deformation model
In the following section, we first present a methodology to evaluate the earthquake occurrence rate from the strain rate. The total moment rates for the Iranian Plateau are then evaluated from two available strain rate models by the implementation of the proposed method. Next, the estimated earthquake occurrence rates are compared with catalog-based seismicity rates. We also carried out a sensitivity analysis to study the effect of seismogenic layer thickness on the estimated earthquake occurrence rates from strain rates.

Methodology
The first step to estimate earthquake occurrence rates is to evaluate the deformation-based total moment rate (Ṁ D 0 ) from strain rates using the following equation (Liu and Bird 2008): where G is the average rigidity of the crust, A is the area of the source zone, and ⟨cz⟩ is the mean coupled seismogenic thickness of the crust.̇1 , ̇2 and ̇3 are the principal components of the strain rate tensor with ordering ̇1 ≤ ̇2 ≤ ̇3 . As presented in Eq. (2), by assuming double truncated Gutenberg-Richter (GR) distribution, the deformation-based occurrence rate of earthquakes with equal to or larger than minimum magnitude, D (m 0 ) , can be given by the following equation (e.g., (Anderson and Luco 1983)): where m 0 and m max are the minimum and maximum magnitude, respectively. M 0 (m) is the seismic moment, which is released by an earthquake with magnitude m , and can be evaluated by the well-known relation of Hanks and Kanamori (1979) as following: where c equals 1.5 and d takes the values 9.05 or 16 to estimate the seismic moment in units of N.m or dyne.cm, respectively. f m (m) represents frequency-magnitude probability density function for double truncated GR distribution: where = b ln 10 is the slope of the GR distribution.
(1, 2, 3, 4), Ṁ D 0 , which is obtained from the strain rates, could be converted to D (m 0 ) given that , m 0 , and m max are available for the study area. Spatial variation of and m max are not as high as seismic activity rate; and m max are mostly presented as regional parameters and can be estimated from the catalog of observed earthquakes for a given region (Kijko et al. 2016). In other words, fault-specific or source-specific -values are rarely used in hazard studies due to the lack of sufficient data. It is very common to construct super zones and evaluate a regional -value and assign it to all inner source zones (see e.g., Field et al. 2014Danciu et al. 2018;). The lower bound or minimum magnitude ( m 0 ) is chosen on the basis that events with lower magnitude are mostly ignored by engineers. A review of practice reveals that values of m 0 generally lie in a range of 4-5. In the present study, the minimum magnitude is set as 4.0. In summary, the regional estimate of and m max from the statistics of observed seismicity can be hence a suitable value to be considered as an input for the proposed methodology and the deformation-based seismicity rate, D (m 0 ) , can be eventually evaluated by using all the aforementioned equations.

Total moment rates for the Iranian plateau
Deformation-Based Strain Rate Model (DBSRM) presented by Khodaverdian et al. (2015) is the most up-to-date and comprehensive model for the Iranian Plateau. The strain rates from that model are converted into the deformation-based total moment rates ( Ṁ D 0 ) (see Fig. 2a) by using Eq. (1) and assuming G = 27.7 GPa, as suggested by Bird and Kagan (2004), and ⟨cz⟩ for a continental convergent boundary from Bird et al. (2010). The highest values of Ṁ D 0 are concentrated along faults, located in major seismically active belts (e.g., the Zagros, Alborz, and Kopeh Dagh mountains). Moment rates in central and eastern Iran (e.g., the CIB and LDB) are generally low while high moment rates can also be seen along the faults surrounding the LDB. Although DBSRM covers the whole Iranian Plateau, the current study focuses on the crustal seismic sources and does not include the Makran subduction region, located in southeastern Iran. Moment rate estimates for that region will be presented in future papers by combining the results of the ongoing research projects aiming to capture a better picture of the seismogenic layer in that region. Another available deformation model covering the whole Iranian Plateau is the Global Strain Rate Model (GSRM), which is the purely geodetic-based model by (Kreemer et al. 2014). Figure 2b shows the spatial distribution of the total moment rates derived from GSRM using Eq. (1) for the whole Iranian Plateau. Comparison of Fig. 2a with b shows that both models evaluate high moment rates in the Alborz, Azerbaijan, Kopeh Dagh, and Zagros seismotectonic provinces. For instance, total moment rate values, derived from DBSRM and GSRM, for the area close to the city of Tehran are 2.31E+13 (Nm/yr/km 2 ) and 1.69E+13 (Nm/yr/km 2 ), respectively. As another example, the rates derived from DBSRM and GSRM, for the area close to the city of Tabriz, located in Azerbaijan, are 13.17E+13 (Nm/yr/km 2 ) and 4.13E+13 (Nm/yr/km 2 ), respectively. Comparison of Fig. 2a with b also shows that DBSRM generally gives higher rates for Central Iran and Eastern Iran since the geological slip rates for active faults in these provinces (e.g., the Dehshir, Anar, Nayband, Nehbandan, and Nosratabad faults (see Fig. 5 for the fault locations)) are also considered in DBSRM.
GSRM and DBSRM are different in terms of scale and inputs. While GSRM is developed based on geodetic data on a global scale, DBSRM is built at a local scale and covers only the Iranian Plateau, which provided the opportunity of profiting from various deformation datasets (e.g., geological and geodetic data and principal stress directions). Considering the finer details used in DBSRM, moment rates derived from DBSRM are selected to be used in the next step of earthquake occurrence rates evaluation.

Deformation-based activity rate and its comparison with observed seismicity
As described in Sect. 4.1, and m max are required seismicity parameters in the evaluation of the deformation-based earthquake occurrence rates. These parameters are calculated by the maximum likelihood procedure of Kijko et al. (2016) from the observed seismicity after removing aftershocks and foreshocks. The homogenized earthquake catalog, compiled by Shahvar et al. (2013) and updated by Khodaverdian et al. (2016b), is used. M w is chosen as a reference magnitude, and the entire earthquake catalog is declustered using Gardner and Knopoff (1974) approach. As the catalog consists of historical and instrumental earthquakes in the Iranian plateau, the historical and instrumental parts are precisely separated using Stepp (1972) method. The instrumental section is also examined to find the various completeness intervals by using the maximum curvature method (Wiemer and Wyss 2000). Generally, the earthquake catalog available for the Iranian plateau can be roughly divided into three periods. The first period contains evidence of historical earthquakes up to the beginning of the twentieth century. The second period (1900 to ∼ 1970) is a transition between the period of historical data and the modern instrumental period, although the installation of the first seismographic instruments in Iran was performed during the earlier decades of the past century. The third period starts in 1970, since which the obtained seismic data is purely instrumentally recorded due to the adequate density of the regional network. Having historical and instrumental parts of the earthquake catalog, together with the completeness magnitude for each time interval from the maximum curvature method, we estimated long-term catalog-based occurrence rates of M w 4.0+ earthquakes for seismotectonic provinces ( P C 4 ). Also, we calculated the corresponding long-term catalog-based seismic moment rates ( PṀ C 0 ) by integrating the seismic moment over the magnitude-frequency distributions, derived from the procedure of Kijko et al. (2016). Those estimates, together with and m max , for all seismotectonic provinces in the Iranian Plateau are summarized in Table 1. Considering the regional estimates of , m max and assuming the minimum magnitude of 4.0 ( m 0 = 4.0 ), we calculated deformation-based occurrence rates of earthquakes with magnitude equal to or larger than M w 4.0 ( P D 4 ) for all provinces using the proposed method in Sect. 4.1. As shown in Table 1, the deformation-based earthquake occurrence rates are higher than the catalog-based occurrence rates; P D 4 ∕ P C 4 varies from 1.34 to 2.35. The lowest amount of P D 4 ∕ P C 4 is observed for the Zagros province. In this province, earthquakes with relatively short recurrence times occur (Talebian and Jackson 2004;Kalaneh and Agh-Atabai 2016) and, as a result, less discrepancy between deformation-based and catalog-based earthquakes rates is expected. The maximum amount of P D 4 ∕ P C 4 corresponds to Central Iran, where infrequent and large-magnitude earthquakes occur with long recurrence times (Berberian and Yeats 2001). For the Central Iran and Eastern Iran provinces, the values of P D 4 ∕ P C 4 are estimated at 2.35 and 1.70, respectively. The values of P D 4 ∕ P C 4 for the Azerbaijan, Alborz, and Kopeh Dagh provinces are varying between the minimum and maximum values, which is consistent with the fact that the recurrence times of the large-magnitude earthquakes in these provinces are less frequent than Zagros and more frequent than Central Iran and Eastern Iran. To have a full picture of the comparison between catalog-based and deformation-based earthquake occurrence rates, the Frequency-Magnitude Distribution (FMD) for all seismotectonic provinces is depicted in Fig. 3.

Sensitivity analysis for seismogenic layer thickness
The seismic moment rate is dependent on the thickness of the seismogenic layer. Using sensitivity analysis of the seismogenic thickness, we investigated the effect of ⟨cz⟩ variations on the P D m 0 and GR distribution. In Table 2, P D m 0 ,⟨cz⟩ , occurrence rates of the M w 4.0+ earthquakes for each seismotectonic province for mean coupled seismogenic thicknesses of 10, 15, and 20 km are presented; these three values are based on the depths of observed seismicity (Engdahl et al. 2006). As an example, the deformation-based GR distributions for the Zagros province with different mean coupled seismogenic layer thicknesses of 10, 15 km, and 20 km are shown in Fig. 4. Also shown in Fig. 4 are the catalog-based GR distribution and its corresponding uncertainty. It can be interpreted that the uncertainty of ⟨cz⟩ is comparable with the uncertainties of and P C 4 , which are derived from observed seismicity. Seismogenic thickness has, therefore, a considerable effect on the estimated activity rate values and its uncertainty should be consequently considered in the seismic hazard analyses.

Seismic hazard assessment
Deformation-Based Hazard Model (DBHM) and its details are presented in the following section. We defined the geometry and characteristics of seismogenic sources using two categories of data (Sect. 5.1) and seismic activity rates for these sources are assigned from long-term estimates we have from the deformation model. In Sect. 5.2, the GMPEs used in the model are presented. Logic tree method for capturing the uncertainties of mean coupled seismogenic thickness, maximum magnitude, -values, and GMPEs is considered and described in Sect. 5.3. Fig. 3 The frequency-magnitude distributions for different seismotectonic provinces

Seismic area sources
From the selected comprehensive deformation model for the Iranian Plateau, the strain rate field is available for a grid of 0.2 • × 0.2 • cells. Square area sources, with the same size and spacing, are considered to cover the whole Iranian Plateau (Fig. 5). The estimate of , and m max for the Iranian Plateau on a grid of 1 • × 1 • degrees is available from Khodaverdian et al. (2016b). Although the regional estimate of , and m max for each seismotectonic province is provided in Sect. 4.3, we used estimated values on a grid of 1 • × 1 • to better capture  . 4 The frequency-magnitude distributions of earthquakes for the Zagros seismotectonic province the spatial variation of these parameters and allocated them to square area sources based on the nearest neighbor rule. Using these values and assuming the minimum magnitude of 4.0 ( m 0 = 4.0 ), we then calculated the deformation-based occurrence rate ( D 4 ) for each square area source, according to the procedure described in Sect. 4.1. It is worth mentioning that uncertainties of the regional estimate of and m max , provided by Khodaverdian et al. (2016b), are also taken into the account (please refer to Sect. 5.3 for more details). D 4 is estimated based on the strain rate, encompasses both seismic and aseismic deformation. With the goal of seismic hazard assessment, it is necessary to estimate the seismic fraction of the deformation-based occurrence rates; the portion of aseismic deformation should be removed from the occurrence rate before the seismic hazard assessment. We assumed seismic deformation fraction ( ) could be estimated from ratio between the long-term catalog-based occurrence rate and the long-term deformation-based occurrence rate when a big area (e.g., seismotectonic provinces) is considered; In other words, for each seismic province is assumed to be equal to P C 4 ∕ P D 4 , which has already estimated and presented in Sect. 4.3. (please see Table 1). In summary, the seismic activity rate for each 0.2 • × 0.2 • area source is calculated by multiplying D 4 by estimated for the corresponding seismotectonic province.
Other characteristics of seismic sources are usually evaluated based on tectonic, geological, and seismological evidence. One of the important parameters of seismic source modeling is the depth of the seismogenic layer. Since the tectonic setting throughout a province is similar (Hong et al. 1996), we assign the identical hypocentral depths to area sources located in the same provinces. We used the map of seismotectonic provinces of the Iranian Plateau presented by Mirzaei et al. (1998) (Fig. 1) and hypocentral depth for the different provinces is evaluated from available studies. For example, Mooney et al. (1998) suggested that seismicity in eastern Iran is restricted to the upper continental crust with a seismogenic thickness of ∼15 km (Maggi et al. 2000), also moderate-to-large magnitude earthquakes occurred with a median depth of 12 ± 5 km (Zare et al. 2014). More recently, Engdahl et al. (2006) have reported precise earthquake depth distribution for the seismotectonic provinces of the Iranian Plateau. For the Alborz province, the range of earthquake depths is 5-30 km, while more earthquakes occur at a depth of 5-15 km and these values for Eastern Iran are 0-20 km and 5-15 km, respectively (Engdahl et al. 2006). Zare et al. (2014) presented a depth of 11 km for most of the earthquakes that occurred in the Alborz and Azerbaijan. A depth of 15 km for most of the earthquakes in the Kopeh Dagh and Central Iran and a depth of 13 km for the earthquakes in the Zagros are reported (Zare et al. 2014). According to the mentioned studies, we considered the hypocentral depth of 10 km, 15 km, and 20 km with the probabilities of 0.333, 0.334, and 0.333, respectively, for the Zagros province. For all other seismotectonic provinces, we considered the hypocentral depth of 7.5 km, 12.5 km, 17.5 km, 22.5 km, 27.5 km with the probabilities of 0.375, 0.375, 0.125, 0.0625, and 0.0625, respectively.
The surrounding zones around the mapped active faults are one of the major seismic sources as the concentration of seismic moment within narrow zones along the faults is shown in Fig. 2a. The characteristics of the corresponding faults should be considered in the seismic source modeling of these zones. Active faults maps of Iran are provided in several studies (e.g., Hessami et al. (2013) and Danciu et al. (2018)) (see Fig. 5) and the fault characteristics are reported in different investigations (e.g., Khodaverdian et al. (2016a) for Azerbaijan, Khodaverdian et al. (2016c) for Eastern Iran, and Bachmanov et al. (2004) for the Zagros and Central Iran). Those studies and the references therein have been used to define near-fault zones with a width of 15 km on each side along the active faults and the predominant faulting plane characteristics (i.e., strike, dip, and rake) of near-fault zones. As shown in Fig. 5, a total of 18, 1, 10, 9, 6, and 19 near-fault zones (shaded areas) in the Alborz, Zagros, Azerbaijan, Kopeh Dagh, Central Iran, and Eastern Iran provinces are considered, respectively. It is worth mentioning that the symmetric near-fault zones may lead to underestimating the hanging wall effect of ground motion in the case of low-dipping faults, which is not common in the Iranian Plateau.
In the Zagros province, due to the presence of several hidden secondary/minor faults (Berberian 1995;Berberian and Yeats 1999), considering separate near-fault zones is not plausible. In addition, due to the different geometrical characteristics of faults in the northwest and southeast parts of this province (Berberian 1995;Berberian and Yeats 2001;Talebian and Jackson 2004;Yamini-Fard et al. 2006;Tavakoli et al. 2008), we defined two large zones for the northwest and southeast part of the province. These two zones together with the near-fault zone for MRF cover the whole province (Fig. 5). If a given square area 0.2 • × 0.2 • source is located within a near-fault zone, the characteristics of the corresponding near-fault zone are assigned to the area source. For area sources located outside of near-fault zones, the faulting characteristics are sampled from probability distribution functions of the strike, dip, and rake, estimated from all faults located in the corresponding seismotectonic province.

Ground motion prediction equations
GMPE is a key element in PSHA as it estimates the intensity level at a particular site given a specified earthquake. Cotton et al. (2006) point out some criteria that could be used for GMPE selection from a list of candidate models. They suggested that the GMPE model should have adequate documents obtained from a relevant tectonic regime and should be published in an international peer-reviewed journal. The model outcome range should be proper for engineering applications and has a suitable functional form. Bommer et al. (2010) modified the Cotton et al. (2006) criteria and provided a more complete set of criteria. For example, they presented that the model must use appropriate definitions for parameters (e.g., magnitude and distance) and has the consideration of site effect with the average shear-wave velocity for the upper 30-m depth, denoted as V s30 .
Based on previous experiences (e.g., Mousavi et al. (2012), Zafarani and Mousavi (2014), Zafarani and Farhadi (2017), and Shoja-Taheri et al. (2010)), in which performance of the local, NGA or regional GMPEs for the Iranian Plateau have been examined, a combination of regional, NGA and local GMPEs is a rational choice for hazard analysis in Iran, from the viewpoint of capturing epistemic uncertainty.
Based on data-driven statistical evaluations, it can be concluded that local models usually provide a good explanation of recorded strong motion data in Iran (e.g., Mousavi et al. (2012), Zafarani and Mousavi (2014), and Zafarani and Farhadi (2017)). However, due to the scarcity of near-field records, their performance in the prediction of near-field motions of large earthquakes is a matter of debate. The database of regional and NGA GMPEs are richer in the near-field region, thanks to gathering data from diverse shallow crustal regions. For example, the NGA-WEST2 contains data from large earthquakes that occurred in Alaska, Mexico, Greece, Turkey, Iran, Taiwan, China, New Zealand, Japan, and Italy.
Finally, we select five local, regional, and global GMPEs for seismic hazard assessment. The first one is Zafarani et al. (2018), a local GMPE, which is developed based on a comprehensive, compiled good-quality strong-motion database of the Iranian earthquakes. The second one is Kale et al. (2015), a regional GMPE, which is developed for Turkey and Iran to investigate the possible regional effects on ground-motion amplitudes in active shallow crust earthquakes. This GMPE is based on a subset of the compiled strong-motion database of the Earthquake Model of the Middle East (EMME) project . The third GMPE is Zhao et al. (2006), which is one of the global GMPEs that was presented based on a large number of strong ground-motion records in Japan. Abrahamson et al. (2014) and Idriss (2014) are global GMPEs that are developed in the NGA-West2 project (Bozorgnia et al. 2014) for shallow crustal earthquakes in active tectonic regions.
Finally, it should be noted that all available GMPEs in Iran are developed based on the simple as the recorded geometric mean of two horizontal components. However, it has been shown by Boore (2010) that as recorded geometric mean of two horizontal components is very similar to the median spectral accelerations over all orientations (RotD50), which NGA-West2 models predict (see Fig. 4 in Boore 2010). As shown in Fig. 6, recorded observed ground motion for the Finn, 2021 dual earthquakes (see Fig. 5, where their epicenter are shown) are compared to five selected GMPEs. General compatibility between the geometric mean of the east and north components of Finn dual earthquakes and selected GMPEs for M6.2 with reverse focal mechanism and Vs30 = 500 m/s is observed.

Uncertainties consideration
For taking into account possible epistemic uncertainties, the logic tree approach is used. A logic tree consisting of four branching levels (mean coupled seismogenic layer thicknesses, seismicity parameters ( and m max ), and GMPEs) is proposed (please see Fig. 7). The branch weights in a logic tree framework represent the degree-of-certainty or degreeof-belief of experts in alternative models and parameter values. According to Bommer (2012), the population of logic trees with models that represent both the best estimates on the basis of what is known and alternative models that represent what could occur in light of our lack of knowledge, requires expertise and experience in both the relevant earth science disciplines and in the characterization of uncertainty.
The first branching level shows the uncertainties of the mean coupled seismogenic layer thicknesses. As shown in Sect. 4.4, variation of seismogenic layer thickness has a considerable effect on the activity rate values and should be considered in the seismic hazard analysis. We considered ⟨cz⟩ min = 10 km , ⟨cz⟩ ave = 15 km , and ⟨cz⟩ max = 20 km with the weights of 0.25, 0.5, and 0.25, respectively. The second branching level demonstrates the uncertainties of regional estimates of . In this branch, ave − , ave , and ave + with the weights of 0.25, 0.5, and 0.25 are considered, respectively. The third branching level belongs to the uncertainties of the maximum magnitude. The weights  Kale et al. (2015), Zafarani et al. (2018), and Zhao et al. (2006). The geometric mean of two horizontal components are shown. Median and ± 1 sigma values from GMPEs are calculated for M6.2 with reverse focal mechanism and Vs30 = 500 m/s of 0.25, 0.5, and 0.25 are considered for m max,ave − m max , m max,ave , and m max,ave + m max , respectively. The fourth branching level belongs to the uncertainties of the GMPEs. We considered the same weight for all GMPEs used. Here, we attempt to assign numbers to our degrees of belief, which are by nature personal and indefinable, and for which there are neither tests nor measurements. However, in addition to the "expert opinion" approach, we also consider a statistically-based scheme to assign the logic tree weights of GMPEs (Mousavi et al. 2012;Zafarani and Mousavi 2014;Zafarani and Farhadi 2017).

Results
The seismic hazard computation is performed with the OpenQuake hazard engine (Pagani et al. 2014), which is an open-source package developed by the GEM for seismic hazard and risk analyses. We estimate PGA at 10% and 2% POEs in 50 years, which correspond to 475 and 2475 years return periods, respectively. Hazard assessment results are presented on bedrock with the V s30 of 760 m/s (V s30 represents the average velocity of shear waves in the top 30 m of soil). It is a common practice to truncate the lognormal distribution of ground-motion residuals at 2 or 3 times of the standard deviation (Bommer and Abrahamson 2006). Truncating above this level has little effect on the hazard curves for the return periods generally used in engineering design; therefore, a truncation at 3 standard deviations is applied. In the following, hazard maps and hazard curves are presented and compared with the results of the other probabilistic seismic hazard analyses for the Iranian Plateau.

Hazard maps
Hazard maps of the mean PGA value for 10% and 2% POEs in 50 years are presented in Fig. 8a and b, respectively. The spatial average values of the PGA with 10% POE for the Azerbaijan, Alborz, Kopeh Dagh, Central Iran, Eastern Iran, and Zagros seismotectonic provinces are 0.40 g, 0.36 g, 0.30 g, 0.24 g, 0.25 g, and 0.26 g, respectively. Also, the standard deviations of the PGA values for 10% POE are 0.1 g, 0.09 g, 0.07 g, 0.07 g, 0.06 g, and 0.05 g, respectively. According to hazard maps derived from DBHM, the highest PGA are observed in the Azerbaijan and Alborz provinces, where the highest values of strain rates are also located. Also, in the Kopeh Dagh province, high PGA values have been seen in a strip, which is extended from the east of the Alborz to Mashhad city. The strip covers Robate-Qarabil, Rivand, Neyshabur, and Binalud faults, for which high strain-based activity rates have been seen (see Fig. 5 for the fault locations).
The EMME is the most up-to-date regional model, in which the whole Iranian Plateau and both geological data and observed seismicity have been used . These facts make the EMME a decent candidate for comparison purposes. To illustrate the spatial pattern of the differences between DBHM and EMME models for various POEs, the ratio of the mean PGA value from DBHM to EMME is calculated and presented in the form of the natural logarithm ( Fig. 9a and b). For Azerbaijan, the average values of PGA for 10% and 2% POEs in 50 years are 23% and 9% higher than the PGA values derived from EMME, respectively. For Alborz, the average values for 10% and 2% POEs are 18% and 1% higher than the PGA derived from EMME, respectively. For large parts of the Alborz and Azerbaijan provinces, where earthquakes with moderate-to-large recurrence times happen, the PGA values estimated from DBHM are higher than the corresponding ones from EMME. According to Berberian and Yeats (1999), the Alborz and Kopeh Dagh provinces are characterized by largemagnitude shallow earthquakes separated by quiescent periods of 3000-5000 years on frontal reverse faults and some strike-slip faults. In contrast, the Zagros province is characterized by rapid uplift, high seismicity, and moderate-magnitude earthquakes on several blind reverse faults with historical recurrence intervals of 300 to 500 years (Berberian and Yeats 1999).
In the Kopeh Dagh province, also the average PGA values for 10% and 2% POEs are 13% and 2% higher than the PGA values derived from EMME, respectively. Although our estimates of PGA for most parts of the Iranian Plateau (e.g., Azerbaijan, Alborz, and Kopeh Dagh provinces) are higher than EMME model estimates, PGA values derived from the EMME for the Zagros province on average are 20% and 28% higher than the values derived from DBHM for 10% and 2% POEs, respectively. It is worth mentioning that the EMME source model consists of two layers of the seismogenic sources (i.e., shallow sources and deep sources with a hypocentral depth of higher than 50 km) in the Zagros province, and two GMPEs (i.e., Youngs et al. (1997) and Lin and Lee (2008)) developed for subduction zones with equal weights are assigned to deep sources. These GMPEs provide higher PGA than our considered GMPEs.
For the Central and Eastern Iran provinces, on average the PGA values for 10% POE in 50 years from DBHM are 19% and 13% higher than PGA values derived from EMME, respectively. The main reason for these differences could be the lack of enough data for large-magnitude events since Central and Eastern Iran are characterized by earthquakes with long recurrence times. In contrast, the averaged PGA values for 2% POE in 50 years from EMME are 12% and 5% higher than PGA values derived from 1 3  DBHM, respectively. In the EMME model, CIB in Central Iran province and the LDB in Eastern Iran province are characterized as stable continental regions. The main reason for the contrast is GMPEs (i.e., Atkinson and Boore (2006), Campbell (2003), and Toro (2002) used for these zones in the EMME model. These GMPEs are developed for stable continental regions and provide higher PGA value for large-magnitude earthquakes in comparison to our considered GMPEs.

Hazard curves
In order to have a full picture of the annual probabilities of exceedance, we depict PGA hazard curves for important cities in Iran by picking up the POEs for 16 hazard levels between 0.001 and 3.5 g in 50 years (Fig. 10). These cities are located in different seismotectonic provinces as follows: Tehran in the Alborz province, Tabriz in the Azerbaijan province, Kermanshah, Ahvaz, Shiraz, and Bandar Abbas in the Zagros province, Kerman in Eastern Iran, Isfahan in Central Iran and Bojnurd and Mashhad in the Kopeh Dagh province. In the following, we made a comparison between the DBHM hazard values and the results of available seismic hazard models for different regions of the Iranian Plateau.
Two catalog-based seismic hazard models for Iran are developed using the smoothed seismicity approach. The first one is Khodaverdian et al. (2016b) (Kh16), which covered the whole Iranian Plateau and the second one is Khoshnevis et al. (2017) (Kh17), which provided for the northern part of Iran including Azerbaijan, Alborz, Kopeh Dagh, the northern part of Zagros, and the northern part of the Central Iran and Eastern Iran provinces. In both models, all the observed earthquakes are distributed in the area sources and are not assigned to specific faults. Hence, the effect of fault proximity is not seen in their hazard assessment results. Moreover, two other models are available for Tehran and surrounding areas: (1) Zafarani et al. (2017) (Za17), in which seismic hazard model is provided based on the characteristic earthquake model, and (2) Jalalalhosseini et al. (2018) (Ja18), providing a time-dependent hazard model. Shabani and Mirzaei (2007) (Sh07) is a catalog-based PSHA model for the Kermanshah region. Eskandarinejad et al. (2018) (Es18) include both traditional catalog-based and Monte-Carlo simulation-based PSHA for Shiraz city. The 16%, 50%, mean, and 84% quantile curves and results of the above-mentioned models for selected cities are shown in Fig. 10.
The results of the mentioned models for Tehran are shown in Fig. 10a. The PGA values for 10% and 2% POEs in 50 years, derived from DBHM are 0.37 g and 0.62 g, respectively.
The PGA values for 10% and 2% POEs, derived from EMME are 0.45 g and 0.87 g, respectively. The differences between PGA values derived from DBHM and EMME for Tehran are lower than the differences we have seen for the large number of sites in the Alborz province. The main reason is that for the region close to Tehran, the geological properties (e.g., geological slip rates) of the faults near the city are studied in various researches (Ritz et al. 2003;Nazari et al. 2009) and implemented in EMME. In general, we may expect the hazard results are in good agreement with our prediction if a geological investigation of the faults near a site is implemented in the hazard model. Estimates from Kh16 (i.e., 0.22 g and 0.45 g for 10% and 2% POEs, respectively), Za17 (i.e., 0.28 g and 0.54 g), Ja18 (i.e., 0.29 g and 0.56 g), Kh17 (i.e., 0.31 g and 0.53 g), in which only the catalog of observed earthquakes are used, are less than the estimate from DBHM. That could be justified by the presence of well-known active faults nearby Tehran (e.g., North Tehran and Mosha faults) while the observed earthquakes are assigned to a wide area instead of near-fault zones, and consequently smoothed seismicity hazard models (i.e., Kh16 and Kh17) resulted in lower PGA values.
Estimated PGA values for Tabriz city (i.e., 0.55 g and 0.87 g for 10% and 2% POEs in 50 years, respectively) are close to the results of the EMME model (i.e., 0.61 g and 1.10 g) (Fig. 10b). The active faults nearby Tabriz city are studied in various researches (Rizza et al. 2013) and the geological slip rates are implemented in the EMME project model, which is the main reason for this consistency. Estimates of Kh16 (i.e., 0.22 g and 0.46 g) and Kh17 (i.e., 0.35 g and 0.61 g) are less than the hazard values from DBHM and EMME models. The reason is the lack of a fault-based approach in those models; for instance, the North Tabriz fault, which is located in the proximity of the city, is not directly considered in both Kh16 and Kh17 models.
For Ahvaz city in the Zagros province, the PGA values for 10% and 2% POEs in 50 years derived from the DBHM model are 0.23 g and 0.41 g, respectively, and are in good agreement with Kh16 estimated values (i.e., 0.21 g and 0.41 g) (Fig. 10c). It shows that distributed seismicity is a plausible assumption for that region. This fact could be supported by the presence of several hidden secondary/minor faults located around the city of Ahvaz (Berberian 1995;Abedi and Oskooi 2015;Sabour et al. 2015). For Shiraz, the estimate of Es18 is ~ 0.22 g for 10% POE, which is close to the DBHM's result (i.e., 0.26 g) (Fig. 10f). For Bandar Abbas, the PGA values for 10% and 2% POEs in 50 years derived from the DBHM model are 0.36 g and 0.63 g, respectively, and are in good agreement with Kh16 estimated values (i.e., 0.35 g and 0.64 g) (Fig. 10e). The PGA values derived from deformation-based and catalog-based models for most selected cities located in the Zagros seismotectonic province (e.g., Ahvaz and Bandar Abbas) are comparable, which could be explained by the relatively short recurrence times of earthquakes that occur in the Zagros. Comparing our model with the EMME, we have seen that our results for all selected cities located in the Zagros (i.e., Ahvaz, Kermanshah, Bandar Abbas, and Shiraz) are lower than EMME's estimates ( Fig. 10c-f), which is supporting the pattern we mentioned and justified in Sect. 6.1.
Bojnurd and Mashhad are two selected cities in the Kopeh Dagh seismotectonic province. The estimated PGA values for Bojnurd (i.e., 0.38 g and 0.66 g for 10% and 2% POEs in 50 years, respectively) are close to the values derived from the EMME (i.e., 0.36 g and 0.69 g) and more than from the PGA values of Kh16 (i.e., 0.27 g and 0.54 g) (Fig. 10g). Moreover, for Mashhad, the estimated PGA from the DBHM (i.e., 0.35 g and 0.62 g) are more than the values derived from the EMME model (i.e., 0.21 g and 0.44 g) and Kh16 (i.e., 0.15 g and 0.34 g) (Fig. 10.h). As mentioned in Sect. 6.1, high activity rates based on strain rates have been seen for the strip, which is extended from the east of the Alborz to Mashhad city. Kh17 estimated 0.22 g, and 0.39 g for 10% and 2% POEs, respectively, which are close to PGA values from the EMME model, as both models are based on observed seismicity.
Our PGA estimate for 10% POE in 50 years for Isfahan city in Central Iran (i.e., 0.15 g) is quite close to the PGA values derived from EMME (i.e., 0.15 g) and Kh16's estimate (i.e., 0.14 g) (Fig. 10i). For 2% POE in 50 years, our PGA value (i.e., 0.27 g) is comparable to the value derived from Kh16 (i.e., 0.32 g); however, our estimate is significantly lower than the EMME's estimate (i.e., 0.43 g) (Fig. 10i). In the EMME model, Isfahan is mostly affected by a zone considered as a stable continental region.

3
GMPEs which have been used for that, provide a higher PGA value for large-magnitude earthquakes in comparison to our considered GMPEs. Figure 10.j shows that our PGA estimates with 10% and 2% POEs in 50 years for Kerman city in Eastern Iran are 0.31 g and 0.56 g, respectively, and more than the estimates of EMME (i.e., 0.25 g and 0.47 g). Kh16's estimates (i.e., 0.26 g and 0.49 g) are close to the estimates of EMME. Earthquakes in Eastern Iran seismotectonic province almost occur with long recurrence intervals; therefore, catalog-based seismic hazard might underestimate the possible expected PGA.
We also build another version of the DBHM model by employing GMPEs from the EMME model to investigate the effect of GMPEs. Hazard analyses are carried out for the cities where significant discrepancies in the outcomes of the DBHM and EMME models have been observed. For cities in regions that are considered as the active shallow crust (e.g., Mashhad (Fig. 11a) and Kerman (Fig. 11b)) in both models, the application of the EMME's GMPE does not lead to substantial variations in hazard curves and the estimated PGA for 475 years return period is changed less than 20%; it can be interpreted that a majority of the difference (more than 80% for 475 years return period) between DBHM and EMME models is due to differences in earthquake occurrence rates. For cities in regions that are considered as stable continental regions (e.g., Isfahan (Fig. 11c)) in the EMME, the results are close to the EMME model, especially for 2475 years return period. Mean hazard curves from DBHM, EMME, and modified DBHM with EMME's GMPEs

Conclusions
Several investigations have shown that available short instrumental earthquake catalogs or limited historical evidence may not be enough for estimation of the long-term earthquake occurrences rate, especially for an area that is characterized by large-magnitude earthquakes. Catalog-based seismic hazard assessment, consequently, could be challenging. Deformation data (i.e., GPS velocities and geological fault slip rates) could be used for the enrichment of a regional seismic hazard model when it is combined with seismic data from earthquake catalogs. We used a methodology to estimate earthquake occurrence rates from strain rates and implemented it by using the outcome of the comprehensive deformation model of the Iranian Plateau. The comparison between the catalog-based long-term earthquake occurrence rates and deformation-based occurrence rates shows the lowest amount of D 4 ∕ C 4 (i.e., 1.34) is observed in the Zagros province, where earthquakes with relatively short recurrence times occur. The maximum amount of D 4 ∕ C 4 (i.e., 2.35) is observed for Central Iran, where large-magnitude earthquakes occur with long recurrence times. Sensitivity analysis shows that the variation of mean coupled seismogenic layer thickness has a considerable effect on driving activity rate from strain rate, and it could not be ignored compared with the uncertainties of the seismicity parameters.
We generated square area sources, for a grid with the same spacing as the strain rate model (0.2 • × 0.2 • ) covering the whole Plateau. The seismicity parameters (e.g., deformation-based activity rates, and m max ), evaluated based on the proposed methodology, are assigned to each 0.2 • × 0.2 • area source. The seismotectonic characteristics of each area source are evaluated based on the characteristics of the nearest fault and the seismotectonic map of the Iranian Plateau. Logic trees of the hypocenter depths, seismicity parameters, and GMPEs are employed in an attempt to capture their uncertainties.
We provided PGA hazard maps at 10% and 2% POEs in 50 years. The highest levels of hazard are found in the Azerbaijan and Alborz seismic provinces, where the highest value of strain rates is also located. Generally, the comparison between the DBHM and catalogbased models shows that catalog-based PGA values are lower than deformation-based PGA values for most selected cities in the Iranian Plateau, which could be justified by the lack of recorded large-magnitude earthquakes due to their long-term recurrence times. Specifically, the comparison between the EMME and DBHM can be summarised in four points: (1) For 10% POE in 50 years, there are up to ~ 23% differences between the spatial average of PGA for provinces. The DBHM resulted in higher PGA values than estimates from EMME for all provinces except the Zagros. In the EMME model, GMPEs, which resulted in high PGA, are used for Zagros. (2) For 2% POE in 50 years, there are up to ~ 28% differences between the spatial average of PGA values. The averaged hazard values from DBHM are higher than the corresponding from EMME for all provinces except the Zagros, Central and Eastern Iran. In the EMME model, CIB in Central Iran province and the LDB in Eastern Iran province are characterized as stable continental regions, and the considered GMPEs for these regions provide higher PGA value for large-magnitude earthquakes in comparison to our considered GMPEs. (3) For the city of Tehran and Tabriz, the results of the EMME and DBHM models are in good agreement. The geological properties (i.e., slip rates) of active faults near these cities are implemented in the EMME model and the consequences of the lack of a long-term catalog have been reduced. Therefore, we expect catalog-based hazard models have a better estimate if the geological properties of the active faults near a site are implemented in a catalog-based model. (4) By employing GMPEs from the EMME model in another version of the DBHM model, we found that a majority of the difference between DBHM and EMME models is due to differences in earthquake occurrence rates in the active shallow crust regions. However, it should be mentioned that for areas, which are considered as stable continental regions in EMME, the differences between DBHM and EMME models are mostly resulted from different GMPE sets used.