# Seismic Hazard Evaluation of the North and Northwest Iranian Plate Considering the Bayesian Perspective

Authors: Mehdi Maleki; Zohreh Sadat Riazi Rad
Journal: Shock and Vibration
 (2022)

---

## Abstract

In the last hundred years alone, there have occurred more than 2000 earthquakes larger than 4 in the northern and northwestern parts of Iran. Because Iran is placed in the heart of the Alpine-Himalayan collision zone, which is one of the world’s greatest seismic locations, massive earthquakes have regularly struck the country. Seismic hazard assessment is a realistic and helpful method for engineering seismology to prevent motions in a location in order to forecast the occurrence of earthquakes. For parameter estimation in seismic hazard assessments, the Bayesian probability theory provides a practical and accurate technique. As a result, Bayesian techniques for analysis are based on an older method. Furthermore, they have the unique capacity to account for the variability of characteristics in probabilistic relationships. In this assessment, we used the Bayesian earthquake risk estimation program developed by Lyubushin. Our research is in the sequence of applications of this program to assess the risk of earthquakes in various parts of the universe. In this method, the basic assumption is that earthquakes follow the characteristic of being Poisson. Its zoning map is drawn in the north and northwest regions of Iran. The highest value was obtained for the city of Tabriz with 0.29 galls, and the lowest value was obtained for the city of Zanjan with the value of 0.06 galls.


---

## Body

## 1. Introduction

Earthquake is one of the most important natural events that have always harmed humans throughout history [1]; hence, earthquake-resistant constructions are unavoidable and required. The most important of them are the earthquakes of June 20, 1990, in Rudbar with a magnitude of 7.2, the earthquake of May 28, 2008, in Kojoor with a magnitude of 6.3, the earthquake of August 11, 2012, in Varzeqan with a magnitude of 6.3, the earthquake of May 8, 2020, in Tehran with a magnitude of 4.5, and the earthquake of October 24, 2020, in Qazvin with magnitude 5.4 as pointed out. The considered area is one of the most active areas with the probability of earthquakes with a magnitude greater than 6. Because of the strong seismicity of Alborz, Azerbaijan, experts are continually looking for new approaches to give earthquake risk analyses. For this purpose, it is necessary to conduct various studies in a densely populated area of economic and political importance in the north and northwest regions of Iran.The application of a Bayesian approach was estimated of seismic hazard parameters in some regions of the circum-pacific belt and local intensity for some cities in Japan [2, 3]. The seismic hazard for a number of locations in Greece was calculated using a Bayesian estimate of peak acceleration [4]. Lyubushin and Parvez created a map of India’s seismic danger using a Bayesian technique [5]. The Yadan’s probabilistic assessment of earthquake hazard parameters was derived using a Bayesian technique in the Northwest Himalayas. In Iran, earthquake risks were estimated using Bayesian methods [6]. Hamzehloo et al. [7] calculated seismic hazard zoning for a 475-year return time. In this region, Yazdani and Kowsari used a nonparametric technique to evaluate seismic danger [8]. In the Iranian Plateau, earthquake hazard indicators were assessed [9]. For diverse places of the Iranian Plateau, a quantitative assessment of earthquake hazard characteristics using a Bayesian technique was carried out [10]. The modified Omori law parameters for the Iranian Plateau were determined using Bayesian methods [11]. Salamat et al. estimated confidence ranges for earthquake maximum magnitudes in several seismotectonic zones in Iran [12]. Salamat and Zare developed Bayesian and frequentist methods for estimating the maximum predicted earthquake magnitude in Iran [13].This study was conducted to obtain the seismic hazard map and peak ground acceleration (PGA) at a 90% probability level with the Bayesian approach for intervals of 50, 100, and 475 years. Considering each earthquake as a potential source has been done for the first time in the north and northwest regions of Iran Figure1 One of the most significant benefits of Bayesian techniques is that they offer users with a strong assumption. In the probabilistic connection, each parameter’s ambiguity, as well as the prognosis, must be taken into consideration [14, 15]. The Bayesian statistical method is very important in seismic hazard analysis for two reasons. [16] First, the Bayesian theory is an accurate way to count seismic prior information. This information is subjective, geological, or statistical [17], combined with historical observations of earthquakes [18]. The method’s second aspect is that it is used to incorporate a statistical uncertainty associated with the estimate of seismicity-related parameters [19]. Possible uncertainty is also related to the inherent uncertainty of the earthquake [20]. One of the advantages of Bayesian statistics over the classic is that in the conventional probabilistic method, a constant value for PGA is provided in different return periods. But in the Bayesian method, PGA values are expressed distributive at different future intervals [10]. Their proximity to large cities and the lack of detailed information on their seismogenic potential call for multiparametric research. Since 2013, the area’s crustal deformation has been monitored by a dense GNSS network (PPGNet), consisting of five stations, equipped with Leica and Septentrio receivers [21]. The group velocities are found from the application of the multiple filter technique in a single-station fashion, while for the phase velocities, slant stacking or linear radon transform are applied in the fashion of multichannel analysis of surface waves (MASW) [22]. Soil’s potential against liquefaction for 45 locations has been carried out using PGASUR results so obtained.Figure 1
The seismicity map of Alborz-Azerbaijan seismic province and seismic locations greater thanMw > 3.In this study, the use of a personal (optional) combination of seismic parameters in seismic hazard assessments was studied utilising the Bayesian method’s described properties. Bayesian also establishes the requirements for accounting for uncertainty in computations. Pisarenko et al. [23] and Pisarenko and Lyubushin [24] describe the Bayesian technique. Lyubushin and Parvez [5] later tweaked maps based on Bayesian estimations of peak acceleration statistics. The method’s core computational code, developed by Lyubushin, is used to predict earthquake risk in various locations of the world [4, 10, 26–31]. These maps show the acceleration in different time periods and the maximum values of the Earth’s acceleration in the cities of Tehran, Tabriz, Rasht, and Zanjan, respectively.


