Spatial Correlation of Peak Ground Motions and Pseudo-Spectral Acceleration Based on the Sarpol-e-Zahab Mw 7.3, 2017 Earthquake Data

Ground motion intensity measures and structural responses are correlated in nearby sites. The value of this correlation relies on some parameters such as the local geology and distance between the two sites and the natural period of structures, particularly, when lifeline systems or distributed structures are of concern, the issue becomes more important. In this study, the spatial correlation of peak ground acceleration and Spectral Acceleration are evaluated as a function of inter-site separation distance for the Mw7.3 Sarpol-e-Zahab earthquake. On 12 November 2017 a large earthquake occurred near the western border of Iran. The epicenter of the earthquake was reported at 34.77 N and 45.76 E. Here, 192 pairs of horizontal components from the above-mentioned event and 35 of its larger aftershocks with magnitude ranging from Mw 4.0 to 7.3 are employed to evaluate the intra-event residual correlation by considering two Ground Motion Prediction Equations (GMPEs) proposed by Akkar and Bommer [2010] and Zafarani et al. [2017]. A correlation analysis is carried out through semivariogram as a powerful geostatistical tool. As a skeleton for correlation modeling, a kind of exponential model is used. According to the proposed model, the results show that the overall trend of correlation Range depends on spectral period. The results demonstrate that there is strong spatial correlation in the proposed model obtained from the Sarpol-eZahab ground motions. The model provided in this study could be employed in earthquake engineering implements such as Shakemaps [Wald et al., 1999] whenever spatial correlation models are required.


Introduction
Ground-motion prediction equation (GMPE) is a remarkable parameter in probabilistic and deterministic seismic hazard and loss evaluations [Crowley et al., 2008a]. The functional form of a ground motion model is typically a function of distance, magnitude, type of soils, style of faulting and other existing parameters [Soghrat and Ziyaeifar, 2017].
GMPEs can be proposed according to the local, regional, or global databases. These models can estimate intensity measures such as peak ground acceleration (PGA), peak ground velocity (PGV) and spectral acceleration (SA) for a particular earthquake scenario. The uncertainties in the GMPEs are usually shown by the inter-event (between-events or earthquake to earthquake) and intra-event (within-event or record-to-record) variability [Sokolov and Wenzel, 2011;. Noted that, in another context uncertainties in GMPEs can be divided into epistemic and aleatory uncertainties. Epistemic uncertainty is due to the limited amount of observed data while aleatory uncertainty is because of the differences between observed and predicted data [Cimellaro et al., 2011]. Inter-event residual is the average shift of the recorded value in an earthquake from the median predicted value. Spatial correlation is usually defined for intra-event residuals. In other words, when an earthquake occurred, it is tried to find correlation in different recording stations. Therefore, the inter-event residuals excluded from analyses in this study. This is due to the fact that we lacked a sufficient number of data points to develop separate intraevent spatial correlation estimates for individual earthquakes [see Goda, 2011 for a fruitful discussion on inter-event variability of spatial correlation].
Moreover, the difference between the recorded ground motion at the specified station and the earthquake-specific median predicted is represented by intra-event residual [Chen and Faccioli, 2013]. The intra-event variability indicates that the intensity measures (IMs) at special sites are scattered around the event median. It is shown that the intraevent residuals are spatially correlated by comparing observed ground motion with a predicted relation. Also, it is resulted that with increasing distance between two stations the spatial correlation will be reduced Goda and Hong, 2008a;Goda and Atkinson, 2010;Cimellaro et al., 2011]. Wagener et al. [2016] developed a model according to the records belongs to small to moderate seismic event for intra-event spatial correlation. They found that the correlation coefficients of short-period SA decrease quickly with increasing inter-station distance.
Estimation of spatial correlation for IMs is essential for assessment of seismic risk of distributed building portfolios, infrastructure systems and lifeline networks (such as gas or electrical networks), water supply and transportation systems. Such studies result in finding the correlation between intensity measures from different earthquakes at different sites. Several researchers investigated the consequences of spatial correlation of ground motion IMs on loss estimation for the networks mentioned above. Crowley et al. [2005] showed that intra-event correlation causes greater variability in the estimation of total earthquake loss because of a particular ground motion scenario. Some studies have represented the significant influence of intra-event spatial correlation in the probability distribution of total seismic losses [Goda and Hong, 2008b;Park et al., 2007]. In other words, disregarding and underestimating the spatial correlation can miscalculate the losses (Bazzurro and Luco 2007; Lee and Kiremidjian 2007).
The loss estimation methods are based on ground motion intensity measures in which the measures can be obtained from GMPEs. Noted that, the GMPEs can also consider the effects of spatial variations if a rich database is available and a modified regression method is essential to reach this aim [e.g. the equations proposed by Kotha et al., 2016 andLandwehr et al., 2016]. Some researchers applied only a single earthquake scenario [Crowley et al. 2008a,b;Goda and Atkinson 2009] by analyzing the intra-event correlation effects, while others studied multiple classes, the spatial correlation of IMs will be stronger. Heresi and Miranda [2018] studied evaluating the event-toevent uncertainty of intra-event spatial correlations and showed high event-to-event variability in the correlation model factors. They showed that magnitude of an earthquake is a significant variable of the spatial correlation model parameter unlike the site classes, tectonic region, and style of faulting. Chen and Baker [2018] described their results according to analogous spatial variation evaluation using physics-based simulations from the Cyber Shake platform. They implied that the geological condition and directivity effect are essential for modeling the spatial correlations. Worden et al. [2018] evaluated the conditional multivariate normal (MVN) distribution, and then they employed the MVN to the particular problem of ground motion interpolation. They proposed an adaptable framework to estimate different intensity measures at random locations.
Taking into account the importance of local spatial correlation models, here we investigate the correlation characteristics of the Mw 7.3 Sarpol-e-Zahab earthquake which occurred on 12 November 2017 with dozens of aftershocks. In this work, spatial correlation relations for spectral accelerations of this earthquake and its larger aftershocks are developed. For this aim, the recorded data obtained from the main event and its aftershocks earthquake are compiled and processed to find a set of corrected ground motions. In the next step, the residuals for each IM are computed using the local and regional GMPEs proposed by  and Akkar and Bommer [2010]. In this work, semivariogram as geostatistical tools is used to compute spatial correlation for eight periods (ranging from zero to one second). Finally, spatial correlation models as a function of structural periods are proposed and compared with some models in the literature.

