Macrozonation of seismic transient and permanent ground deformation of Iran

Iran is located on the Alpide earthquake belt, in the active collision zone between the Eurasian and Arabian plates. This issue makes Iran a country that suffers from geotechnical seismic hazards associated with frequent destructive earthquakes. Also, according to the rapid growth of population and demands for construction lifelines, risk assessment studies which should be carried out in order to reduce the probable damages are necessary. The most important destructive effects of earthquakes on lifelines are transient and permanent ground displacements. The availability of the map of the displacements caused by liquefaction, landslide, and surface fault rupture can be a useful reference for researchers and engineers who want to carry out a risk assessment project for each specific region of the country. In this study, these precise maps are produced and presented by using a considerable number of GIS-based analyses and by employing the HAZUS methodology. It is important to note that a required accuracy for risk assessment is approximately around the macro scale. So, in order to produce a suitable map for risk assessment goals, in terms of accuracy, the GISbased analyses are employed to map all of Iran.


Introduction
As the plate tectonic configuration of Iran in Fig. 1 shows, the country, which is located on the Alpide earthquake belt, is one of the most highly earthquake-prone zones of the world (Taherian and Kalantari, 2019). In addition, Iran has various geographical conditions that can cause landslide and liquefaction, in mountainous and coastal areas, respectively. The first earthquake effect, which can damage lifelines and infrastructure, is the transient ground displacement (TGD), which is caused by seismic wave propagation. The second one is the permanent ground displacement (PGD), which may result in liquefaction, landslide, and ground failure. For risk assessment of lifelines and infrastructure which are widespread throughout the country, investigating the TGD and PGD is of vital importance. Many studies have proposed technical methods for evaluating TGD and PGD and for specific cases in different regions of the country, some of which are discussed in the following paragraphs.
While landslides are considered as one of the most disastrous natural hazards in Iran, there is a lack of precise information on them for most parts of the country, and only a small percentage of the country's area has been specifically investigated to provide data for landslide susceptibility maps. Tangestani (2004) investigated the landslide susceptibility mapping using the fuzzy gamma approach in a GIS basis for the Kakan catchment area, south-west Iran. Babakan et al. (2009) proposed a seismo-geotechnical zonation mapping of the southern Caspian Sea coastline. Daneshvar and Bagherzadeh (2011) evaluated the landslide hazard zonation using GIS analysis at Golmakan Watershed, north-east of Iran. Moradi et al. (2012) implemented a GIS-based landslide susceptibility mapping by employing the analytical hierarchy process (AHP) method for Dena City. A landslide hazard zonation was carried out by employing statisticalbased methods for Pishkuh region in Fereydonshahr by Shirani and Seif (2012). Aghda and Bagheri (2015) evaluated an earthquake-induced landslide hazard zonation method for the Sarein earthquake in 1997. A landslide hazard zonation and risk analysis in Goloord region, north of Iran, was carried out using the AHP method by Adib and Afzal (2018). Arjmandzadeh et al. (2019) presented a GIS-based landslide susceptibility mapping for Qazvin Province in Iran. Mokhtari  (Alavi, 1994). and Abedian (2019) investigated the spatial prediction of landslide susceptibility in the Taleghan basin. Vakhshoori et al. (2019) studied the landslide susceptibility mapping of Bandar Torkaman by employing GIS-based data mining algorithms.
There have also been investigations on landslides using remote sensing tools. Esmali and Ahmadi (2003) evaluated a mass movement hazard zonation using GIS and remote sensing (RS) in Germichay Watershed, Ardebil. A monitoring of the massive slow Kahrod landslide in the Alborz range was implemented using GPS and synthetic aperture radar interferometry by Peyret et al. (2008). Akbarimehr et al. (2013) assessed the slope stability of the Sarcheshmeh landslide, north-east Iran, by using interferometric synthetic aperture radar (InSAR) and GPS observations. Mirzaee et al. (2017) evaluated three InSAR time-series methods to assess the creep motion of the Masouleh landslide in north Iran. Pirasteh et al. (2018) used a lidar-derived digital elevation model (DEM) and a stream length-gradient index approach for investigating the landslides in the Zagros Mountains. A landslide hazard mapping using a radial basis function neural network model was performed for a case study in Semirom, Isfahan, by Yavari et al. (2019).
From a different view, liquefaction is also one of the seismic geohazards which can significantly affect the performance of lifelines during or after earthquakes. Studies have addressed liquefaction through different methods for different regions of Iran. Askari et al. (2006) evaluated the liquefaction potential of the south of Tehran using the standard penetration test and the shear wave velocity measurement. Naghizadehrokni et al. (2018) presented liquefaction maps in Babol City using probabilistic-and deterministic-based approaches. Risk assessment of existing structures due to the liquefaction potential of Astaneh-ye Ashrafiyeh City was performed by Ziabari et al. (2017). Liquefaction assessment using micro-tremor measurements and artificial neural networks was carried out by Rezaei and Choobbasti (2014) for Babol City. Sakvand et al. (2011) investigated liquefaction risk zoning in the Silakhor plain. Liquefaction-induced lateral spreading displacement was evaluated probabilistically for a site in the south of Iran by Kavand and Haeri (2009). Koike et al. (2004), Mousavi et al. (2014), andFarahani et al. (2020) also evaluated liquefaction-induced displacement of Tehran, Azerbaijan, and Asaluyeh, respectively, in order to assess the risk of the gas pipelines.
The majority of large earthquakes are associated with surface ruptures, which cause secondary hazards to arise. Fault rupture hazard is defined as a displacement imposed by fault rupture on structures and objects during an earthquake (Perrin and Wood, 2003). There are empirical equations which are established based on the global and regional records of seismic events and are used to predict geometrical and kinematic characteristics of the potential ruptures along active faults, including surface rupture length (SRL), maximum displacement (MD), and average displacement (AD) (e.g. Öztürk et al., 2018;Manighetti et al., 2007;Dowrick and Rhoades, 2004;Mason, 1996;Wells and Coppersmith, 1994). SRL and MD are correlated with each other and earthquake magnitudes and provide the most well-known equations for deterministic evaluation of earthquake hazards imposed by faults as significant sources of seismic energy. Stramondo et al. (2005) investigated the surface displacements and source parameters of the 2003 Bam earthquake using Envisat advanced synthetic aperture radar imagery. Surface displacement and fault modelling for the 2003 Bam earthquake were evaluated using the InSAR method by Stramondo et al. (2005).
However, there are few studies that have addressed all the ground displacements caused by earthquakes for all the regions of Iran. Moreover, there is no comprehensive study presenting a map of surface-rupture-induced displacement of Iran. Some studies proposed only empirical relations between different parameters of Iran's faults. However, these parameters have never been calculated for all of Iran's faults in order to estimate the rupture-induced displacements in a widespread zone of the country. In this study, PGD is calculated and mapped using the HAZUS methodology (FEMA, 2012). Also, a map of ground displacement due to surface rupture is produced via a GIS-based approach and the HAZUS methodology. Hence, the novelty of this study is not only the macro zonation of the PGD caused by earthquakes all over Iran, but also the presentation of the first map of fault displacement, which can affect the lifelines it is near to or crosses. As well, all mapping of displacements is carried out on a macro scale. This is due to the fact that from a risk assessment perspective, macro zonation is useful enough and that there is no need to study the issues using a micro-scale approach. Therefore, the HAZUS methodology is employed here in order to take advantage of its straightforward equations and fragility curves, which were obtained by a huge number of analytical and experimental studies worldwide. It is important to note that although the HAZUS methodology can be used as a functional tool for loss estimation in almost all countries, it is rare to find a comprehensive study to present macrozonation of seismic ground deformation of any country around the world. The majority of the previous studies were carried out for specific small regions in the countries. The exhaustive maps presented in the current study can facilitate future seismic risk assessment projects, which will study any lifeline systems in all areas of the country. Hence, the authors of this study believe that such a comprehensive endeavour to present practical maps can be useful for other countries which suffer from seismic hazards.
As we know, a rough estimation of the severity and geographic distribution of lifeline seismic risks can help risk managers to allocate critical resources in disaster management, mitigation, and preparedness processes. Since investigating and quantifying the seismic hazards is the first step to evaluating the seismic-induced damages to any lifelines, the output of the current study, in the form of PGD and TGD maps, can be considered as one of the most important input data for the lifeline seismic risk assessment projects. As shown in Fig. 2, in order to produce a peak ground velocity (PGV) map, a classified soil map was produced using the USGS ShakeMap method and by employing an available slope map of the country. According to the Iranian seismic Table 1. Correlations between topographic gradient and VS30 using the NED 9c digital elevation models for the National Earthquake Hazard Reduction Program (NEHRP) site classes (Allen and Wald, 2009). NEHPR V S30 range 9 arsec gradient range 9 arsec gradient range Modified 30 arsec site code, a reflection factor map and after that a spectral acceleration map was produced in order to generate the PGV map of the country. In addition, an available liquefaction and landslide susceptibility map was employed to produce the PGD maps using HAZUS methodology and GIS-based analyses. Finally, the fault maximum displacement map is produced using a fault maximum magnitude map, which was presented previously by Karimiparidari (2014).