## 2. Method of Estimate

In research, the approach employed is detailed in depth [2, 4, 23, 24]. As a result, we will just present the basic hypotheses and crucial equations in this section. Let R be the decimal logarithm of the maximum seismic peak ground acceleration at a certain site, computed as a series through using the attenuation law during a preceding time frame (−τ, 0)(1)⟶Rn=R1,…,Rn,Ri≥R0,Rτ=max1≤i≤nR1,…,Rn.where R0 is a minimum cutoff value, i.e., a number that is determined by registration system capabilities or was chosen as the lowest value from which the value sequence (1) is statistically representative. Consider that values (1) are distributed according to the Gutenberg–Richter type law(2)PrR<x=Fx|R0,ρ,β=e−βR0−e−βxe−βR0−e−βρ,R0≤x≤ρ.In this case,x is an unknown parameter with a maximum potential value of R. When the dependency (2) is represented on doubly logarithmic axes, the unknown parameter is sometimes referred to as the “slope” of the Gutenberg–Richter type law at small values of x. The second assumption is that sequence (1) is a the Poissonian process with an intensity value that is likewise unknown. It is important to note that the initial earthquake series is not a Poissonian process. As a result, before the estimations, a preparatory procedure comparable to aftershock removal is performed. This procedure is detailed below [32]. As a result, the whole unknown parameter vector is as follows:(3)θ=ρ,β,λ.Let us now propose an error for which we have values (1) that are specified by a formula(4)R¯=R+ε,and let nx|δ| be a density of probabilistic distribution of the error ε, where δ is a given scale parameter of the density. We now utilise a normal distribution function, which is more usual in the seismic hazard assessment, instead of a uniform probability distribution function (to speed up computations)(5)nx|δ=12δ,x≤δ,0,x>δ.The Bayesian approach for estimating the parameter is based on the Bayes formula [33](6)fpθ|R⟶n,δ=gcR⟶n|θ,δfaθg0R⟶n|δ,where faθ is an a priori probability distribution function (p.d.f) of parameter vector, which contains our understanding well about the statistical model, which is dependent on some of the most basic assumptions. Here, faθ is an a priori probability distribution function (p.d.f) of parameter vector θ which carries our information about the statistical model that is based on some of the most general assumptions. Taking into consideration the vector of observations R⟶n, we wish to make this knowledge more exact. The more precise information about the vector θ is contained within its a posteriori p.d.f. fpθ|R⟶n, δ. When the number of observations n⟶∞ increases, the variance of the a posteriori p.d.f. decreases, implying a more exact estimate. Function gcR⟶n|θ,δ is a p.d.f. of observations R⟶n , which is following from the used statistical model. Function gcR⟶n|θ,δ is a p.d.f. of observations R⟶n, that is regarded independent of any model. It might be determined using the complete probability formula(7)g0R⟶n|δ=∮ngcR⟶n|ϑ,δfaϑdϑ,where п is an a priori domain where vector of parameters θ can take its values.When p.d.f. fpθ|R⟶n,δ is computed, the estimate of the vector θ could be calculated as a mean value over a posteriori distribution(8)θ^R⟶n|δ=∫Πϑfpϑ|R⟶n,δdϑ.Furthermore, the Bayesian estimate’s variance might be determined by calculating(9)Varθ^R⟶n|δ=∫Πϑ−θ^R⟶n|δ2fpϑ|R⟶n,δdϑ.The expression for p.d.f.gcR⟶n|θ,δ for uniform distribution of errors has a rather long although not very complex derivation, and it is given in the papers [4, 23, 24]. A priori p.d.f. faθ is taken as uniform within a rectangular 3-dimensional domain(10)λmin≤λ≤λmax,βmin≤β≤βmax,ρmin≤ρ≤ρmax.The rules for obtaining boundary valuesλmin , λmax, βmax,βmin, and ρmin, ρmax for maximum peak ground acceleration problems are given in the papers [4, 24]. Besides the Bayesian estimate of vector (3), the following problem could be solved by the Bayesian approach as well. Let us denote by RT a maximal value of R on the time interval [0, T]. Then, ProbRT<X=exp−λ1−FX|θT. However, a situation where there are no occurrences on [0, T] is also included in this likelihood. Let ϑT be the number of events with RC R0 on the interval [0, T]. Then,(11)PrϑT=0=e−λT;PrϑT≥1=1−e−λT.That is why,(12)∅Tx|θ=PrRT<x|ϑT≥1=exp−λT1−Fx|θ−exp−λT1−exp−λT=exp−λTFx|θ−1expλT−1,where function Fx|θ is given by formula (2). Let us consider a priori quantile YTα|θ of probability α for maximal values of R on time interval [0, T].(13)φTYTα|θ|θ=α,0≤α≤1.Applying the Equations (7)-(9) to any given values of T may compute a priori the PDF, mean value and variance for quintile YTα|θ. All 3D integrals in equation (6) were determined by breaking a priori domain (10) into small cubes by a grid of 30 × 30 × 30 nodes and taking the values of integrated functions inside the centres of these small cubes using a simple integral sum.