Spatial correlation Model
The computation of spatial correlation will be discussed briefly in this section. For more details, it can be referred to Jayaram and Baker [2009], Goda and Atkinson [2010] and Markhvida et al. [2018]. To model the spatial correlation, we need to describe the GMPE and Intra-event residuals which they are defined in the following. Generally, GMPEs can estimate the intensity measures in site i because of earthquake j as shown in the following form: (1) Where the predicted median of ground motions is shown by . It depends on Magnitude (M), source-tosite distance (R) and other available parameters (θ) such as type of soils and style of faulting. is the observed ground motion parameter of interest such as PGA, PGV or SA. Also, and are intra-event and inter-event residuals with zero mean and standard deviations of and , respectively. Noted that the overall standard deviation ( ) is computed by equation (2): (2) Generally, GMPEs have been proposed according to different database. So, three categories can be considered for these equations including models proposed specifically for the Iranian plateau, the ones proposed for the Middle East and Europe region and global ground motion relations proposed by the "Next Generation of Ground-Motion Attenuation Models" (NGA). In this study, two models proposed by  and Akkar and Bommer The semivariogram, which applied to find the dissimilarity or decorrelation between data [Cressie, 1993], can be computed assuming the hypothesis of second-order stationary and isotropy using equation (3).
where a and (ℎ) is the sill and the correlation coefficient between intra-event residuals at two different sites separated by h distance, respectively. An assessment of the experimental semivariogram for given IM (e.g. SA) based on dual empirical observations at two adjacent sites x and x+h is obtained as follows: (4) Here (ℎ) is the estimated stationary semivariogram, ( ) -( + ℎ) represents the difference between intraevent residuals at sites with distance h; (ℎ) is the set of pairs of adjacent sites separated by the same distance ℎ, and | (ℎ)| is the cardinal of (ℎ). The formula is resulted from the method-of-moments and is referred to as the classical estimator [Matheron, 1962]. Because the summation is over a squared term, the estimator will be seriously affected by typical observations (LJ Young and J Young, 2013). Therefore, Cressie and Hawkins [1980] developed a more robust estimator of semivariograms, which is less sensitive to outliers, as follows: (5) Note that in the above formulation a common semivariogram in variant through earthquakes is assumed for different events. However, individual events are treated separately in computing the empirical semivariogram, that is, ( ) -( + ℎ) differences are taken considering only intra-event residuals from the same earthquake without mixing up the events [see Figure 4 in Esposito and Iervolino, 2011].
The recorded logarithmic intensity measures follow a lognormal distribution that may be shown as ln = (ln ; ). As a result, the normalized intra-event residuals () can be calculated as follows: (6) Noted that, the normalized intra-event residuals are approximated by equation (6). We know that inter-event residuals are constant during a seismic event for each site, this approximation is reasonable and does not change the results.
The spatial correlation of the normalized intra-event residuals can be studied using semivariogram.
Semivariogram is a geostatistical tool used to quantify the degree of divergence between observations as a function of distance. In other words, the semivariogram is used to find the dissimilarity or decorrelation between data.
Because normalized intra-event residuals have unit variance, we keep away from the evaluation of the sill, as it should be equal one [Jayaram and Baker, 2009]. Consequently, equation (3) becomes equation (7), where the superscript shows an empirical computation: After calculation of semivariograms in many discrete points, for practical usage it is important to fit a continuous function to the resulted measures, because of calculating semivariogram values for any probable separation h.
Several correlation models have been proposed in the literature which they can be considered Gaussian model, exponential model and spherical model as shown in equation 8.
where b is Range of the semivariogram function. The exponential model, as a most common one in literature is used for further analysis in this study.