Hazard analysis of ground shaking
For estimating the TGD caused by seismic waves propagation (ground shaking), PGV is needed. As HAZUS proposed, for obtaining PGV, the first step is to calculate the spectral acceleration by having a soil classification of a region in terms of dynamic properties. According to the ShakeMap (Wald et al., 2005) method, for regions lacking Vs30 maps, including most of the globe, the approach of Allen and Wald (2007), which was revised by Allen and Wald (2009) and provides estimations of V s30 as a function of more available topographic slope data, can be employed. It is important to note that some validation studies were performed by researchers like Shahvar (2013), in order to evaluate the merit of the Allen and Wald (2009) study for the geological situation of Iran. In this study, soil classification is carried out using a topographic gradient map. As shown in Fig. 3a, a global 1 arcsec (30 m) SRTM DEM of Iran is used for producing a slope map (stated in Fig. 2(1a)) as shown in Fig. 3b. After that, the soil classification map (stated in Fig. 2(1c)) is produced as shown in Fig. 4 and using Table 1, which presents correlations between topographic gradient and V S30 . According to the Iranian seismic code (also known as the Standard No. 2800) (BHRC, 2015), for calculating spectral acceleration, a reflection factor should be obtained. Reflection factor (known as B factor) is considered to account for the resonating effect of soft soil on ground movement at bedrock level; its value increases as the soil gets softer. The value of the reflection factor is relevant to two main param-  eters consisting of B1, spectrum shape factor, and N, spectrum modification factor. The parameters mentioned are correlated to the soil type and level of seismicity. According to the Iranian seismic code, Iran is divided into four seismic zones, including low, moderate, high, and very high seismicity levels. Also, the soil types consisting of types B-E are presented for the country. Hence, by merging the zonation of seismicity level and the soil classification map, the soil and seismic hazard class map is produced as shown in Fig. 5.
The value of B is obtained in eight different combinations of soil type and seismicity level by using the reflection factor spectrum (see Fig. 6) in order to calculate the PGV inferred from 1 s spectral response. The results are shown in Table 2. Therefore, the map of the reflection factor for the 1 s period (stated in Fig. 2(1e)) is obtained, as shown in Fig. 7a. Finally, by multiplying the reflection factor to a peak ground acceleration (PGA) map, the 1 s spectral acceleration (stated in Fig. 2(1f)) is produced as shown in Fig. 7b. PGV is inferred from 1 s spectral acceleration using Eq. (1).
The constant value of 1.65 in Eq.
(1) represents the amplification assumed to exist between peak spectral response (1 s) and PGV. This value is based on the median spectrum amplification, as given in Newmark (1982), for a 5 % damped system whose period is within the velocity-domain region of the response spectrum. A PGV map (stated in Fig. 2(1g)) of Iran is presented in Fig. 8.

