Shallow crustal radial anisotropy beneath the Tehran basin of Iran from seismic ambient noise tomography

Physics of the Earth and Planetary Interiors 231 (2014) 16–29

Contents lists available at ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Shallow crustal radial anisotropy beneath the Tehran basin of Iran from seismic ambient noise tomography Taghi Shirzad a, Z. Hossein Shomali b,c,⇑ a

Department of Physics, Damavand Branch, Islamic Azad University, 39715-194 Damavand, Iran Institute of Geophysics, University of Tehran, 14155-6466 Tehran, Iran c Department of Earth Sciences, Uppsala University, 75236 Uppsala, Sweden b

a r t i c l e

i n f o

Article history: Received 14 September 2013 Received in revised form 5 April 2014 Accepted 6 April 2014 Available online 2 May 2014 Keywords: Tehran basin Ambient seismic noise Tomography Shear wave velocity Radial anisotropy

a b s t r a c t We studied the shear wave velocity structure and radial anisotropy beneath the Tehran basin by analyzing the Rayleigh wave and Love wave empirical Green’s functions obtained from cross-correlation of seismic ambient noise. Approximately 199 inter-station Rayleigh and Love wave empirical Green’s functions with sufficient signal-to-noise ratios extracted from 30 stations with various sensor types were used for phase velocity dispersion analysis of periods ranging from 1 to 7 s using an image transformation analysis technique. Dispersion curves extracted from the phase velocity maps were inverted based on non-linear damped least squares inversion method to obtain a quasi-3D model of crustal shear wave velocities. The data used in this study provide an unprecedented opportunity to resolve the spatial distribution of radial anisotropy within the uppermost crust beneath the Tehran basin. The quasi-3D shear wave velocity model obtained in this analysis delineates several distinct low- and high-velocity zones that are generally separated by geological boundaries. High-shear-velocity zones are located primarily around the mountain ranges and extend to depths of 2.0 km, while the low-shear-velocity zone is located near regions with sedimentary layers. In the shallow subsurface, our results indicate strong radial anisotropy with negative magnitude (VSV > VSH) primarily associated with thick sedimentary deposits, reflecting vertical alignment of cracks. With increasing depth, the magnitude of the radial anisotropy shifts from predominantly negative (less than 10%) to predominantly positive (greater than 5%). Our results show a distinct change in radial anisotropy between the uppermost sedimentary layer and the bedrock. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Seismic anisotropy in the Earth’s crust provides new insight on past and present deformation of the crust (Cochran et al., 2006; Guo et al., 2012). The mechanisms controlling crustal anisotropy are debated among seismologists; however, two potential processes that are often used to explain crustal anisotropy are the alignment of microcracks due to the presence of a stress field in the uppermost crust and the alignment of large rock or mineral fabrics in the middle to lower crust (e.g., Stein and Wysession, 2003; Schubert, 2007). In the shallower part of the crust, the primary sources of the anisotropy are often thought to be pervasive wavelength-scale deformation and the presence of microcracks, which are often fluid-saturated, in places where pressures are ⇑ Corresponding author at: Department of Earth Sciences, Uppsala University, 75236 Uppsala, Sweden. Tel.: +46 18 471 25 19; fax: +46 18 50 11 10. E-mail addresses: [email protected] (T. Shirzad), [email protected], Hossein. [email protected] (Z.H. Shomali). http://dx.doi.org/10.1016/j.pepi.2014.04.001 0031-9201/Ó 2014 Elsevier B.V. All rights reserved.

low enough to allow them to remain open (Mainprice and Nicolas, 1989). Previous studies using shear wave analysis (e.g., Tadokoro et al., 1999; Crampin et al., 2004; Peng and Ben-Zion, 2004) showed that the fastest shear-wave splitting directions are oriented toward crack alignments and that the delay time in shear-wave splitting is proportional to the density of the cracks. An anisotropic layered structure with a vertical symmetric axis and hexagonal symmetry is conventionally introduced to account for the period-dependent discrepancy between the speeds of horizontally (Love) and vertically (Rayleigh) propagating waves (Dziewonski and Anderson, 1981; Xie et al., 2013). Therefore, shear-wave splitting resulting from the surface wave data can provide detailed information on the radial anisotropy at various depths and even on sub-crustal structures (Fry et al., 2010; Endrun et al., 2011; Liu and Niu, 2012). However, short-period surface waves, which are primarily sensitive to the shear wave structure and thickness of the crust, are difficult to obtain from earthquake waveforms. This is because short-period teleseismic surface waves are always diffracted and are therefore associated

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

with multi-pathing. Recently, empirical Green’s functions (EGFs) were extracted from ambient seismic noise collected by a dense seismic array. Shapiro and Campillo (2004) stated that the dispersion curves obtained from inter-station measurements of the noise field are similar in character to those estimated from earthquakes. Thus, seismic ambient noise tomography (ANT hereafter) provides a new way to study shallow crustal heterogeneity at a higher resolution than that of traditional surface wave tomography methods. ANT has been applied to various regions around the world, primarily to map lateral heterogeneity in the crust and mantle using vertical (e.g., Shirzad et al., 2013; Ward et al., 2013) and transverse components (e.g., Lin et al., 2008; Bensen et al., 2008). Recently, seismic ambient noise was used to image radial seismic anisotropy at regional (e.g., Huang et al., 2010; Moschetti et al., 2007; Guo et al., 2012) and global (e.g., Fry et al., 2010; Pawlak et al., 2012) scales. The main points we address in this study are as follows: (1) obtain the quasi-3D shear wave velocity, (2) attempt to constrain the radial anisotropy, and (3) explain the physical cause of the anisotropy under the Tehran sedimentary basin. Tehran, the capital of Iran, is located in a tectonically active region situated between the Central Alborz Mountain belt and the Central Iran plateau (see Fig. 1). Because of the low seismicity and the high population of the study area, traditional passive earthquake tomographic or active seismic exploration methods are not well suited for determining the shallow crustal heterogeneity. The lateral resolution of regional seismic tomographic methods is also limited to approximately 60–100 km (e.g., Shad Manaman et al., 2011) in the study area. Previous works found the prominent crustal and upper mantle anisotropy in the study area; e.g., Sadidkhouy et al. (2008) and Sadidkhouy et al. (2009), who used shear-wave splitting of P-to-S converted phases on receiver functions and SKS phases. Most previous studies of anisotropy are restricted by the lateral and vertical resolutions in the shallow crustal structures (upper 15 km). The Tehran basin is characterized by soft sediments in the uppermost kilometers that overly more rigid bedrock. This leads to a large velocity contrast for high- and low-frequency surface waves, which are primarily sensitive to sub-surface and deeper structures, respectively. Such structures are prone to resonance phenomena. Detailed knowledge of the shallow (upper 3 km) crustal structure of the Tehran basin is essential for a variety of reasons, including improving seismic hazard mitigation policies by using more realistic strong ground motion simulations. In this study, phase velocity dispersion curves of the Rayleigh and Love fundamental modes in a broad period range (1–7 s) were calculated using inter-station measurements of the seismic noise field to determine the radial anisotropy in the uppermost crust. For this purpose, we extracted EGFs of the Rayleigh and Love waves from ambient noise in the study area following the common low-frequency technique developed by Bensen et al. (2007) and Lin et al. (2008) and applied the RMS-stacking method developed by Shirzad and Shomali (2013). Subsequently, to obtain a quasi-3D model of shear velocity, we applied two techniques. First, we performed a 2D tomographic inversion using the observed EGFs to obtain Rayleigh and Love phase velocity maps at each period, which we then used to extract Rayleigh and Love synthetic phase velocity dispersion curves at each evenly spaced geographic grid point. Second, the Rayleigh and Love synthetic phase velocity dispersion curves at each grid point were separately inverted to obtain the local 1D shear wave velocity; these data were mapped onto the original grid to construct the quasi-3D shear wave velocity. Finally, shear-wave polarization revealed the lateral and vertical distribution of radial anisotropy beneath the Tehran basin area in the upper 3 km of the crust.

