Ambient noise surface wave tomography of the Makran subduction zone , south-east Iran : Implications for crustal and uppermost mantle structures

Seismic ambient noise of surface wave tomography was applied to estimate Rayleigh wave empirical Green’s functions (EGFs) and then to study crust and uppermost mantle structure beneath the Makran region in south-east Iran. 12 months of continuous data from January 2009 through January 2010, recorded at broadband seismic stations, were analyzed. Group velocities of the fundamental mode Rayleigh wave dispersion curves were obtained from the empirical Green’s functions. Multiplefilter analysis was used to plot group velocity variations at periods from 10 to 50 s. Using group velocity dispersion curves, 1-D vS velocity models were calculated between several station pairs. The final results demonstrate significant agreement to known geological and tectonic features. Our tomography maps display low-velocity anomaly with SW-NE trend, comparable with volcanic arc settings of the Makran region which may be attributable to the geometry of Arabian Plate subducting beneath the overriding the Lut block. The northward subducting Arabian Plate is determined by high-velocity anomaly along the Straits of Hormuz. At short periods (\20 s), there is a sharp transition boundary between lowand high-velocity transition zone with the NW trending at the western edge of Makran which is attributable to the Minab fault system.


Introduction
The Iranian plateau is subject to several tectonic episodes, including active stages of intense folding (e.g., in the Zagros region), faulting and different types of tectonic domains. The Makran subduction zones is located in the SE of Iran, extending from the Main Zagros Thrust (MZT) to the western end of the Makran wedge and to the Ornach-Nal and Chaman fault zones in southwestern Pakistan (see Fig. 1). The transition between the Zagros continent-continent collision zone and the western Makran subduction zone is marked as the Minab fault (see Fig. 1b). This complex was formed by subduction of the oceanic crust of the Arabian Plate under the Eurasian Plate (Fig. 1a). The Arabian Plate convergence rate, which occurs in a north direction towards Eurasia, is *23-25 mm/a according to GPS measurements McClusky et al. 2003;Vernant et al. 2003;Masson et al. 2007). The Makran subduction zone subducts at a higher convergence rate compared with the continental collision portion of the Arabian Plate boundary in the Zagros Suture Zone. In the Makran region, the convergence rate increases from west to east . Geodetic data provide selected evidence that nearly all of the convergence of the Arabian Plate is accommodated in the eastern boundary of Makran, whereas only half of the convergence is accommodated in the west. At the eastern boundary of Makran, the convergence rate is *42.0 mm/a (Fig. 1a), as estimated by DeMets et al. (1990).
Notably little seismic tomography has been performed in the Makran subduction zone, especially on the structure of the upper mantle. Compared with the remainder of Iran, relatively little is known about the structure of the Makran subduction zone. Most of the tomographic studies performed in the Makran region are limited to global tomography surveys with low resolution. Few shallow seismic investigations have examined the sedimentary structure of the Makran belt to reveal the heterogeneity of the crust and upper mantle, and few seismic tomography studies have been conducted on the crustal and upper mantle structure of the Makran subduction zone. Giese et al. (1984) studied the Moho depth using a refraction profile consisting of sparse recordings along a line from central Iran to the Straits of Hormuz and indicated a crustal thickness of 40 km beneath central Iran. Using gravity measurements in the seismic results of Giese et al. (1984), Dehghani and Makris (1984) prepared a Moho map of the Iranian plateau and found that the crust beneath the Lut depression is less than 40 km thick. Snyder and Barazangi (1986) used the same data and found that the Moho depth was nearly 40 km beneath the Persian Gulf (Maggi and Priestley 2005). The crustal thickness of the Makran region is less well known. Few studies of the deep structure of the upper mantle exist for this area. Recent surface waveform tomography (Shad Manaman et al. 2011) indicated that the crustal thickness beneath the Oman seafloor and Makran fore-arc setting is *25-30 km and is increasing to the volcanic arc. The Moho depth increases up to *48-50 km under the Taftan-Bazman volcanic arc (see Fig. 1b), where the subducting plate bends. Similarly, from the fore-arc setting to the volcanic arc in eastern Makran, the Moho depth increases to *40 km. Regional tomographic studies of the Iranian plateau do not provide detailed information on the structures in the crust and upper mantle due to the lack of welldocumented earthquakes. The lateral resolution is usually limited in these studies, e.g., on the order of 200-km horizontal grid spacing (Maggi and Priestley 2005) and 4 (Shad Manaman et al. 2011) in the Iranian plateau. In the Makran region, the resolution is even more limited because the path coverage is not dense due to the inconvenient distribution of stations and limited seismicity. Compared with our understanding of active tectonics, much less is known about the Makran subduction zone, and the seismic behavior of the Makran subduction zone has remained largely unknown. Additionally, basic limitations exist in the earthquake-based surface wave tomography because of limited seismicity and poor station coverage.
Recent studies use ambient noise to extract the surface wave empirical Green's functions (EGFs). Rayleigh and Love waves obtained from these EGFs can provide important information on the 3D shear wave velocity structure in the crust and uppermost mantle on both a . Black and yellow triangles show Tajikistan National Seismic Network (TJ) and the Virtual European Broadband Seismic Network (VEBSN). Ridge, trench, and transform boundaries are indicated by red, green, and blue lines, respectively. Plate motions are calculated in www.unavco.org based on APKIM2005 plate motion model (Drewes, 2009). b Topography map of the study area. WM Western Makran, EM Eastern Makran, SSZ Sanandaj-Sirjan Zone, UDMA Urumieh-Dokhtar Magmatic Arc, MZT Main Zagros Thrust, JAZ M. Jaz Murian., SH Straits of Hormuz Earthq Sci global and regional scales. Shapiro and Campillo (2004) inferred Rayleigh wave dispersion curves from ambient noise data recorded in station pairs separated by distances from *100 to [2000[ km. Shapiro et al. (2005 and Sabra et al. (2005) extracted time-domain Green's functions from ambient seismic noise and estimated the group velocity of Rayleigh waves at a local scale. Certain studies applied ambient noise tomography at short-and intermediate-band (e.g., 5-40 s) periods Yao 2012). However, other studies performed ambient noise tomography over longer periods, e.g., up to 60 s in China (Sun et al. 2010) and up to 70 s in the US (Bensen et al. 2008). Lin et al. (2008) presented seismic ambient noise tomography (ANT) for Love waves. Nishida et al. (2008) also analyzed Love waves to study seismic ambient noise tomography, including the crustal overtones.
Recent theoretical studies demonstrate that under the assumption, that the seismic noise wavefield is diffused, the empirical Green's functions between two stations can be estimated by correlating the noise recordings from these two sites (Weaver and Lobkis 2001;Derode et al. 2003;Snieder 2004). If sufficient stations exist, the study of ambient noise seismic waves yields results with a resolution higher than that of traditional surface wave tomography methods (e.g., Yang et al. 2007). The earthquake-based surface wave tomography has certain limitations, such as difficulty in obtaining high-quality short-period (\20 s) dispersion measurements from teleseismic events and inaccuracy in earthquake hypocentral locations and moment tensors, especially for small events. These limitations affect the inversions of seismic surface waves. Additionally, seismic surface waves only sample certain azimuths due to the uneven distribution of earthquakes. The purpose of this study is to produce Rayleigh wave group velocity maps with higher resolution than previous surface wave maps produced in the Makran region to study the crust and the upper mantle structures in the Makran subduction zone.
In this study, the ANT is applied using data from periods of *10-50 s collected from the recordings of 41 stations between January 1, 2009 and January 1, 2010. Dispersion curves were measured using the fundamental mode of Rayleigh waves extracted from the ambient noise and inverting them to obtain 2D group velocity maps for the crustal and upper mantle structures of the region. We also investigated the directionality and seasonal variations of the noise sources. The difference between the causal and acausal components of the cross-correlation results of station pairs was examined to measure the main direction of the energy flux across the region. Using the group velocity dispersion curves, 1-D v S velocity models were calculated between several station pairs. Finally, the resulting group velocity maps were used to infer crust and upper mantle structures in the Makran region.