Hazard analysis of ground failure
The ground failure is divided into the three main categories: liquefaction, landslide, and faulting. Each of these types of ground failure is quantified by PGD. Methods and alternatives for determining PGD due to each mode of ground failure are discussed below.

Liquefaction
Liquefaction is the most important hazard due to ground failure that often threatens infrastructures. Liquefaction is a soil behaviour phenomenon in which a saturated soil loses a substantial amount of strength due to high excess pore-water pressure generated by and accumulated during strong earthquake ground shaking (FEMA, 2012). In this study, in order to consider the failure caused by soil liquefaction, the Iran liquefaction susceptibility map (stated in Fig. 2(2a)) is used. This map is provided by the International Institute of Earthquake Engineering and Seismology (IIEES) and based on previous studies by Komakpanah and Farajzadeh (1995), as shown in Fig. 9a. The likelihood of experiencing liquefaction at a specific location is primarily influenced by the susceptibility of the soil, amplitude, duration of ground shaking,  and depth of groundwater. Based on the HAZUS methodology, the probability of liquefaction for a given susceptibility category can be determined using Eq. (2): where P [Liquefaction|PGA = pga] is the conditional liquefaction probability for a given susceptibility category at a specified level of PGA, K M is the correction factor for moment magnitudes other than M = 7.5 as presented in Eq. (3), K W is the groundwater correction factor, and P ml is the proportion of the map unit susceptible. In terms of the values of K W and K M factors, the k W parameter is ignored due to the lack of a groundwater-level map of the coun-try. However, the K M factors are calculated using the moment magnitudes of seismic provinces of the country, which was presented by Karimiparidari (2014). The values of P ml and P [Liquefaction|PGA − pga] are assigned according to HAZUS recommendations for each susceptibility class, as presented in Table 3. Zonation of the probability of liquefaction for all susceptibility categories is carried out, as shown in Fig. 9b-d.
The expected value of PGD conditioned to the occurrence of liquefaction can be stated as a function of PGA (Sadigh et al., 1986), as presented in Eq. (4).
where PGA(t), which is presented in Table 4, is the threshold ground acceleration corresponding to zero probability of liquefaction, and K is the displacement correction factor given by Eq. (5). Mapping of the threshold ground acceleration (stated in Fig. 2(2b)) is shown in Fig. 10. As a final result, Fig. 11 presents the liquefaction-induced displacement map (stated in Fig. 2(2e)) of Iran.