17

2. Data and processing methodology 2.1. Dataset In this study, we analyzed three-component recordings of continuous data from 30 stations in the study area depicted in Fig. 1. The dataset originates from 18 digital narrow-band seismometers (SH) with SS1 seismometer sensors (corner frequency P1 Hz) operated by the Iranian Seismological Center (IRSC) at the University of Tehran, 10 digital accelerometers (HG) with CMG-5T sensors manufactured by the Guralp company and operated by the Tehran Disaster Mitigation and Management Organization (TDMMO), and two digital broadband instruments (BH) with a CMG-3T sensor manufactured by the Guralp company and operated by the International Institute of Earthquake Engineering and Seismology (IIEES). For the TDMMO acceleration network, continuous data spanning 2009 October 1 and 2011 May 25 were used, whereas for the IRSC and IIEES seismic networks only continuous data from 2010 were analyzed.

2.2. Preprocessing and extracting EGFs In general, we followed the time-domain procedure introduced by Bensen et al. (2007) and Lin et al. (2008) to extract Rayleigh and Love wave EGFs, respectively. We note that the processing steps required to analyze the ambient noise depend on various parameters, including the quality of data and its frequency content. Thus, we cannot define a single set of instructions to optimally process data acquired under different conditions. The horizontal components (E, N) were first rotated into the radial (R) and transverse (T) components for each station pair. These were corrected for trend and mean, and a 5-point, zero-phase band-pass Butterworth filter was applied over a period range of 1–7 s. The data were subsequently down-sampled to 10 Hz, and one-bit normalization and spectral whitening were performed on the remaining data to reduce the effects of transient signals and to broaden the available frequency band. To calculate the cross-correlation, the data were split into 10-min segregated time segments, then all possible combinations of the station pairs with a match sensor type were calculated (Cho et al., 2007; Shirzad et al., 2013). We cross-correlated waveforms with the matched sensor types and therefore did not correct for the instrumental response, primarily due to the ambiguity in the instrumental response of the short-period sensors deployed within the IRSC network. In addition, as reported by Tibuleac et al. (2011), the level of coherency of the noise crosscorrelation function (NCF) degrades for non-similar sensor types. Previous studies (e.g., Roux et al., 2005; Sabra et al., 2005) suggest that it is more appropriate to take the time-derivative of the NCF to get the EGF. However, the higher frequency content of the EGFs is enhanced by the negative time-derivative. Therefore, we applied the negative time-derivative to the NCFs before stacking to extract suitable EGFs. We followed the RMS-stacking method initially introduced by Picozzi et al. (2009) for high-frequency data and developed by Shirzad and Shomali (2013), which should always increase the root mean square (RMS) of the obtained signals within the signal window. The number of NCFs (threshold value) used in the RMS-stacking process is a function of the change in the gradient of the RMS-curve, as shown in Fig. 2a. The threshold values of all components among all inter-station pairs varied between 2000 and 3000. The threshold value is a function of various parameters, including inter-station distances, frequency content and signal coherency among all stations. In general, stabilization of the EGFs is reached faster for low frequencies and for tangential components (see Fig. 2a). The strongest lag between causal and acausal

18

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

Fig. 1. (a) The locations of the various seismic networks used in this study: IRSC (white triangles), IIEES (inverted white triangles), seismic and TDMMO (black triangles) accelerometer networks. (b) The geological maps of the study area from Shafiee and Azadi (2007). The known active faults are indicated by solid black lines. The black circles (points P1 and P2) show the positions referred to in Fig. 7. The two dashed lines show the location of the anisotropy cross-sections shown in Fig. 11. The Tehran region and province are marked by black polygons in (a) and (b).

components was used to measure Rayleigh and Love wave dispersion phase velocity curves (Yang and Ritzwoller, 2008). Fig. 2b shows an example of the resulting EGFs between the station pair AFJDMV. There is a clear difference between the

arrival times of the EGFs on the transverse–transverse (TT), the radial–radial (RR) and the vertical–vertical (ZZ) cross-correlation functions. A fast Love wave arrival is seen in the transverse– transverse (T–T) cross-correlation, whereas similar EGFs seen in

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

19

Fig. 2. (a) RMS curve of transverse (T–T; black circles) and vertical (Z–Z; gray circles) components of the station pair AFJDMV. (b) The extracted Green’s functions of station pair AFJDMV at the period band of 3–5 s. TT, RR and ZZ are the transverse, radial and vertical components, respectively.