Strong motion data
The Iranian plateau is one of the tectonically active regions which located along the Alpine-Himalayan orogenic belt where many large seismic events occurred in this region .
The Iranian strong motion network (ISMN) of the Building and Housing Research Center (BHRC) records the strong ground motion data. This network was equipped with Kinemetrics SMA-1 analog accelerographs, which have been gradually replaced and densified with SSA-2 and CMG-5TD digital instruments. Since 1975, more than 1100 stations are available across the country which recorded more than 10,000 accelerograms until now.
On 12 November 2017, an earthquake with moment magnitude of 7.3 occurred in the Zagros seismotectonic zone, near the western border of Iran (see Figure 1). This event happened near Sarpol-e-Zahab city. The epicenter was reported to be at 34.77 N and 45.76 E by the Iranian Seismological Center (IRSC) at depth of 18 km. Noted that this event is the largest earthquake over the last 100 years in the Zagros region [Farzanegan et al., 2017]. This event caused numerous landslides and building damages which leads to more than 600 deaths and 7000 homeless. People sensed this earthquake at distances about several hundred kilometres such as Tabriz and Arak. The maximum horizontal PGA of 0.69g was recorded at the Sarpol-e-Zahab station. In this study, we considered the main earthquake and its larger aftershocks (35 aftershocks with Mw more than 4.0 with at least two high-quality records) which their characteristics are shown in Table 1.
In the further analysis, the site categorization is considered according to the characteristics of the top 30 m