Data processing
This study is based on a variety of seismic sensor data collected from digital broadband instrument (BH) recordings from the International Institute of Earthquake Engineering and Seismology (IIEES), which is equipped with a CMG-3T broadband sensor (0.01-100 s), as well as seismic data from the Global Seismic Network (GSN), the Virtual European Broadband Seismic Network (VEBSN), the GEOFON program of GFZ Potsdam (GEOFON), and the Tajikistan National Seismic Network (TJ), as depicted in Fig. 1a. The EGFs were calculated for one year of continuous vertical component seismograms recorded from January 1, 2009 to January 1, 2010. Use of vertical component seismic data implies that the resulting cross-correlations contain only Rayleigh wave signals. We followed the data processing procedure described in detail by Bensen et al. (2007). First, the continuous noise data were cut into one-day data files, and those data with gaps of less than 10 s were chosen. The instrumental responses were subsequently removed from all data, followed by decimation of the data to one sample per second. Decimation reduced the amount of storage space and computational time required. The next step involved removing the trend and mean value, zero-phase Butter-worth high-pass filtering with a corner frequency of 0.01 Hz, and whitening and bandpass filtering around the target frequency (period from 10-50 s) as a function of the inter-station distance (Cho et al. 2007;Pedersen and Krüger 2007). The next processing step included temporal or time-domain normalization to remove further contaminating effects of earthquakes on the noise correlations (Bensen et al. 2007). The results were stacked over the total time period available for each pair to produce the resulting time series.
To quantitatively evaluate the quality of the stability of the stacking process, we calculated the signal-to-noise ratio (SNR) for each cross-correlation. The SNR is defined as the ratio of the peak amplitude within a time window surrounding the expected arrival time of the fundamental mode Rayleigh waves at a given period to the root-mean-square of noise trailing the signal arrival window (Bensen et al. 2007). The signal window is determined using the arrival times of Rayleigh waves at the minimum and maximum periods of the chosen pass-band frequency. All empirical Green's functions in the 10-to 50-s period bands are plotted in Fig. 2 to evaluate the quality of the cross-correlation functions. Based on cross-correlations of broadband seismic records obtained at stations within or adjacent to the Pacific Basin, Lin et al. (2006) indicated that broadband ambient noise propagates coherently between island stations and between island and continental stations. The Green's functions observed for fundamental mode Rayleigh waves with high SNR establishes the physical basis for ANT across the Pacific, and any non-physical arrivals caused by long ocean-crossing paths were successfully rejected. Both positive and negative correlation lags show clear surface waves with an average apparent velocity of *3.0 km/s. The group velocities used to predict arrival times were calculated from the AK135 velocity model, (Kennett 1995), and the RMS noise level was measured in a 500-s noise window at the end of the signal.