the Z–Z and R–R cross-correlations are dominated by Rayleigh waves. The Rayleigh wave EGFs observed in the Z–Z signal are larger in amplitude than those observed in the RR signal. Fig. 3 shows all Rayleigh and Love wave EGFs in the 3–5 s period band with signal-to-noise ratio (SNR) greater than 15 at stations from the TDMMO acceleration network. The SNR is defined as the ratio of the peak within the signal window (proportional to the 0.6–1.5 km/s velocity window) and root-mean square (RMS) noise in the noise window (proportional to the 0.2–0.4 km/s velocity window). Linear trends corresponding to the Rayleigh and Love wave packages are denoted on the Z–Z and the T–T cross-correlations as depicted in Fig. 3, with apparent phase velocities of the order of 0.95 and 1.10 km/s, respectively. 2.3. Data selection To stabilize the tomographic results, we applied several data selection criteria prior to performing the tomographic analysis, following Bensen et al. (2008). First, we applied an inter-station constraint of a minimum three wavelengths (3k) to the EGFs to stabilize the measurements at shorter distances and to satisfy the farfield assumption of the cross-correlation method. Second, we applied a period-dependent SNR, which is defined as the ratio of the peak in the signal window to the RMS noise in the noise window filtered with a specific central period. All EGFs of the vertical and traversal components with SNR larger than 15 were preserved

for further processing. By applying this constraint, the number of station pairs corresponding to Rayleigh and Love wave signals was reduced by 4 and 2, respectively. The third data selection criterion is based on the variety of the phase velocity on temporally seasonal subsets of the data. To extract seasonal EGFs, we stacked 3-month cross-correlation functions using the RMS stacking method with a 2-month overlap period. Our analysis found that the standard deviations of the dispersion measurements of both Rayleigh and Love waves in the segregated seasonal subsets of data were almost in the range of the 12-month phase velocity measurement error bars. The final data selection criterion is based on the tomographic residuals. To stabilize the tomographic results, we rejected data with residuals larger than three standard deviations (3r) in the observed phase velocity for both Rayleigh and Love waves at each period. The number of rejected signals for each data selection criteria is listed in Table 1. 2.4. Directionality A fully diffusive noise field or perfectly azimuthal noise sources are usually assumed in ANT applications. This assumption implies that the strength and frequency content of the ambient seismic noise observed on causal and acausal lags of EGFs are equal, or, in other words, the signal is symmetric. Any deviation by the two-sided EGFs from this symmetric behavior may cause a bias in the dispersion measurements and consequently in the

Fig. 3. (a) Love and (b) Rayleigh wave empirical Green’s functions of the TDMMO data in the 3–5 s period band with SNR >15. The dashed lines indicate 1.10 and 0.95 km/s in (a) and (b), respectively.

20

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

Table 1 The number of inter-station paths rejected by the data selection criterion at each period and the input background velocity and damping regularization parameter of the 2D tomography for the vertical (Z–Z) and horizontal (T–T) components as a function of period. Period (s)

Data selection criterion

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 7.0

Tomographic information 3r

3k

Average velocity (m/s)

Damping

ZZ

TT

ZZ

TT

ZZ

TT

ZZ

TT

3 5 6 9 8 12 14 20 20 18 25 23

5 4 4 3 8 14 14 33 33 34 33 27

5 4 7 3 5 5 9 10 12 16 19 23

3 5 6 2 5 5 5 8 12 10 14 14

850 858 862 865 871 877 885 894 905 920 932 950

970 963 965 977 985 992 1002 1010 1015 1030 1032 1075

4.5 3.0 1.7 1.6 1.5 1.26 1.3 1.28 1.25 0.9 0.9 0.7

5.5 4.0 3.5 2.0 2.0 1.5 0.95 0.9 0.7 0.7 0.65 0.58

Fig. 4. (a) The azimuthal distribution of inter-station paths calculated for the 20° azimuthal bin. The directional dependence of high SNR (>15) Rayleigh (b) and Love (c) wave empirical Green’s functions (EGFs) signals plotted in the period band of 1–7 s. The average fraction of empirical Green’s functions was calculated by normalizing the number of path with a SNR >15 in a given 20° azimuthal bin to the total number of paths in that given bin.

tomographic results. We therefore studied the strength and azimuthal distribution of the ambient seismic noise in the Tehran region. The azimuthal distribution of inter-station paths calculated for the 20° azimuthal bin, shown in Fig. 4a, indicates that the dominant inter-station direction is ENE–WSW. Following Bensen et al. (2008), we used the average fraction of the EGFs to limit the influence of array-specific geometry. The average fraction of Rayleigh and Love EGFs were calculated by normalizing the number of paths with SNR >15 in a given 20° azimuthal bin to the total number of paths in that bin. Fig. 4(b) and (c) shows the azimuthal distribution of Rayleigh and Love wave EGFs with SNR >15 in the period band of 1–7 s. Note that the SNR was calculated for each lag separately and that Fig. 4b and c include two-sided EGFs for Rayleigh and Love waves. The average fraction of Rayleigh and Love EGFs with SNR >15 is approximately 0.70 and 0.81, respectively. Note that even though the coherency level of the ambient seismic noise is more pronounced for the Love TT EGFs, it is at least 45% in all period bands used in this study for the Rayleigh (ZZ and RR; Fig. 4b) EGFs; similarly, the coherency is always greater than 50% for the Love (TT; Fig. 4c) EGFs. Consequently, a useful coherency level of Rayleigh and Love wave signals exists in ambient seismic noise at most azimuths in the Tehran basin, and this information can be used for further analysis.

2.5. Phase velocity measurements Following Yang and Ritzwoller (2008), we used the EGFs-lag with the largest amplitude to calculate the phase velocity

Fig. 5. Example of selected Love (black line) and Rayleigh (gray line) phase dispersion velocity curves between station pairs D01–D20 (solid line) and D15–D22 (dashed line).

dispersion curves of the Rayleigh and Love waves separately. The uncertainties associated with the phase velocity dispersion measurements are smaller than those associated with the group velocity dispersion curves (e.g., Bensen et al., 2008; Lin et al., 2008); thus, in this study, phase velocities were estimated using a modified far-field approximation of an image transformation analysis technique developed by Yao et al. (2006). The phase velocities were measured using band-pass filtered EGFs in the central

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

21

Fig. 6. (a) Percentage of Love and Rayleigh wave phase velocity variance reduction as a function of period. (b) The RMS final travel-time residuals in millisecond as a function of period. (c) The RMS final velocity residual in km/s as a function of period.