Landslide
Earthquake-induced landslide of a hillside slope occurs when the static plus inertia forces within the slide mass cause the factor of safety to drop below 1.0 temporarily (FEMA, 2012). The value of the PGA within the slide mass required to cause the factor of safety to drop to 1.0 is denoted by the critical or yield acceleration (a c ). This value of acceleration is determined based on pseudo-static slope stability analyses and/or empirically based on observations of slope behaviour during past earthquakes. The landslide hazard evaluation requires the characterization of the landslide susceptibility of the soil and geologic conditions of a region or sub-region. For this purpose, the Iran landslide susceptibility map (stated in Fig. 2(3a)), provided by Geological Survey and Mineral Explorations of Iran (GSI, 2018), is used as shown in Fig. 12. Also, critical acceleration (stated in Fig. 2(3b)) at any location proposed by HAZUS for susceptibility categories is presented in Table 5 and Fig. 13.
The permanent ground displacements are determined using Eq. (6): where E[d/a is ] is the expected displacement factor, a is is the induced acceleration (in a decimal fraction of g's), and n is the number of cycles of ground shaking. A relation derived from the results of Makdisi and Seed (1978) is used to calculate downslope displacements. In this relation, shown in Fig. 14, the displacement factor d/a is is calculated as a function of the ratio a c /a is using the upper bound values, in order to be conservative. Also, Eq. (7), which represents the relationship between the number of cycles and earthquake moment magnitude based on Seed and Idriss (1982), is used for calculating the number of cycles of ground shaking (n). As a result, maps of the ratios of critical to induced acceleration (a c /a is ) (stated in Fig. 2(3d)) and displacement factor (d/a is ) (stated in Fig. 2(3e)) are presented in Fig. 15. Finally, the zonation of landslide-induced displacement (stated in Fig. 2(3f)) is carried out using GIS-based analyses and presented in Fig. 16.   Table 5. Critical acceleration at any location proposed by HAZUS for susceptibility categories. None  I  II  III  IV  V  VI  VII  VIII  IX

Surface fault rupture
Active faulting in Iran is a direct indicator of active crustal deformation due to the convergence between Arabia and Eurasia, which occurs at 2.1-2.5 cm yr −1 . During the last 500 years surface ruptures associated with large earthquakes have appeared or been documented in various places in Iran.
Most of these ruptures occurred along the active faults which moved repeatedly in the Quaternary period, thus constituting evidence that these active faults have the potential of reactivating in the future (Hessami and Jamali, 2006).  The most recent seismic hazard map of Iran was developed by Karimiparidari (2014) using the available data and based on PSHA approach. This covers a wide time span of earthquakes history and contains uniform scaled magnitudes. Karimiparidari (2014) has also developed new seismic source models and seismotectonic zoning maps of Iran. The seismotectonic models were developed based on the latest data of active tectonics, topography, magnetic intensity, and seismicity catalogue. These new maps divide the country into 27 seismotectonic zones and demonstrate two models for linear and regional seismic sources. As shown in Fig. 17, seismicity parameters of 104 seismic regions, presented in 27 seismotectonic zones, are assigned to the faults. The mentioned parameters are considered to estimate the most prob-  able maximum magnitude of each fault in order to calculate the rupture-induced displacement.
By using the database of the surface ruptures of Iran, empirical relations are established for moment magnitude and maximum displacement (MD), as given in Table 6. Coefficients of the relations are separately calculated for the thrust, strike-slip faults, and all of the fault types. It is worth noting that active normal faults are rare in Iran, and surface ruptures associated with this kind of earthquake faulting are even more scarce (Ghassemi, 2016). As a result of the surface fault rupture study and using the empirical equation (Eq. 8) the map of surface-rupture-induced displacement (stated in Fig. 2(4b)) is produced by employing GIS-based analyses as presented in Fig. 18.  where M w is moment magnitude. Also, the regression coefficients used, a and b, are presented in Table 6.

Conclusion
Being located in the active collision zone between the Eurasian and Arabian plates, Iran is a country that suf- fers from hazards associated with frequent destructive earthquakes. The susceptibility assessment of infrastructures is crucial in the modern era due to the very rapid growth of population and major cities, which are mostly located on or in the vicinity of earthquake faults, and also demands the construction of infrastructure that is susceptible to earthquake hazards. The geotechnical seismic hazard which can affect the serviceability of lifelines during or after earthquakes can be classified into two categories: transient ground displacement (TGD) caused by seismic wave propagation (ground shaking) and permanent ground displacement (PGD), which refers to liquefaction, landslide, and surface fault rupture. There are many theoretical, experimental, and numerical methods for evaluating earthquake-induced displacements, which can affect lifelines significantly. For example, in order to investigate the landslide and liquefaction potential of a specific limited region, geotechnical-based field experimental studies and finite-element-based methods can be implemented. However, from a risk assessment point of view, empirical-theoretical-based methods are even more useful for macro-scale regions. This is because the number of required parameters for empirical equations is lower than the number of parameters which are required for numerical analyses. Hence, from a risk assessment point of view, the zonation of earthquake-induced displacements can help researchers and engineers to carry out their research more rapidly by using the prepared map of displacements in the country. Therefore, the main goal of this paper is to produce and present maps of earthquake-induced displacements.
For creating the precise maps, GIS-based analyses were carried out by employing the HAZUS methodology. A peak ground velocity (PGV) map of Iran was produced using soil classification estimation based on topographical data, spectral acceleration calculation, and the HAZUS equations. Although the PGV can be obtained using attenuation relationships, the proposed method by HAZUS is selected for being employed in this study. Investigating the liquefaction-induced displacements, the probability of liquefaction for each susceptibility category was calculated using the HAZUS equations, and a map capable of presenting the most probable displacements was produced. GIS-based analyses, Makdisi and Seed's equation, and a landslide susceptibility map were used for preparing the landslide-induced displacement maps. Also, a seismotectonic zoning map was employed to estimate the most probable maximum magnitude of each fault and to evaluate the surface fault rupture based on displacement. A map of the surface-rupture-induced displacements was also produced.
In this study, there are some limitations which the authors faced. The first one is the accuracy of the available DEM of the country. As was discussed, the accuracy of the used DEM is around 1 arcsec (30 m) and can affect the produced PGV map of the country. The other limitation is the Iran liquefaction susceptibility map, which is respectfully old fashioned (1996). The Iran liquefaction susceptibility map should be updated periodically because the level of the groundwater has continuously varied in recent decades due to the severe climate changes. Consequently, having a more accurate DEM and employing up-to-date liquefaction susceptibility zonation can help produce a cutting-edge version of the result of this research in the future.
Appendix A: Reflection factor calculation based on the Iranian seismic code (BHRC, 2015) B is the reflection factor that is obtained from a smoothedelastic design response spectrum in the following form: where B 1 is the spectral shape factor, which is calculated by Eqs. (A2) to (A4), and N is the spectral modification factor, which is evaluated by equations presented in Table A3.
In these equations, T is the structural fundamental period of vibration (in s); T 0 , T s , and S are parameters related to the site soil conditions and seismic potential of the region as given in Table A2.    Table A3. N (spectral modification factor) calculation equations (BHRC, 2015).
High and very high Low and medium relative hazard area relative hazard area