## 3. Seismic Hazard Estimate for North and Northwest Regions of Iran

Estimating the maximum expected ground acceleration is fundamentally different from the maximum magnitude value. Reference [34]. Direct measurement of seismic acceleration is very limited and independent [5]. Therefore, no catalog is available that contains the maximum acceleration values for the desired site. But there are a number of reduction rules that represents the functions between the logarithms of the maximum ground acceleration with the magnitude of the earthquake (M), and there are different distances to the earthquake site and other factors at the site of the earthquake. These relationships are shown in Formula equation (13).(14)R=LogAmax=φM,D.Usually these functions collect data from a specific area and are obtained experimentally using regression laws (Joyner and Boore [35]; Fukushima and Tanaka [36]; Theodulidis and Papazachos [37]). Initially using the catalog [38] and also other seismic database centers (IGUT [39]; IIEES [40]; ISC [41], and BHRC [42]), an earthquake catalog was prepared and used for use in this study. Because the existing catalog covers a long period of time, it includes temporal and spatial heterogeneities. The step method is used to calculate the completeness of the catalog. In this method, it is assumed that the earthquakes in a period of time are stable and that the rate of events in a period and large in the period is constant. In this study, the complete start time table provided by [43] was used to apply the step method. For the Poisson process of events, which has occurred as the main hypothesis in this method, to be realized, time-dependent events must be removed from the sequence (1). To do this, the aftershocks are removed using the modified Gardner and Knopoff methods. In other words, in this method, all aftershocks are not removed. As we always have a sequence of major earthquakes and aftershocks, in this sequence, an earthquake will remain that will produce more acceleration in the intended construction. Sometimes, one of the aftershocks is able to produce more acceleration due to its proximity to the target site. In this particular case, the aftershock is not removed, Gardner and Knopoff write [44].After providing the earthquake catalog with the specifications mentioned above, for measuring the seismic hazard index, such as the maximum acceleration of the Earth’s surface, reduction curves are used. Global examples of these relationships include the New Generation Reduction (NGA) relationships that have been introduced. It should be added that the relations are based only on Iran’s records. They cannot provide accurate estimates of acceleration or other seismic indicators near the fault. Because there are a few records near the fault in Iran and if these records are also classified into different soils, there will not be enough to perform a nonlinear regression. Hence, now for the north and northwest regions of Iran, this has several active faults in and around this region. There is no choice but to use relationships New Generation Reduction (NGA). In general, the relationships between Campbell and Bozorgnia [45] and Boore and Atkinson [46] are more consistent with Iranian data. It must be considered that the Boore and Atkinson (2008) relationship is more coordinate with Iran earthquakes data especially in the Alborz-Azerbaijan tectonic seismic zone [47], Figure 2.Figure 2
The Alborz-Azerbaijan province [43] which is the studied state in this article is marked in yellow.Therefore, to reduce the seismic hazard in the north and northwest part of Iran, the reduction relationship of Boore and Atkinson [46] has been used. It is defined for magnitudes M > 5 and is shown in Figure 3. Now, with the catalog of earthquakes in the study area and the appropriate reduction relationship appropriately, the selection method of which is completely above, we form the sequence of equation (1). The values of sequence (1) are log (Amax). This sequence consists of the calculated values of equation (13). As mentioned, sequence (1) is the result of a reduction relationship.Figure 3
Graphics for the attenuation law.In reduction relationships, they always make a certain mistake and an inherent error results from statistical fits. In general, the errors resulting from the reduction relationship are divided into two parts, the first part contains an inherent error related to statistical fluctuations in the data, and the second part is related to the incompleteness and inadequacy of the function class selected for these relationships. In this study, these general errors are considered uniform and zero. And in contrast the value ofδ for the uncertainty of the value of Log (Amax) in this method, we put the reduction relationship almost twice the amount of the standard deviation.Thissdudy takes the parameter ρmax, the maximum possible value of Log (Amax) and with probability α = 90% for the time interval T = 475, 100, 50 years at 35° ≤ Lat ≤ 40°N; and 44° ≤ Lon ≤ 56°E estimated. For each station of the network, a sequence of Log (Amax) has been calculated using the catalog list of earthquakes in the north and northwest part of Iran in the area mentioned above, which includes historical earthquakes. The next step is to remove aftershocks in order to provide a random earthquake sequence, and 20 major earthquakes were considered for a maximum log value (Amax) for analysis. Therefore, for each station of the network in formula (1), we have the value n= 20. But the different values R0Ri=max1≤i≤nRi are the prediction limits for ρ. ρmax=Rτ+0.5. If the number of events is considered large, the experimental tail distribution function approaches the Gutenberg–Richter function.