period of 1–7 s with an interval of 0.5 s. The lower bound of the filter was fixed to 1.0/(T + nT) and the upper bound to 1.0/(T  nT), where T is the central period and n was fixed to 0.2 (personal communication). The phase velocity dispersion curves of the Rayleigh and the Love waves for station pairs D01–D20 and D15–D22 running N–S and NW–SE across the Tehran basin are shown as examples in Fig. 5. 2.6. 2D tomography We constructed 2D phase velocity distribution maps for each period from 1 to 7 s using the measured Rayleigh and Love wave phase velocity dispersions and a non-linear 2D tomographic inversion technique (Rawlinson, 2005; Rawlinson and Sambridge, 2005). The inversion algorithm was performed using a fast-marching method (FMM; Rawlinson and Sambridge, 2005) that tracks rays for the initial and updated velocity distributions before every step of the inversion. This accounts for propagation effects caused

by rapid changes in the velocity field. Thus, by an iterative application of FMM using a subspace inversion method, non-linearities between the travel time and velocity are accounted (Rawlinson, 2005; Rawlinson and Sambridge, 2005). For each period, an average of the observed phase velocities satisfying the 3r condition was calculated and used as the starting model for the inversion. We also fixed the spatial smoothness operator to 0.9 km. The optimal damping parameter was determined from standard L-curve analysis. The average phase velocities as well as the damping parameters for each period are listed in Table 1. The study area was divided into 0.050°  0.050° geographical grids for periods from 1.0 to 4.5 s and 0.075°  0.075° geographical grids for longer periods, i.e., 5–7 s. To obtain 2D phase velocity maps with evenly spaced geographic grid points (0.075°  0.075°) for all periods, Bspline interpolation was applied separately to the phase velocity maps of shorter periods (1.0–4.5 s). The variance reductions, data residuals and model parameter residuals for the full period range (1–7 s) are shown in Fig. 6. In general, the higher variance

22

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

reduction at shorter periods is due to a greater number of paths, whereas the larger residuals for longer periods indicate greater scatter and more limited number of ray paths in these periods (see Fig. 6). After obtaining the 2D Rayleigh and Love phase velocity maps, we extracted the phase velocity dispersion curves for 0.075°  0.075° geographical grids in the period range 1–7 s to construct synthetic phase velocity dispersion curves. The extracted synthetic phase velocity dispersion curves at the grid points P1 (35.70° N, 51.22° E) and P2 (35.72° N, 51.49° E) are shown as an example in Fig. 7(a) and (b), respectively. 3. Results 3.1. Shear wave velocity model The synthetic phase velocity dispersion curves at each evenly spaced geographical grid point (0.075°  0.075°) were inverted using an iterative damped least-squares inversion procedure called the surf96 program (Herrmann and Ammon, 2002) to produce a 1D shear wave velocity model at each grid point. The initial VS-model was parameterized with multiple layers of constant thickness (0.1 km) over a half-space with depths of 0.0–3.0 km. The initial VS-model contained no low-velocity layer and was constructed using the average VS model obtained by Shabani et al. (2008) for the study area. The P-wave velocity was updated by the Vp/Vs ratio of the initial model, which was fixed to the value reported by Shabani et al. (2008). Finally, these 1D shear velocity models were

inserted into the original geographic grid points to construct a quasi-3D shear wave velocity model. The inversion procedure was repeated separately for the synthetic Rayleigh and Love wave dispersion curves. The final results of the inversions for VSV and VSH at the points P1 and P2 are shown in Fig. 7c and d, respectively, as an example. The difference between VSV and VSH at these two points indicates a significant amount of radial anisotropy at depths less than 3.0 km. 3.2. Uncertainty of shear wave velocities In this paper a stochastic measure of the model uncertainties is used using parametric test to determine the influence of the initial model and the uncertainties in dispersion measurements. The parametric tests guarantee the robustness and stability of the tomographic results. Thus, we performed parametric tests to determine the influence of the initial VS-model and regularization parameters, i.e., the weighting factor and smoothness parameter, on the non-linear iterative damped least squares inversion. In the first test, we perturbed the initial VS-model 1000 times using a normal random distribution with a standard deviation of 0.1. We then inverted the synthetic Rayleigh and Love dispersion curves using these perturbed initial models 1000 times. The results of the inversion initialized with different VS-models at 1000 m depth at the grid points P1 and P2 are shown in Fig. 8(a) and (b), respectively. The results are shown as a function of depth in Fig. A(a) and (b) in Appendix A. The distribution of the calculated shear wave velocities is normal, as shown in Fig. 8(a) and (b). Note that in these

Fig. 7. Synthetic phase velocity dispersions (solid line) and synthetic phase velocity dispersions from VS (dote line) at grid points P1 (35.70° N, 51.22° E) (a) and P2 (35.72° N, 51.49° E) (b). (c) and (d) indicate the corresponding inverted shear wave velocity (left panel) and radial anisotropy (right panel). Synthetic Love and Rayleigh dispersion curves are shown by the black and gray solid lines, respectively. Synthetic Love and Rayleigh dispersion curves from VSH and VSV are shown by the black and gray dotted lines, respectively.

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

results, the errors associated with VSV are generally smaller than VSH errors (see Fig. 8a and b). The second parametric test consisted of applying the bootstrapping method to calculate the standard error associated with the 1000 hypothetical test models constructed from the observed Rayleigh and Love phase velocity dispersion curves (e.g., Adam and Lebedev, 2012; Guo et al., 2012). In this part, the observed Rayleigh and Love phase velocity measurements were perturbed by 5%. The 1D shear wave velocity calculation step, which includes 2D tomography, interpolation of phase velocity to a 0.075°  0.075° geographical grid, extraction of synthetic phase velocity dispersion curves and 1D shear wave velocity inversion, was repeated for the 1000 sets separately using the perturbed Rayleigh and Love dispersion measurements. We applied the same 2D phase velocity tomography regularization parameters (e.g., damping, smoothness parameters) as we used for the observed data. The results are shown in Fig. 8(c) and (d) as histograms at 1000 m depth and are also shown as functions of depth in Fig. A(c) and (d) in Appendix A. The distribution of the calculated shear wave velocity is Gaussian, as depicted in Fig. 8(c) and (d). As before, VSV has smaller errors than VSH does. Spatial distribution of shear wave velocity standard deviation is depicted by white contours in Fig. 9 at different depth levels. In order to calculate the spatial distribution of shear wave velocity model uncertainties, the initial VS-model was firstly perturbed for 1000 times using a normal random distribution with a standard

23