Directionality
To obtain a better understanding of the distribution of the noise source in space and time, the seasonal variability of the relatively continuous noise was studied, a process necessary for optimization of the noise-based seismic tomography. Seasonal variations of the cross-correlations for periods of 10-20 and 25-50 s were computed for two station pairs located perpendicular to each other (BNDS-KRBR with the NS direction and KRBR-ZHSF with the EW direction) for different seasons (Fig. 3). A positive time delay indicates waves propagating from the coastlines to the continent (BNDS to KRBR). For the other pair of stations, the positive lags indicate signals propagating from KRBR to ZHSF. Considering a period band of 10-20 s, the cross-correlations exhibit a clear seasonal variation (Figs. 3). For the NS station pair, the apparent asymmetry of the data indicates that the energy originating from the coastlines is much larger than that from the continent. This result shows that an important contribution of the noise observed in the Makran region originates from the south and the coastlines, likely from the Persian Gulf and Oman Sea. For another pair of stations (KRBR-ZHSF), the directionality is significant but with different characteristics. In spring and summer, a preferred directionality does not appear; however, a clear pattern in autumn indicates that the energy flux flows from the east to the west, and this trend is completely the opposite in winter.
At longer period bands (25-50 s), the behavior of the noise varies in different seasons. The cross-correlation in this band is symmetric in spring and summer for BNDS-KRBR ( Fig. 3b), whereas in the period band of 10-20 s, the amplitude of the positive lag of the cross-correlations is much larger (Fig. 3a). Again, for KRBR-ZHSF, contrary to the observations in the period band of 10-20 s, the resulting using cross-correlation is not symmetric in spring (Fig. 3c,d). Moreover, in autumn, Fig. 3d indicates a wave propagation direction opposite to the results calculated in the period band of 10-20 s (Fig. 3c). This result suggests that the primary microseism is not recorded by the same process that generates the longer period noise (Stehly et al. 2006). The primary microseism has periods similar to those of the main swell (10-20 s). Therefore, it is believed that the primary microseism is related to the interaction of the sea waves with the coast (Gutenberg 1951). The long-period noise or ''the hum'' has been attributed to the so-called infragravity waves, i.e., the ocean wave mode that exists at long periods (Webb et al. 1991). Thus, the results of our analysis show that the sources of the primary microseism exhibit seasonal variability that differs from that of the long period noise (hum).
In practice, for different frequency bands, the observed distribution of ambient noise can be far from homogeneous (Stehly et al. 2006). Therefore, it is necessary to determine whether the ambient noise is sufficiently distributed in azimuth to return unbiased dispersion measurements for use in tomography. To quantify the effect of the strongly anisotropic background noise source distribution, Yang et al. (2008) performed synthetic experiments and found that in the presence of low level and homogeneously distributed ambient noise, \0.5 % of the measured phase velocities are affected by much stronger ambient noise in an off-axis direction. Therefore, we must show that in all period ranges studied, the useful amount of ambient noise signals in all azimuths is greater than 50 %.
To investigate the directions of the incoming ambient noise, we plotted the azimuthal distribution of the SNR for the positive and negative components of each cross-correlation for the four period bands of 10-20 s, 20-30 s, 30-40 s, and 40-50 s in the northern winter (Oct. to Mar.) and northern summer (May to Sept.) of 2009 (Fig. 4). The length of each line is the amplitude of signal and the angle points in the direction from which the energy arrives. Each 20 azimuth bin shows the number of paths for both the inter-station azimuth (causal) and back-azimuth (acausal) components of the cross-correlation functions. Following the work of Bensen et al. (2008), the average Rayleigh wave EGFs with SNR[10 was computed at all four periods, and to compute the average fraction of yearly EGFs, the number of paths with SNR[10 in a given 20 azimuth bin was divided by the total number of paths in that bin. The average results over all azimuths at the four period bands of 10-20 s, 20-30 s, 30-40 s, and 40-50 s were 0.53, 0.64, 0.69, and 0.51, respectively. In other words, these values reveal that the fractions of relatively high SNR paths in all azimuths are greater 50 % in all period ranges studied, and hence, the useful amount of ambient noise signals is sufficiently distributed in different azimuths.

Group velocity measurements
In the next step, multiple-filter analysis (Herrmann and Ammon 2013) was used to measure group velocity dispersion curves. Each of the frequency components of the surface wave is sensitive to a different depth interval. In general, wave components with longer wavelength that propagate deeper will travel faster than the shallower ones because the seismic velocity of the Earth increases radially downwards. A technique known as phase-matched filtering was applied to determine the correct dispersion curve. The waveforms were narrow-bandpass-filtered with the Gaus- where f c is the center frequency. A trade-off exists between resolution in the time and frequency domains that is caused by filtering, i.e., larger values of a enhance the resolution in the frequency domain, whereas it decreases resolution in the time-domain (Herrmann 1973;Levshin et al. 1989). For distances between 200 and 3000 km and 3000-5000 km, a = 25 and 50 s were found suitable for our measurements, respectively. Selected dispersion curves are plotted in Fig. 5, and the corresponding paths are depicted in Fig. 1a, b.
The estimated uncertainties for the group velocities are based on seasonal variability because dispersion measurements from cross-correlations of ambient noise are naturally repetitive. To analyze the uncertainty, we selected 12 overlapping three-month time series for each station pair. The three-month time windows are reliable for obtaining dispersion measurements and also contain the seasonal variation. Fig 6 shows group velocity measurements for four station pairs (BNDS-ZHSF, GHIR-UOSS, MSEY-SHRT, and NASN-KRBR) obtained in twelve three-month cross-correlations that were bandpass filtered from 10 to 50 s periods. The paths for the station pairs are depicted in Fig. 1a, b. The one-year measurement is plotted as a black line with error bars, indicating the computed standard deviation. The standard deviation is computed on all sequential three-month stacks, and following the work of Yang et al. (2007), it was computed for a station pair if four SNR values of the three-month stacks exceed the criterion (in our study, SNR[10). The measurement was rejected if the standard deviation was either so large or so small that it could not be obtained in the three-month stacks. The uncertainty tends to increase with the period, possibly because of the decreasing amplitude of ambient noise at periods greater than 20 s ). If the standard deviation was more than three times the average of the standard deviations taken over all measurements, it was rejected because this indicates instability in the measurement (Bensen et al. 2008).
The 1-D crustal and upper mantle v S structures were subsequently calculated between different station pairs using the surf96 package (Herrmann and Ammon 2013). The initial v S model was parameterized into 4-km-thick   Fig. 1a, b Earthq Sci layers from surface to 240 km and then a half space below 240 km. Group velocity dispersion curves were input data for the inversion procedure. The reference v S model refers to the AK135 model (Kennett 1995). For the surf96 package, during the v S structure inversion, it was also necessary to input other parameters. The input density was assigned using the parameter of the background model (AK135), which gradually increases from 2.3 g/cm 3 in the shallow crust to 3.3 g/cm 3 in the deeper crust and is set to 3.3 ± 0.1 g/cm 3 in the upper mantle. In addition to the shear wave speed (v S ), the P-wave speed (v P ) and density (q) are also required for inversion. These parameters were assigned based on parameters of the background model (AK135). An example of the 1-D v S velocity models calculated between two station pairs (GHIR-ZHSF and KRBR-ZHSF) from the current dataset is presented in Fig. 7a next to the observed and synthetic dispersion curves (Fig. 7c). The paths related to two station pairs are depicted in Fig. 1b. The sensitivity kernels for different periods were also calculated from the AK135 velocity model (Fig. 7b). Because Rayleigh waves are more sensitive to shear wave velocities than compressional velocities, except in the uppermost crust (Bensen et al. 2009;Lin et al. 2007), we fixed the v P /v S ratio of each layer in accordance with the Ak135 velocity model, and the P-wave velocity was subsequently updated by the v P /v S ratio of the initial model in successive iterations.
The parametric test was used to determine the influence of the initial velocity model. The initial v S models were perturbed using a normal random distribution with a standard deviation of 0.1, and 200 v S models were produced. The Rayleigh dispersion curves were inverted using 200 perturbed initial models. Selected v S models for the station pair BNDS-  Fig. 8. The parameter setting did not change. It can be observed that the histograms of inverted results tend to have a near normal distribution that congregates at the mean value with some deviation. The distributions reflect how well v S is constrained at each depth (Fig. 8).

Rayleigh wave tomography
A 2-D tomographic inversion technique was applied for the calculation of the group velocity variations derived from the dispersion measurements of Rayleigh waves from the one-year cross-correlations recorded from January 1, 2009 Fig. 7 a S-wave velocity model for two station pairs GHIR-ZHSF and KRBR-ZHSF. The starting v S model is shown by red lines. The blue line is the final v S model. b Sensitivity kernels of group velocity as a function of period for the v S velocity measurements. c Dispersion curves related to the two station pairs mentioned above with their predicted dispersion curves (syn.) and their measured dispersion curves. The paths for GHIR-ZHSF and KRBR-ZHSF are depicted in Fig. 1b Earthq Sci to January 1, 2010. Before conducting ANT, we checked the data availability of all broadband stations inside and around Iran and finally decided to use data from January 1, 2009 to January 1, 2010 because of its continuous data from most stations existed during this time period.
Fast marching surface wave tomography (FMST), the iterative non-linear inversion package developed by Rawlinson (2005) and Rawlinson and Sambridge (2005), was used for analysis. This method includes the forward calculation and inversion procedures. Subspace method was used in the inversion, which is based on an assumption of local linearity that can reduce the perturbation of the model parameters and then make the inversion result approach to the current model. The inversion step allows both smoothing and damping regularization to limit the non-uniqueness of the solution. The Fast Marching Method (FMM) is a grid-based numerical algorithm based on the eikonal equation that is formulated to determine the first arrival phase of surface waves rather than the group time. However, to describe the dissipation of the group energy, an eikonal solver can be used if multi-pathing is not included. In this case, the interfering waves cause the group energy to follow notably different paths. Therefore, if the phase and group velocities have a similar geographic pattern, comparable results can be obtained (Arroucau et al. 2010;Saygin and Kennett 2010;Young et al. 2011;Saygin and Kennett 2012). Young et al. (2011) obtained similar group and phase velocity maps using FMM in southeastern Australia. Generally, the FMM method is an iterative nonlinear approach with an assumption of local linearity in the inversion step. The non-linear relationship between the travel time and the group velocity could have been explained by repeated applications of FMM and subspace inversions (Rawlinson 2005;Rawlinson and Sambridge 2005). The damping value was estimated based on a trade-off curve between the data misfit and model roughness. Combination of the FMM for calculation of the forward problem and the subspace method for inversion provides tomographic imaging. Fig. 8 The influence of the initial velocity model on the non-linear iterative damped least squares inversion was investigated using the normal randomly distribution with standard deviation 0.1 at 200 runs. Selected v S models for station pair BNDS-GHIR and at 200 runs and at depths of a 10 km, b 20 km, c 30 km, d 40 km, e 60 km, and f 80 km are shown. The path for BNDS-GHIR is depicted in Fig. 1b The potential resolution of the tomographic results was evaluated using the checkerboard synthetic tests for 16, 20, 30, and 40 s periods that used actual path distribution. Two sets of tests with pattern sizes of 2 9 2 and 1 9 1 were conducted. The general recoveries for both tests were good; however, higher resolution for central Iran was attained for the 1 9 1 cell size (see the Appendix Fig. 12). The maximum error was 5 % noise signal with a perturbing velocity of 2.8 ± 0.1 km/s and superimposed alternating high-and low-velocity anomalies (as shown in the Appendix, Fig. 12). The number of iterations in the subspace inversion depends on the frequency value, which varies from 2 for the higher frequencies to 5 for the lower frequencies. The checkerboard test results of the observed data and the ray-path coverage for periods of 16, 20, 30, and 40 s are shown in Fig. 9a-d, respectively. The checkerboard results suggest that the resolution is fairly good in most periods in central Iran, indicating that the pattern and absolute amplitude values were recovered well. However, for the eastern and southeastern locations of the region, due to rare distribution of stations, the path coverage is not dense, and most waves travel in parallel; therefore, the resolution is so low that smearing effects are apparent in the eastern region, which make the inversion process not recover the absolute amplitudes well. The tomography maps at  [4.0, during 1974-2008, is from EHB catalog (Engdahl et al. 1998) different periods indicate the general features of the structure at different depths is influenced by their sensitivity kernels (e.g., Yang et al. 2007;Huang et al. 2010;Tibuleac et al. 2011) (Fig. 9). To guide the interpretation, the sensitivity kernels for different periods were used (Fig. 7b). The Rayleigh wave with the shortest period of 16 s has fair sensitivity to the top 10 km, and the wave with the longest period of 40 s has peak sensitivity at *60 km depth and fair sensitivity up to *80 km. Thus, the use of the dispersion curves from the 16 to 40 s periods allows us to constrain the shear velocities at the depth range of 10 to 80 km. Each tomography map is illustrated with the corresponding seismicity in each period (Fig. 9). Events with focal depths between 15 and 30 km were plotted on the 16 and 20 s group velocity maps, and events with focal depths between 30 and 55 km were plotted on 30 and 40 s group velocity maps.

Discussion
The level of seismicity in Makran is low and increases from west to east. The western and eastern portions of the subduction zone of Makran have different seismic and tectonic characteristics (Byrne et al. 1992;Zarifi 2006). Three N-S profiles (AA 0 , BB 0 , and CC 0 ; profile locations are  Fig. 9a). The topography along the profiles is also shown in the top panels. The black dots are the earthquakes (with magnitude C4) from EHB catalog (Engdahl et al. 1998) occurred along the profiles with width 0.5 in each side shown in Fig. 9a) are illustrated in Fig. 10, in which the focal depth distribution is shown. The western portion of Makran has an abnormally low level of deep seismicity. The intermediate depth seismicity characteristics related to western Makran within the downgoing plate are different from the dominant shallower seismicity of the Zagros region (Fig. 10). Across the Sistan Suture Zone, this seismicity pattern changes to a low seismicity condition compared with the Zagros region. Recent and active deformation in Sistan is dominated by right-lateral strikeslip and thrust faults related to the indentation of Iran by the Arabian shield (Berberian et al. 2000). The seismic activity in the mountain ranges, including the Taftan-Bazman volcanic arc, is quite weak (Fig. 9). In 1979, several right-lateral moderate-sized earthquakes occurred between the Lut and Helmand blocks, and inside these two blocks, there is little seismicity. This seismic activity increases the possibility that the Sistan Suture Zone plays a role in the segmentation between eastern and western Makran, and therefore, the continuity of this structure could be defined as a boundary between western and eastern Makran (Byrne et al. 1992). To the east, the distance of the volcanic arc and fore-arc setting increases, and this suggests that the slab dips shallower as it moves eastward (Byrne et al. 1992;Zarifi 2006;Shad Manaman et al. 2011). The eastern portion of Makran has relatively lower dips compared with the western portion (Zarifi 2006) (see Fig. 10b, c). Eastern Makran has experienced most of its seismic activity near Chaman and Ornach-Nal Faults (Zarifi 2006).
According to our tomographic maps, shallow earthquakes are expected to occur in the location of high-velocity anomaly, where Arabian Plate begins to subduct beneath the Straits of Hormuz (Fig. 9a, b). This high-velocity anomaly far beneath the Straits of Hormuz emerges more clearly at shorter periods, reflecting the thin crust under the Oman sea floor and Makran fore-arc setting (25-30 km). Most of the earthquakes that occur in this region are expected to be shallow, and as a confirmation of our outcome, nearly all of the seismicity associated with this region occurs at depths of \30 km (Jackson and McKenzie 1984), see Fig. 10a. Within the downgoing plate towards the north, where the low-velocity anomaly is located, we expect events would occur at intermediate depths due to the down-dip elongation of the subducting slab (see Fig. 10b, c). The deeper events occur along the downgoing slab where the subducting plate bends below the Taftan-Bazman volcanic arc (see Fig. 10b, c). The deepest earthquakes of the Makran region concentrate near the Taftan volcano due to the accommodation of the final component of the motion between Arabia and Eurasia (Byrne et al. 1992). We investigated the crustal thickness using the latest Moho map obtained for the same area using a different approach and the data reported by Shad Manaman et al. (2011). For better accuracy in the analysis, a high-resolution version of the Moho map in Shad Manaman et al. (2011) was used, as illustrated in Fig. 11.
The 16 and 20 s maps are sensitive to the upper crust at *25 to 30-km depth based on the sensitivity curves in Fig. 7b. A low-velocity anomaly in the Oman Sea south of the Makran region is observable (Fig. 9a). Although low-velocity sediments with 6-to 7-km thickness in front of the Makran deformation front (White 1982;Fowler et al. 1985) might not be visible at this period, the Rayleigh waves of the 16 s period are still likely to sample the accretionary prism. At the same time, these waves might be influenced by the uppermost section of the subducted crust. The northward subducting Arabian Plate is shown as the high-velocity anomaly along the Straits of Hormuz in the 16 and 20 s map (see Fig. 9a, b). This high-velocity anomaly reflects the influence of subduction of the descending slab that is older, denser, and colder than the continental crust next to it. The Strait of Hormuz is considered a transition between the Zagros collision and the Makran oceanic subduction (Regard et al. 2010). A sharp transition boundary between the low-and high-velocity zones with a northwest trend is clearly depicted at the Minab fault system (see Fig. 9a, b). This transition indicates the boundary between the Zagros region with high seismicity in the northwest and the Makran region with low seismicity to the east. The earthquake distribution surrounding the Minab fault is restricted to the west of the Jaz Murian depression (Fig. 9a, b). The Jaz Murian depression has been interpreted as a fore-arc basin (Farhoudi and Karig 1977). The large trench-arc distance (*400-600 km) suggests that the angle of subduction is notably low and is consistent with thermal modeling in this region (e.g., Smith et al. 2013). An aseismic region extends from the deformation front for nearly 200 km in the western Makran (Fig. 9a, b). Most of the shallow seismicity is related to the Zagros Mountain. A pronounced low-velocity anomaly extends to the SW-NE in the east of Minab fault, which is attributable to volcanic arc and back-arc settings of the Makran region and the Bazman and Taftan volcanoes (see Fig. 9a, b). This low-velocity anomaly has its origin in thicker crust caused by a warm lithosphere wedge overlying the subducting Arabian Plate and might be coincident with the source of the volcanic magmas. The subducting slab could cause partial melting when volatiles are released and rise into the overlying upper mantle. This effect is consistent with previous seismic observations in this region (e.g., Shad Manaman et al. 2011). Another low-velocity anomaly is observable at the Sultan volcanic setting; however, this anomaly occurs at the edges of the area with acceptable resolution (see Fig. 9a). The trend of low-velocity anomaly indirectly suggests geological and geophysical evidence for the geometry of slab. A transition from low to high velocity is observable in central Makran between the Sistan Suture Zone and the Lut block. Although Byrne et al. (1992) assumed that this suture zone separates the Lut and Helmand blocks, our results show that this suture does not appear to segment different blocks.
The 30 and 40 s maps are most sensitive to a depth of *25-55 km based on the sensitivity curves in Fig. 7b; although these maps are of lower resolution than others, resolution within the area is still reasonable. The 40 s tomographic map displays its maximum sensitivities at a depth of *55 km, as Fig. 7b illustrates. The low-velocity anomalies beneath the volcanic arc on the maps are similar to those at 16 and 20 s and reveal that the crustal thickness below the Taftan volcano is at least 50 km deep, which is compatible with the latest Moho Map obtained for the same area using a different approach and data by Shad Manaman et al. (2011) (Fig. 9d). The Moho map was produced using a surface wave tomography method to image the S-velocity structure of the upper mantle and Moho depth. The earthquake distributions surrounding these volcanoes and illustrated in the 30 and 40 s maps (Fig. 9c, d) indicate that the number of intermediate earthquakes increases with the downgoing plate. These earthquakes have normal faulting focal mechanisms (Jackson and McKenzie 1984;Laane and Chen 1989) with predominantly down-dipping T-axes, which indicate that the subducted slab is in tension. The high-velocity anomaly compatible with the descending slab moves northward in these two maps, which demonstrates the northward direction of Arabian Plate subduction. Because few studies exist on the deep structure of the upper mantle in the Makran region, the use of ANT within this area provides new images on the crust and uppermost mantle. Results obtained from recent studies in this region have poor resolution compared with the results of the checkerboard test with input models consisting of anomalies of 4 spike size (e.g., Shad Manaman et al. 2011); however, in this study, the results demonstrate the higher resolution of *1 9 1 for central Iran and limited resolution on the order of 2 9 2 towards the southeastern region of the study area. Many of the prominent features in our results are consistent with the known geological structures. The lateral resolution of the tomographic maps obtained by ANT greatly depends on various parameters, including the path coverage and inter-station distances. In this study, for the period range of 10-50 s, we retained only paths longer than 250 km ([3 wavelengths at 10 s). Therefore, we selected a grid spacing of *110 km. Additionally, ANT is most powerful when ambient noise exists over a broad azimuthal range. Thus, we demonstrated that the fraction of high SNR paths in all azimuths are greater 50 % at four period bands of 10-20 s, 20-30 s, 30-40 s, and 40-50 s and were approximately 0.53, 0.64, 0.69, and 0.51, respectively. Consequently, the distribution of the useful amount of ambient noise is sufficient for use in tomography. Fig 7a shows the v S velocity models obtained for two different inter-station distances at the study area. The sensitivity kernels as a function of period are also presented in Fig. 7b. The v S velocity models presented in Fig. 7 demonstrate significant agreement with tomography results according to the sensitivity curves (see Fig. 7b). The v S begins to increase gradually with depth. Comparing the v S structures of the two station pairs, the most obvious v S discrepancy is the low velocity of KRBR-ZHSF at a depth of 18-30 km. This low velocity is quite similar to the tomographic results along the path of KRBR-ZHSF at 16-20 s (Fig. 9a, b), which are sensitive down to a depth of *20-30 km, according to the sensitivity curves in Fig. 7b.

Conclusions
The following results were obtained from this study: (1) In this research, we explained how sufficient noise energy can be recorded in the Makran region at periods of 10-50 s from which the empirical Green's functions were extracted. It was also shown that the Rayleigh wave Green's functions were extracted by computing cross-correlations between records using observations collected over 12 months at pairs of seismic stations. (2) Our results exhibit seasonal variability in the study area. This seasonal variation indicates that the Green's functions reconstructed using cross-correlation can show differences in quality during the summer and winter. Although we showed that coherent Rayleigh wave signals exist in all periods and most azimuths across the Makran region, they are sufficiently isotropically distributed in azimuth to deliver accuracy in dispersion measurements if integrated over a long time period, such as a year.
(3) In conclusion, various resolution tests showed that our data and methods are sufficient to provide highresolution tomographic images of surface wave group velocities in the region. The group velocity maps show low-velocity anomalies beneath the volcanic arc correlated with large Moho depths, a finding compatible with the Moho Map obtained by Shad Manaman et al. (2011) (Fig. 11). The high-velocity anomalies along the Straits of Hormuz in the group velocity maps determine the northward subducting Arabian Plate. Finally, these group velocity maps are significantly improved in lateral resolution over those of previous studies, which have relied on traditional earthquake-based surface wave tomography.
Acknowledgments Seismicity of each period map is plotted from EHB catalog (Engdahl et al. 1998). The seismic data used in this study were obtained from the GSN/IRIS Global Seismographic Network (http://www.iris.edu/hq/programs/gsn), the Virtual European Broadband seismic Network (VEBSN) and international Institute of Earthquake Engineering and Seismology (IIEES). Many of the figures in this paper were prepared using GMT (Wessel and Smith 1998;www.soest.hawaii.edu/gmt, last accessed April 2015). We would like to thank Dr. N. Shad Manaman (from Sahand University of Technology, Iran) for generously providing us with his research results on the Moho depth and also Dr. N. Mirzaei (from University of Tehran, Iran) for his valuable suggestions in this study. A special note of thanks goes to Dr. T. Shirzad (from Islamic Azad University, Tehran, Iran) for his help and support. We also appreciate the research council of Tehran University for their support of this research.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix
See Fig. 12. Fig. 12 Input checkerboard test models with velocity perturbation of about 2.8 ± 0.3 km/s and grid size 2 9 2 (a) and 1 9 1 (c), and the corresponding recovered models are given in (b) and (d) Earthq Sci