## 4. Results and Discussion

At the regional level, there are numerous methods for determining seismic danger which usually investigate the possibility of these events. In this study, using the Bayesian approach, we have predicted the maximum acceleration of the Earth’s surface and the areas with the highest seismic risk. Inputs in this have been the method catalog of past earthquakes, which includes historical earthquakes and instrumental earthquakes in the region. Select the appropriate reduction relationship for the region, which here is the reduction relationship of Bohr and Atkinson [46]. This approach is an accurate tool for integrating fundamental seismic activity information with previous seismic event data. Given that the Bayesian theory provides a precise approach for integrating seismic history information, such information may be used to complete seismic data sets. Whether this information is subjective, geological, or statistical, it can be combined with historical observations of earthquakes. For this reason, the method used in this study is much better than other conventional methods and does not require intermediaries to obtain the maximum acceleration of the Earth and its values. Figures 4(a)–4(d) show a graph of experimental tail probability functions for 4 stations of this design network.Figure 4
Empirical tail probability functions for log (Amax) values within 4 nodes of the grid. (a) Tehran station. (b) Tabriz station. (c) Sari station. (d) Rasht station.
(a)(b)(c)(d)The best-fit exponential law curves are the dotted lines. As a result, the required maps may be created after generating estimations for each network station. The map of 90% quintile of distribution of maximum values of log (Amax) on the future time interval of the length T= 50, 100, and 475 years Figures 5–7.Figure 5
Map of 90% quintile of distribution of maximum values of log (Amax) on the future time interval of the length T= 50 years.Figure 6
Map of 90% quintile of distribution of maximum values of log (Amax) on the future time interval of the length T= 100 years.Figure 7
Map of 90% quintile of distribution of maximum values of log (Amax) on the future time interval of the length T= 475 years.Table1 shows the PGA value for some important cities in the north and northwest part of Iran in terms of gall. The city of Rasht and Zanjan, despite the existence of active faults, have a small amount of instrumental and historical data. Historical earthquakes are reported only in the cities of Tabriz and Tehran. And according to the introduced attenuation relation, they have the possibility of creating more acceleration than instrumental earthquakes in that area. Due to the long return period of earthquakes, these two cities are more affected by historical data in their constructions.Table 1
Estimated earthquake values in some cities of the Alborz-Azerbaijan province.
NameLongLatPGA (g)50 yr100 yr475 yrArdabil48.32838.2170.0870.1180.175Gorgan54.38536.8390.0850.1090.201Karaj50.9635.820.1150.1380.218Orumiyeh45.0537.530.0780.0940.117Qazvin50.1137.2920.1080.1210.203Rasht49.59137.2920.0750.0830.125Sari53.0636.600.0870.1170.211Semnan53.4435.590.0830.1130.196Tabriz46.2838.050.2060.2860.389Tehran51.4035.8270.1170.1470.238Zanjan48.4936.670.0490.0610.085


## 5. Conclusions

In this research, a 90% probability level of PGA value distribution function for the next time intervals of 50, 100, and 475 years with the effect of historical earthquakes has been obtained by the Bayesian statistical method and by coding Lyubushin as a probability. In general, when we have relatively appropriate data in an area, the Bayesian statistical method estimates a higher value than the conventional probabilistic method. But for cities such as Rasht and Zanjan, where we do not have enough data, this method estimates a lower value than the modified probabilistic method. It should be noted that the small number of events in cities such as Zanjan and Rasht is due to the ability of Bayesian statistics. At higher probability levels, due to the greater uncertainty of the PGA obtained, it approaches the values obtained from the classical probabilistic method, according to hazard maps with 90% probability levels for the next 100 years. The Bayesian method is a manual based on the catalog of earthquakes and the selection of a suitable attenuation (reduction) relationship at the desired site. This work is the estimation of the maximum ground acceleration at the probabilistic level of 90% (galls) in the bedrock. Its probable distribution is for a period of 50, 100, and 475 years.



---

*Source: 7046067-2022-07-12.xml*