deviation of 0.1 km/s. Then synthetic phase velocity dispersion curves were inverted using an iterative damped least-squares inversion procedure at each evenly spaced geographical grid point. As shown in Fig. 9, lateral variation of model uncertainties is in order of 0.04–0.1 km/s, and increases slightly with depth. Moreover lateral velocity anomalies are, on average, approximately ten times of the corresponding standard deviation values. 3.3. Radial anisotropy In general, radial anisotropy is controlled by five parameters: VPH, VPV, VSH, VSV and the ellipticity g (Dziewonski and Anderson, 1981; Becker et al., 2008). Thus, differences between VSV and VSH at any depth can be attributed to the presence of radial anisotropy at that depth. We used the information available in the calculated VSV (bottom maps of Fig. 9) and VSH (top maps of Fig. 9) to determine the distribution of radial anisotropy with depth beneath the Tehran basin. We used the formula given by Yanovskaya et al. (1998) to address the presence of radial anisotropy: c = 2(VSH  VSV)/(VSH + VSV)  100, where c is the radial anisotropy coefficient. Maps showing the distribution of radial anisotropy estimated from the discrepancy between Rayleigh and Love waves at depths of 500, 900, 1250, 1600, 2000 and 2450 m are shown in Fig. 10. Analysis of the calculated anisotropy coefficients c shows that in the uppermost crustal layers, e.g., 500 m, the anisotropy is dominated by the relatively high VSV (c < 0, VSV > VSH); however,

Fig. 8. The left panels of (a) and (b) indicate the influence of the initial VS model on the obtained VSH and VSV models for the two grid points P1 (35.70° N, 51.22° E) and P2 (35.72° N, 51.49° E) at 1000 m depth, respectively. The left panels of (c) and (d) indicate the influence of the regularization parameter on the obtained VSH and VSV models for the two grid points P1 (35.70° N, 51.22° E) and P2 (35.72° N, 51.49° E) at 1000 m depth, respectively. The corresponding radial anisotropy is shown in the right panels.

24

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

Fig. 9. Horizontal slices through the VSH model at depths of 0.5 km (a), 1.25 km (b) and 2.0 km (c). (d), (e) and (f) show horizontal slices through the VSV model at these same depths. Gray contours indicate the calculated shear wave velocities in km/s. The estimated model uncertainties of VSH and VSV models derived from the parametric test are also depicted by white contours. The Tehran region see Fig. 1 is marked by blue polygon, and also the known faults are indicated by solid black lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

at increasing depths, it gradually becomes dominated by the relatively high VSH (c > 0, VSV < VSH; Fig. 10). Two radial anisotropy profiles are also shown in Fig. 11 (refer to Fig. 1b for their positions). Prominent large-scale anisotropic zones with positive magnitude (c > 0, VSV < VSH) can be observed in the mountain belts and at greater depths. The sedimentary basin is, in general, dominated by a layer with negative radial anisotropy (c < 0, VSV < VSH).

4. Discussion Incorporating seismic anisotropy into velocity models is a necessary step for obtaining a more realistic parameterization of the Earth and a better understanding of geodynamic development. In the Earth’s crust, anisotropy may be caused by preferentially aligned microcracks, by layered bedding in sedimentary formations, or by highly foliated metamorphic rocks (Crampin and Chastin, 2003). Surface waves generated by earthquakes are often used to identify radial anisotropy in the crust and mantle. Because of the uneven distribution of seismic sources and seismograph stations, the intrinsic attenuation and scattering along the ray paths and the (usually) unknown characteristics of the sources, the resolution of traditional tomographic methods is poor, especially at short periods (<20 s), which are primarily sensitive to structures in the upper crust. However, seismic ambient noise provides a valuable way to obtain short-period surface wave EGFs between station pairs. ANT not only allows higher lateral resolution but also provides a broader frequency content, with resolution as good as 1 s, as presented in this study. In this study, observed phase velocity dispersion curves extracted from ANT were applied to

determine the vertical distribution of anisotropic fabric beneath the Tehran basin. Our region of interest is characterized by soft sediments in the uppermost 3.0 km. Although it is well known that sedimentary rocks can be seismically anisotropic, there have been few detailed investigations of the underlying cause of such anisotropy. A strong data selection criterion was applied to ensure that only the most stable EGFs were used in our analysis. The robustness of the tomographic inversion was also improved significantly by posing a smoothness regularization parameter. Consequently, the lateral resolution of the 2D phase velocity maps and the amplitudes of their anomalies are underestimated, and in the following only a qualitative discussion of the anisotropic maps as function of depth is presented. Generally speaking lateral resolution of the 2D tomographic Rayleigh/Love results are mainly controlled by ray coverage. A general view of available inter-station ray coverage at four selected periods, 1, 3, 5 and 7 s, for the 2-D tomography inversion is provided in Appendix B. At all depths shown in Fig. 9, standard deviation of the calculated shear wave velocities (depicted by white contours) is lowest near the northwestern edges of the study area, which is due to spars ray coverage (see Appendix B). At depth of 1250 m (see Fig. 9b and e), the average values of the standard deviations are smaller than the corresponding values for other depths. At 2000 m depth (Fig. 9c and f), the average values of the standard deviations are larger than shallower depths and are more uniform, although the lateral variation of the standard deviations changes slightly at this depth due to poorer sensitivity at greater depths. Below 2900 m depth, the model is very poorly constrained. Previous studies (Berberian and Yeats, 1999; Abbasi and Shabanian, 1999) showed that the Tehran region can be divided

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

25

Fig. 10. Horizontal 2D maps showing the magnitude of radial anisotropy (c) at depths of 500, 900, 1250, 1600, 2000 and 2450 m. The corresponding model uncertainties for magnitude of radial anisotropy (%) are shown by white contours.