Evaluation of Spatial correlations and results
The predicted measures for the horizontal component of ground motions were calculated at different periods using the model developed by  and Akkar and Bommer [2010]. The geometric mean of horizontal components in the GMPEs is considered in the both proposed models. According to previous sections, computing the intra-event residuals is the first step in the analyses. Figure 3 shows the intra-event residuals for the relations proposed by  and Akkar and Bommer [2010] at four different periods including T = 0 (or PGA), 0.3, 0.6 and 1s.
After calculation of the residuals, the Ranges of semivariograms computed using the earthquakes introduced in the strong motion data section. Journel and Huijbregts [1978] have been suggested half of the maximum distance between sites in the dataset and they recommended to consider the number of bins so that there are at least 30 pairs per bin. Different distance bins are tried in this study for the analyses. Different bin widths from 1 to 7 km were tried, and finally the bin of 7 km and the maximum distance of source-to-site of 120 km was considered in the present study. Figure 4 demonstrates the number of pairs in each bin according to separation distances. Noted that, for fitting experimental values, the exponential model was chosen for the analyses that were used in similar studies [e.g., Jayaram and Baker, 2009;Iervolino, 2011, 2012].
7 Spatial Correlation of Spectral Acceleration By assuming the nugget effect in the modeling, the Range b (equation (8-b)) is required to determine. To estimate the b some methods such as least square fit, weighted least square fit, and the manual fitting method has been proposed in the literature [see Jayaram and Baker, 2009;Du and Wang, 2013;Wang and Du, 2013]. In this study, we used the least square method to fit the exponential model (equation (8-b)) on provided data obtained from the robust estimator. Table 2 shows the model parameters for the exponential model (equation (8-b)) obtained from using the robust estimator and the classic estimator results. As observed in previous studies, such as Esposito and Iervolino [2012], the results of two estimators are not so different and therefore the only the robust estimator results are considered.  . Thisobservation (no clear trend between Ranges and periods) has been consistent with past studies in the topics of spatial correlation of spectral acceleration [Esposito and Iervolino, 2012;Zerva and Zervas, 2002]. As Jayaram and Baker [2009] expressed in their study, there are significant differences in the obtained Ranges depending on the ground-motion time histories Hamid Zafarani et al.
8 Figure 4. The number of data pairs as a function of site-to-site separation distance for different bin width a) 1 km, b) 2 km, c) 3 km, d) 4 km, e) 5 km, f) 6 km, g) 7 km. used at short periods (below 1 s). There is no certain trend in the study for whole periods less than 1 sec which the most structures in the studied region are in this range of period. According to the results, three segments can be considered to investigate the trend of correlation Ranges including 0 ≤ T ≤ 0.1s, 0.1s < T ≤ 0.4s, and 0.4s < T ≤ 1.0s. The Ranges in the first part is increasing up to T = 0.1 s which is in harmony with Du and Wang [2013], and Esposito and Iervolino [2012]. The Ranges in the second segment has a decreasing trend with the increasing period which has been observed in recent studies such as Du and Wang [2013], Jayaram andBaker [2008, 2009], and Zerva and Zervas [2002]. The trend of Ranges in third part is incremental similar to the model developed by Du and Wang [2013], and Jayaram and Baker [2009].
Experimental semivariogram values using classical estimator and exponential fitted model at different period by considering the model proposed by Akkar and Bommer [2010] shown in Figure 6.
Moreover, experimental semivariogram values using classical estimator and exponential fitted model at different period by considering the model proposed by  shown in Figure 7.
It seems that the values estimated by classic and robust estimators are similar except some cases (as shown in Figures 6 and 7). The classic estimator can be affected by undesirable effects while the robust estimator proposed by Cressie and Hawkins [1980] is less sensitive to outliers. Furthermore, an exponential model is fitted to the results of Robust estimator according to the study of Esposito and Iervolino [2012].
9 Spatial Correlation of Spectral Acceleration Finally, in order to show the regional dependency of the correlation model, a comparison with some correlation models in California, Italy, Europe, Taiwan and Turkey is done in Figure 8. Figure 8 shows a comparison between correlation model fitted to the Sarpol-e-Zahab earthquake data with using two different GMPEs and those models from pan-European data [Esposito and Iervolino, 2012], Italy [Esposito and Iervolino, 2012], California and Taiwan [Jayaram and Baker, 2009], California [Goda and Hong, 2008] and Turkey [Wagener et al., 2016].
10  The ESD data in this Figure, correspond to records used in the study of Esposito and Iervolino [2012] and the Akkar and Bommer [2010]. ITACA dataset considered by Esposito and Iervolino [2012] contains 763 ground motions from 97 events with Mw 4.0-6.9. Goda and Hong [2008] used 592 ground motions from 39 events in Western North America. Jayaram and Baker [2009] employed the ground motion records of seven earthquakes including Anza earthquake, Alum Rock earthquake, Parkfield earthquake, Chi-Chi earthquake, Northridge earthquake, Big Bear City earthquake and Chino Hills earthquake. The data used in Wagener et al. [2016] included 372 records of eight earthquakes that were obtained from the Istanbul earthquake which occurred between 2003 and 2013 with Mw 3.5-4.8. Generally, California, Italy, Europe, Taiwan and Turkey data shows less correlation than Sarpol-e-Zahab earthquake data which attenuate more rapidly. This may be due to the different characteristics of seismic sources, e.g. different stress drops, and also different crustal properties.   Figure 9a, the result of proposed models in this study [using GMPEs proposed by Akkar andBommer, 2010, andZafarani et al., 2019] are consistent with Du and Wang [2103] including ground motions from nine other earthquakes occurred in Taiwan, California and Japan. Comparison between the results of our models with Wang and Takeda [2005] in Figure 9b shows that the model proposed by  is in the upper limit while the model proposed by Akkar and Bommer (2010) is in lower limit. From the epistemic uncertainty point of view, we recommend to use both correlation models in practical applications.
The results of this study show that there is strong spatial correlation in the proposed models obtained from the Sarpol-e-Zahab earthquake ground motions. Therefore, these models should be used for hazard analyses and risk assessment studies in this region. Noted that, ignoring the spatial correlations in hazard evaluations for distributed structures (such as lifelines) can lead to significant distortion of results.
12 Figure 8. Comparison between correlation function (ρ(h))fitted to the Sarpol-e-Zahab earthquake data with using two different GMPEs and those models from Esposito and Iervolino [2012], Jayaram and Baker [2009], Goda and Hong

Conclusions
The spatial correlation models are proposed for western Iran region based on the Sarpol-e-Zahab earthquake and its aftershocks data (36 seismic events). Evaluation of the spatial correlation for the horizontal component of spectral acceleration (SA) at different spectral periods has been performed using geostatistical tools. Based on the results, an exponential function was selected for the provided models. A similar trend of correlation Ranges for the horizontal component of ground motion has been observed in comparison with other studies. The results show that the obtained Ranges can be considered in three distinct parts. The first and third segments have an increasing trend although the trend in the second one decreases with increasing spectral period. The Ranges increase up to 0.1s and more than 0.5s while the Ranges decrease in the period of 0.1 to 0.5 s. As it is clear, there is the strong regional dependency in the models, which implies that the local correlation models should use for investigation of damage patterns caused by the Sarpol-e-Zahab earthquake. According to the model developed in this study, there is strong spatial correlation along with the Sarpol-e-Zahab earthquake ground motions. For that reason, this model can be employed in local risk assessments studies for lifeline systems or distributed structures.