into three main geological units: (1) the northern Central Alborz Mountain belt, (2) the eastern Sepayeh Mountain belt, and (3) the Tehran basin, which includes central, southern and western Quaternary alluvial deposits (see Fig. 1b). The Tehran basin is also surrounded by many major active faults, such as the North Tehran Fault, which is the most prominent tectonic structure in the immediate vicinity of the city. Improved knowledge of the shallow crustal structure of the Tehran basin provided by a quasi-3D absolute shear wave velocity model is of fundamental importance for understanding the geodynamic evolution of this large and important tectonic structure. The VS maps provide a first approximation of shallow crustal properties that can be used to calculate seismic amplifications and seismic hazards. The VS models may also benefit to estimate rapid assessment of earthquake damage (i.e., Wald et al., 1999). The distribution of shear wave velocities, especially VSV, delineates several distinct low- and high-velocity zones separated mainly by geological boundaries (see Fig. 9). The lowest shear wave velocities are located at shallow depths, i.e., 500 m, where low-velocity sedimentary alluvial deposits clearly outline the triangular shape of the Tehran basin (Fig. 9d). In the Tehran basin, the low-velocity zone (VSH and VSV are both approximately 850 m/s) deepens toward the south-southwest at approximately 1250 m depth, which shows the presence of thicker sediments in the southern part of the basin than in the north (Fig. 9b and e). Previous studies independently found that the thickness of the sedimentary layer is approximately 1.0–1.5 km in the Tehran basin (Jafari, 2002; Hamzehloo et al., 2007). At a depth of approximately 2.0 km, the North Tehran Fault separates a strong high-velocity anomaly associated with the Central Alborz Mountain belt from the low-velocity layers associated with the Tehran basin (see VSH in Fig. 9c and VSV in Fig. 9f). The Sepayeh Mountain belt in the east-

ern part of the study region appears to lack a deep high-velocity root, unlike the Central Alborz Mountain belt (Fig. 9c and f), whose high-velocity root has been reported previously (Jafari, 2002). The results obtained in this study show very low velocity (<875 m/s) for the sedimentary layers and high velocity (1150 m/s) for the bedrock in the Tehran basin. No low-velocity zone was detected in our analysis. Although resolving the layering of anisotropy in fine detail would require an investigation of the trade-offs of various anisotropic and isotropic parameters, our radial anisotropic results (Fig. 10) display patterns that clearly reveal the presence of radial anisotropy in the uppermost 3.0 km of the crust beneath the Tehran basin. As shown in Fig. 10, the Quaternary alluvial deposits of the Tehran basin are resolved by prominent radial anisotropy with a negative magnitude (c < 0, VSV > VSH) at all depth intervals, reflecting vertical alignment of cracks. However, for the northern Central Alborz Mountain belt, the sense of the radial anisotropy changes from negative to positive at approximately 2000 m depth. We also observe a major anomaly with a VSV < VSH signature at 500–1250 m depth resolved at approximately 35.66° N, 51.45° E bounded by the North Tehran Fault to the north and the North Rey Fault to the south. In contrast to the Central Alborz Mountain belt, which exhibits positive radial anisotropy at depths greater than approximately 2000 m, the eastern Sepayeh Mountain belt is characterized by negative anisotropy at all depths. An interesting finding of this study is the separation between the thick sedimentary deposits of the Tehran basin at a depth of 0.4– 1.5 km and the bedrock, identified from a change in magnitude of the radial anisotropy (Fig. 11a). According to the available geological, geo-technical and borehole data from the study area (e.g., Abbasi and Shabanian, 1999; Hamzehloo et al., 2007), the bedrock deepens from 0.4 km to 1.4 km across the Tehran basin in the

26

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

Fig. 11. Vertical cross-sections through the radial anisotropy model in North–South (a) and East–West (b) directions. The cross-section positions are shown in Fig. 1b. White contours indicate the estimated model uncertainties for magnitude of radial anisotropy (%) as a function of depth.

north–south direction. This is similar to what our results resolve, as shown in profile A–A0 in Fig. 11a. Comparison of Fig. 11(a) and (b) indicates that the variation in bedrock depth beneath the Tehran area along the E–W profile (200 m) is much less than that along the N–S profile (1000 m). In fact, the bedrock is nearly flat in the E–W direction. Unlike checkerboard tests, the parametric tests carried out in this research provide a stochastic measure of model uncertainties. Part of the apparent radial anisotropy resolved at different parts of the model could be caused by the different resolutions of Rayleigh and Love waves i.e., sparse ray-coverage. Although, the errors associated with VSV are generally smaller than VSH errors (see Fig. 8), however, model uncertainties associated with the calculated shear wave velocities, VSV and VSH, show similar pattern for well constrained parts of model e.g., depths less than 500 m (see Fig. 9). Therefore the risk that part of apparent radial anisotropy in the Tehran region might have been caused by different resolutions of VSH and VSV is not significant. A qualitative discussion of the tomographic maps presented in this paper was also guided by model uncertainties associated to VSV and VSH and the corresponding uncertainties for magnitude of radial anisotropy as depicted in Figs. 10 and 11. Improved knowledge of the shear velocity structure of the Tehran basin is of fundamental importance for understanding the geodynamic evolution of this large and important tectonic structure. Aligned cracks, either open or filled with a different material, are an important influence on velocity at shallow depths (less than

3.0 km) in the crust and are the most likely cause of the radial anisotropy observed beneath and bordering the Tehran basin. Anisotropy with horizontally polarized shear waves propagating faster than vertically polarized ones (VSH > VSV) is observed in the lower crust and mantle lithosphere and indicates horizontally oriented fabric. Anisotropy with VSV > VSH is observed in the upper crust in a number of locations and suggests the occurrence of vertically oriented fabric. The fastest direction of seismic-wave propagation aligns with the direction of maximum compressive stress where the seismic anisotropy is predominant and where vertical cracks are oriented with the ambient stress field (Crampin, 1978).

5. Conclusions (1) Our results show that a useful level of coherent noise energy exists at periods 1–7 s in the Tehran basin, even under high local noise and poor path coverage that limits the resolution in this region. Quasi-3D shear wave velocity and seismic radial anisotropy (transverse isotropy with a vertical symmetry axis) models were obtained over the Tehran basin by analyzing the Rayleigh and Love wave EGFs recovered from seismic ambient noise recordings. The average Rayleigh and Love wave phase velocity curves for the Tehran basin were inverted to investigate the presence of radial anisotropy in the shallow sub-crustal structures. The results show clear lateral variation in the upper crust VS structure in the study

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

region. Strong signals of low VSV appear in the central part of the Tehran basin, which signifies the relatively narrow triangular shape of the basin. (2) Shear wave velocity polarization, VSV and VSH, resolved using Rayleigh and Love waves reveals the presence of significant radial anisotropy in the upper 3 km of the crust beneath the Tehran basin. (3) The results and the best-fitting models show evidence of strong (up to 20%) negative (VSV > VSH) upper crustal anisotropy and a slight but significant positive anisotropy in the bedrock. Radial anisotropy in the uppermost sedimentary layer is only poorly resolved, but the data suggest that VSV > VSH. (4) From the shear wave velocity results, it can be concluded that the depth of the bedrock varies between 400 and 1400 m from north to south beneath the Tehran basin with an average VS of approximately 1150 m/s. An interesting finding of this study is the change in magnitude of the radial anisotropy (Fig. 11a), which suggests a clear separation between a sedimentary layer overlying rigid bedrock.

27

study would not have been possible without the digital data sets provided by TDMMO. We would also like to thank Dr. M. A. Riahi and M. Etedal (University of Tehran, Iran) for their constructive comments and useful suggestions. We appreciate Dr. Mallory Yoang from Research School of Earth Sciences, the Australian National University is thanked for her help and support during the calculation of phase velocity. We thank Dr. Nick Rawlinson (ANU) for making his Fast Marching Surface Tomography Package available to us. The editor, V. F. Cormier and an anonymous reviewer are specially thanked for their invaluable comments during the review process. The second author appreciates the research council of Tehran University for their support of this research. Some plots were made using the Generic Mapping Tools (GMT) version 4.2.1 (Wessel and Smith, 1998; www.soest.hawaii.edu/ gmt, last accessed March 2014). Appendix A See Fig. A.

Acknowledgments Appendix B We express our sincere thanks to IRSC, IIEES seismic and TDMMO accelerometer networks for providing us the data. This

See Fig. B.

Fig. A. (a) and (b) indicate the influence of the initial velocity on the obtained 1D shear wave velocity. The initial velocity was perturbed 1000 times using a normal distribution with a standard deviation of 0.1. (c) and (d) show the effects of different inversion parameters on the obtained shear wave velocity. The corresponding radial anisotropy is shown in the right panels as a function of depth. The reference velocity model is depicted by solid black lines.

28

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29

Fig. B. Ray-paths coverage for Rayleigh (upper panels) and Love (lower panels) waves phase velocity maps at four selected period: (a and e) 1, (b and f) 3, (c and g) 5 and (d and h) 7 s. The black triangles indicate the seismic stations and the Tehran region is depicted by blue polygons. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

References Abbasi, M.R., Shabanian, E., 1999. Fault Map of North Tehran. Publication by International Institute of Earthquake Engineering and Seismology, Tehran, Iran. Adam, J.M.-C., Lebedev, S., 2012. Azimuthal anisotropy beneath southern Africa from very broad-band surface-wave dispersion measurements. Geophys. J. Int. 191, 155–174. Becker, T.W., Kustowski, B., Ekström, G., 2008. Radial seismic anisotropy as a constraint for upper mantle rheology. Earth Planet. Sci. Lett. 267, 213–227. http://dx.doi.org/10.1016/j.epsl.2007.11.038. Bensen, G.D., Ritzwoller, M.H., Shapiro, N.M., 2008. Broadband ambient noise surface wave tomography across the United States. J. Geophys. Res. 113, B05306. http://dx.doi.org/10.1029/2007JB005248. Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M., Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys. J. Int. 169, 1239–1260. Berberian, M., Yeats, R.S., 1999. Pattern of historical earthquake rupture in the Iranian plateau. Bull. Seism. Soc. Am. 89 (1), 120–139. Cho, K.H., Herrmann, R.B., Ammon, C.J., Lee, K., 2007. Imaging the upper crust of the Korean Peninsula by surface-wave tomography. Bull. Seism. Soc. Am. 97, 198– 207. Cochran, E.S., Li, Y.-G., Vidale, J.E., 2006. Anisotropy in the shallow crust observed around the San Andreas fault before and after the 2004 M 6.0 Parkfield earthquake. Bull. Seism. Soc. Am. 96, S364–S375. http://dx.doi.org/10.1785/ 0120050804. Crampin, S., Chastin, S., 2003. A review of shear-wave splitting in the crack-critical crust. Geophys. J. Int. 155, 221–240. Crampin, S., 1978. Seismic wave propagation through a cracked solid: polarization as a possible dilatancy diagnostic. Geophys. J. R. Astr. Soc. 53, 467–496. Crampin, S., Peacock, S., Gao, Y., Chastin, S., 2004. The scatter of time-delays in shear-wave splitting above small earthquakes. Geophys. J. Int. 156, 39–44. Dziewonski, A., Anderson, D., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297–356. Endrun, B., Lebedev, S., Meier, T., Tirel, C., Freiderich, W., 2011. Complex layered deformation within the Aegean crust and mantle revealed by seismic anisotropy. Nat. Geosci. 4, 203–207. Fry, B., Deschamps, F., Kissling, E., Stehly, L., Giardini, D., 2010. Layered azimuthal anisotropy of Rayleigh wave phase velocities in the European Alpine lithosphere inferred from ambient noise. Earth Planet. Sci. Let. 297, 95–102. Guo, Z., Xing, G., Wang, W., Yao, Z., 2012. Upper- and mid-crustal radial anisotropy beneath the central Himalaya and southern Tibet from seismic ambient noise tomography. Geophys. J. Int. 189, 1169–1182. http://dx.doi.org/10.1111/j.1365246X.2012.05425.x. Hamzehloo, H., Vaccari, F., Panza, G.F., 2007. Toward a reliable seismic microzonation in Tehran, Iran. Eng. Geol. 93, 1–16. Herrmann, R.B., Ammon, C.J., 2002. Computer Programs in Seismology-Surface Waves, Receiver Functions and Crustal Structure. Saint Louis University http:// www.eas.slu.edu/People/RBHerrmann/ComputerPrograms.html. Huang, H., Yao, H., van der Hilst, R.D., 2010. Radial anisotropy in the crust of SE Tibet and SW China from ambient noise interferometry. Geophys. Res. Let.. http:// dx.doi.org/10.1029/2010GL044981.

Jafari, M.K., 2002. Seismic Geotechnical Microzonation of South-East of Tehran. International Institute of Earthquake Engineering and Seismology, Tehran, Iran (in Farsi). Lin, F.-C., Moschetti, M.P., Ritzwoller, M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps. Geophys. J. Int. 173, 281–298. Liu, H., Niu, F., 2012. Estimating crustal seismic anisotropy with a joint analysis of radial and transverse receiver function data. Geophys. J. Int. 188, 144–164. http://dx.doi.org/10.1111/j.1365-246X.2011.05249.x. Mainprice, D.M., Nicolas, A., 1989. Development of shape and lattice preferred orientations: application to the seismic anisotropy of the lower crust. J. Struct. Geol. 11, 175–189. Moschetti, M.P., Ritzwoller, M.H., Shapiro, N.M., 2007. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh wave group velocity maps. Geochem. Geophys. Geosyst. 8, Q08010. http://dx.doi.org/ 10.1029/2007GC001655. Pawlak, A., Eaton, D.W., Darbyshire, F., Lebedev, S., Bastow, I.D., 2012. Crustal anisotropy beneath Hudson Bay from ambient noise tomography: evidence for post-orogenic lower-crustal flow. J. Geophys. Res. 117. http://dx.doi.org/ 10.1029/2011JB009066. Peng, Z., Ben-Zion, Y., 2004. Systematic analysis of crustal anisotropy along the Karadere-Düzce branch of the North Anatolian fault. Geophys. J. Int. 159, 253– 274. Picozzi, M., Parolai, S., Bindi, D., Strollo, A., 2009. Characterization of shallow geology by high-frequency seismic noise tomography. Geophys. J. Int. 176, 164– 174. http://dx.doi.org/10.1111/j.1365-246X.2008.03966.x. Rawlinson, N., Sambridge, M., 2005. The fast marching method: an effective tool for tomographic imaging and tracking multiple phases in complex layered media. Explor. Geophys. 36, 341–350. Rawlinson, N., 2005. FMST: Fast Marching Surface Tomography Package. Research School of Earth Sciences, Australian National University, Canberra ACT 0200. Roux, P., Sabra, K., Gerstoft, P., Kuperman, W.A., Fehler, M.C., 2005. P-waves from cross correlation of seismic noise. Geophys. Res. Let. 32, L19303. http:// dx.doi.org/10.1029/2005GL023803. Sabra, K.G., Gerstoft, P., Roux, P., Kuperman, W.A., Fehler, M.C., 2005. Extracting time-domain green’s function estimates from ambient seismic noise. Geophys. Res. Lett. 32, L03310. http://dx.doi.org/10.1029/2004GL021862. Sadidkhouy, A., Javan-Doloei, Gh., Siahkoohi, H.R., 2008. Seismic anisotropy in the crust and upper mantle of the Central Alborz region, Iran. Tectonophysics 456, 194–205. Sadidkhouy, A., Javan-Doloei, Gh., Siahkoohi, H.R., 2009. Mantle anisotropy in the Central Alborz obtained from SKS analysis. Iran. J. Geophys. 2 (2), 1–12. Schubert, G., 2007. Treatise on Geophysics. Elsevier, Amsterdam, Boston. Shabani, E., Cornou, C., Hagshenas, M., Wathelet, M., Bard, P.Y., Mirzaei, N., & Eskandari-Ghadi, M., 2008. Estimating shear-waves velocity structure by using array methods (FK and SPAC) and inversion of ellipticity curves at a site in south of Tehran, The 14th Word Conference on Earthquake Engineering, October 12– 17, 2008, Beijing, China. Shad Manaman, N.S., Shomali, H., Koyi, H., 2011. New constraints on upper-mantle S-velocity structure and crustal thickness of the Iranian plateau using partitioned waveform inversion. Geophys. J. Int. 184, 247–267. http:// dx.doi.org/10.1111/j.1365-246X.2010.04822.x.

T. Shirzad, Z.H. Shomali / Physics of the Earth and Planetary Interiors 231 (2014) 16–29 Shafiee, A., Azadi, A., 2007. Shear-wave velocity characteristics of geological units throughout Tehran city, Iran. J. Asian Earth Sci. 29, 105–115. Shapiro, N.M., Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett. 31, L07614. http:// dx.doi.org/10.1029/2004GL019491. Shirzad, T., Shomali, Z.H., Riahi, M.-A., 2013. An application of ambient noise and earthquake tomography in the Rigan area southeast of Iran. Seis. Res. Let. 84 (6), 1014–1020. http://dx.doi.org/10.1785/0220130044. Shirzad, T., Shomali, Z.H., 2013. Shallow crustal structures of the Tehran basin in Iran resolved by ambient noise tomography. Geophys. J. Int. 196, 1162–1176. http://dx.doi.org/10.1093/gji/ggt449. Stein, S., Wysession, M., 2003. An Introduction to Seismology, Earthquakes and Earth Structure. Blackwell Publishing, Oxford, ISBN 0-86542-078-5. Tadokoro, K., Ando, M., Umeda, Y., 1999. S-wave splitting in the aftershock region of the 1995 Hyogo-ken Nanbu earthquake. J. Geophys. Res. 104, 981–991. Tibuleac, I.M., von Seggern, D.H., Anderson, J.G., Louie, J.N., 2011. Computing green’s functions from ambient noise recorded by accelerometers and analog, broadband, and narrow-band seismometers. Seis. Res. Let. 82, 661–675. Wald, D.J., Quitoriano, V., Heaton, T.H., Kanamori, H., Scrivener, C.W., Worden, C.B., 1999. TriNet ‘‘ShakeMaps’’: rapid generation of instrumental ground motion

29

and intensity maps for earthquakes in southern California. Earthquake Spectra 15, 537–556. Ward, K.M., Porter, R.C., Zandt, G., Beck, S.L., Wagner, L.S., Minaya, E., Tavera, H., 2013. Ambient noise tomography across the Central Andes. Geophys. J. Int. 194, 1559–1573. http://dx.doi.org/10.1093/gji/ggt166. Wessel, P., Smith, W.H.F., 1998. New, improved version of the generic mapping tools released. EOS Trans. Am. Geophys. Un. 79, 579. Xie, J., Ritzwoller, M.H., Shen, W., Yang, Y., Zheng, Y., Zhou, L., 2013. Crustal radial anisotropy across Eastern Tibet and the Western Yangtze Craton. J. Geophys. Res. Solid Earth 118. http://dx.doi.org/10.1002/jgrb.50296. Yang, Y., Ritzwoller, M.H., 2008. Characteristics of ambient seismic noise as a source for surface wave tomography. Geochem. Geophys. Geosyst. 9, Q02008. http:// dx.doi.org/10.1029/2007GC001814. Yanovskaya, T.B., Kizima, E.S., Antonova, L.M., 1998. Structure of the crust in the Black Sea and adjoining regions from surface wave data. J. Seismol. 2, 303–316. Yao, H., Van Der Hilst, R.D., De Hoop, M.V., 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis – I. Phase velocity maps. Geophys. J. Int. 166, 732–744. http://dx.doi.org/10.1111/j.1365246X.2006.03028.x.