- Email: [email protected]
Crust and uppermost mantle structure of Tehran region from analysis of teleseismic P-waveform receiver functions
Tectonophysics 364 (2003) 115 – 133 www.elsevier.com/locate/tecto
Crust and uppermost mantle structure of Tehran region from analysis of teleseismic P-waveform receiver functions Javan Doloei a,b,*, Roland Roberts a a b
Department of Earth Sciences and Geophysics, Uppsala University, Villava¨gen 16, Uppsala S-752 36, Sweden International Institute of Earthquake Engineering and Seismology (IIEES), P.O. Box 19395/3913, Tehran, Iran Received 7 August 2001; accepted 3 February 2003
Abstract Teleseismic P-waveform receiver function analysis has been used to estimate the main features of the crust and uppermost mantle structure of the Tehran region. In this study, 20 teleseismic earthquakes that were recorded at the seven stations of Iranian Long Period Array (ILPA) were used to calculate receiver functions. The radial receiver functions are analysed through inverse modelling and forward modelling. The results of modelling show a good match between observed radial receiver functions and synthetic radial receiver functions up to 30 s. Based on this study, the crustal structure of this area is characterised by three main layers: The upper crust has a P-wave velocity between 4 and 5.8 km/s and a 14-km thickness. The results representing the uppermost crustal structure indicate a strong dependency on the lithology around each station. The middle crust has a positive P-wave velocity gradient from 6 to 6.4 km/s down to f 30 km depth. A P-wave velocity gradient from 6.4 to 7.5 km/s characterises the lower crust with Moho at 46 F 2 km at this area. A low-velocity mantle, with velocity gradually increasing with depth down to f 60 km depth, has been revealed beneath the Tehran region. D 2003 Elsevier Science B.V. All rights reserved. Keywords: Crust; Iran; Receiver function; Uppermost mantle structure; Seismology; Tehran
1. Introduction The major crustal features and the seismic velocities of the uppermost mantle are important constraints in understanding the tectonic evolution of an area. Pwaveform teleseismic receiver function analysis is becoming a common method to estimate this information. In this study, the main features of the crust and
* Corresponding author. International Institute of Earthquake Engineering and Seismology (IIEES), P.O. Box 19395/3913, Tehran, Iran. E-mail address: [email protected] (J. Doloei).
uppermost mantle structure beneath the Iranian Long Period Array (ILPA) are investigated using the receiver function technique. The ILPA network consists of seven stations, which are located southwest of the metropolis of Tehran. This area is situated at the southern border of the central Alborz mountains, which are the branches of the Alpine – Himalayan orogeny in Iran (Fig. 1, top). Iran lies within the continental collision zone between Africa and Eurasia. Active faulting, recent volcanism, and high surface elevations are important characteristics of this area. Geomorphic and trenching investigations (De Martini et al., 1998) along the Kahrizak fault, which is
0040-1951/03/$ - see front matter D 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0040-1951(03)00049-0
116
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
117
Table 1 Information on the 20 teleseismic earthquakes recorded at the ILPA network Event number
Date
Origin time (hhmmss)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1977/2/16 1977/2/18 1977/2/19 1977/3/2 1977/3/8 1977/3/18 1977/3/19 1977/3/19 1977/5/12 1977/5/28 1977/9/19 1977/9/23 1977/10/16 1977/11/18 1977/12/23 1978/1/13 1978/1/14 1978/6/14 1978/7/23 1978/12/23
104020.90 205129.80 223404.10 095323.20 231728 214352.40 105625.10 193508 111753.10 055148.10 051309.20 055755.60 210917.70 101741 210207.50 200303.80 032439 113420 144236.90 112312
Latitude (j) 0.50 33.07 53.57 6.77 0.45 16.77 44.20 16.81 39.27 1.73 1.98 11.21 9.73 4.35 39.13 44.69 34.81 38.32 22.28 23.25
Longitude (j)
Depth (km)
Magnitude
GCARC (j)
Backazimuth (j)
125.98 140.82 170.03 123.74 100.02 122.33 148.20 122.35 117.71 120.52 126.63 118.22 117.12 102.02 143.16 149.70 139.26 142.36 121.51 122.07
33 42 33 52 22 37 70 39 22 54 33 33 33 33 19 52 14 40 17 33
6.1 6.0 6.7 6.1 6.0 7.0 6.0 5.8 5.8 6.0 6.2 6.0 6.0 5.9 6.0 6.0 6.7 6.3 7.4 7.2
78 72 77 72 57 65 71 65 52 74 79 78 77 62 70 71 70 70 62 62
98 62 32 94 116 86 48 86 64 103 100 112 112 118 55 48 61 56 81 80
The GCARC and backazimuths are related to station IR1.
located south of Tehran and west of the ILPA network, provided some new constraints on the activity of this region. Previous studies of seismic structure of Iran give us a general view of the crust and uppermost mantle across the Iranian plateau. The earliest study of Hedayati et al. (1976) utilised a microearthquake survey to show a 31-km crustal thickness below the Tehran region. A general study of seismic P-wave and S-wave velocities by using arrival time data of regional earthquakes (Chen et al., 1980) provided seismic structure information beneath Iran with average Pn velocity of 8.0 F 0.1 km/s. Surface wave analysis of a few events (Asudeh, 1982) suggested a one-layer model crust with 45 km thickness beneath the eastern part of Alborz mountain range. Asudeh (1982) showed that the crust and upper mantle in the eastern part of Iran are characterised by low body wave velocities. The gravity field and crustal structure of Iran, in the study by Dehghani and Makris (1984),
showed that the Bouguer gravity along the Alborz mountain range varies between 100 and 120 mgal and that the crustal thickness is less than 35 km. That is, these mountains have no root. The first and only deep seismic soundings in Iran were conducted in 1978 between the Lut Block in the east of the country and the Zagros mountain range in the west. The results of this study identified a 40-km crustal thickness for the central Iranian plateau (Giese et al., 1984). Numerical modelling of the deformation of the Iranian plateau (Sobouti and Arkani-Hamed, 1996) identified that crustal thickening and shear deformation tend to localise near the southwest of the country. Their results specified a 45-km crustal thickness along the Alborz mountain range. Based on receiver function analysis of teleseismic data, Mangino and Priestley (1998) provided more information on the crustal structure of the southern Caspian region in the north of the Alborz mountain range. To give an
Fig. 1. (Top) Map of major structural divisions of Iran. Rectangle = research area (bottom). (Bottom) Location map of stations (triangles) for ILPA network in the Tehran area. Q = unfolded Quaternary cover alluvium; Qpl = Paleocene – Eocene, deposited on pre-Paleocene folded substratum, folded by middle and late Alpine events, conglomerate, limestone, and volcanic rocks; A3 OM = sediments deposited on pre-Oligocene substratum, folded during late Alpine event; Ev = extrusive volcanic rocks, mainly andesites and pyroclastics, acidic and intermediate basic.
118
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
improved picture of the Iranian part of the tectonically active Alpine – Himalayan orogenic belt, a receiver function study of available teleseismic ILPA data is necessary and useful. The purpose of this paper is to determine the crustal and uppermost mantle seismic velocity structures of the Tehran area based on a receiver function analysis of 20 teleseismic earthquakes, which occurred in the range of 52– 79j from centre of the ILPA network with magnitude greater than 5.8. Table 1 gives information on the teleseismic earthquakes used. The information of teleseismic earthquakes of Table 1 is extracted from PDE catalogue of US Geological Survey.
2. Receiver function analysis The time domain receiver function technique has been described in detail (Langston, 1977, 1979; Owens, 1987; Owens et al., 1984; Ammon et al., 1990) for use with single-station three-component data. Deconvolving the vertical seismogram from the radial component gives us the radial receiver function. A similar deconvolution of the vertical component from the transverse component provides the tangential receiver function, which is useful to test the data regarding the assumptions of a laterally homogeneous, isotropic Earth structure. It is well known that a teleseismic P-wave incident to a velocity discontinuity (e.g., the Moho) will generally produce P-waves and S-waves [both transmitted (Pp, Ps) and Moho reverberation phases such as PpPms and PpSms + PsPms]. Because the S-wave will travel more slowly than the Pwave, the time difference in the arrival of the two phases Pp and Ps is a direct measure of the depth of the original interface. The time interval between the Ps converted phase and the Moho PpPms reverberation phase gives the ratio of P-wave and S-wave velocities (Zandt et al., 1995). The seismic crustal and uppermost mantle structure of the Tehran area is investigated using teleseismic Pwaveforms recorded at the ILPA network for 2 years (1977 – 1978). The continuous waveform database was obtained from the Data Management Centre (DMC) of the Incorporated Research Institution for Seismology (IRIS). The ILPA network consists of seven stations IR1 –IR7. Six stations (IR2 –IR7) are placed around a
circle with a radius of about 30 km, and one station (i.e., IR1) is located in the centre (Fig. 1, bottom). This configuration allows investigation of the nonuniqueness of the deduced velocity models through comparing and combining the results of receiver function inversions from all stations. Due to the low sampling rate of Long Period data (1 Hz), a time window with a length of 180 s of Pwaveforms is used to calculate true amplitude receiver functions (Ammon, 1991). To separate actual teleseismic P-waveforms from noise and tapering effects, all time windows of three-component seismograms are selected 60 s before and 120 s after P-wave onset. The maximum frequency content of this type of data is 0.5 Hz, which is sufficient for resolving the main features of the crust and uppermost mantle structure around the three-component single station. The mean and a trend are removed from the selected time window of threecomponent observations and a 20% cosine taper is applied to the signal. The gains of horizontal and vertical components of borehole sensor are the same at each station of ILPA network. Therefore, it was not necessary to correct instrument gains. The true amplitude receiver function technique of Ammon (1991) is used to calculate the radial and tangential P-waveform receiver functions. To stabilise deconvolution, a water level parameter of c = 0.000001 is used. This value was selected after various trials in the range c = 1.0 10 8 to 1.0 100. The value chosen produces acceptable noise levels in most of our receiver functions. The averaging function is defined by a deconvolution of the vertical component from itself at each water level c to find out the narrowest possible Gaussian curve. Due to the lowlevel noise in our observations, a small value of water level gives minimum distortion effects. Applying a Gaussian filter (Ammon, 1991) controls the frequency content of the receiver functions. A Gaussian filter width parameter a = 1.015 is used, which lowpasses the waveforms at f 0.48 Hz. As an example, Fig. 2 shows a three-component seismogram of a teleseismic earthquake (December 23, 1978), with a blow-up of the P-waveform and corresponding receiver functions. This event was recorded at station IR1 with a surface wave magnitude of 7.2 and an epicentral distance of 62j. Although the amplitude of PpPms phases is weak due to surface topography, fortunately, a few teleseismic earthquakes in this study
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
119
Fig. 2. (a) Three-component seismogram of the earthquake of December 23, 1978 that was recorded at station IR1 in BAZ = 80j. Origin time is fixed at zero time so the time of different phases shows the travel time of each phase from source to receiver. (b) Windowing P-waveform; a lag window of the P-waveform with a length of 180 s is selected by taking 60 s of the signal before and 120 s after P-onset. (c) The corresponding radial and tangential receiver functions. The Pp phase is plotted with a 10-s delay.
clearly show PpPms phases on radial receiver functions at the other stations.
3. Velocity structure of Tehran area from teleseismic receiver functions The ILPA network is situated in a complex tectonic setting in northcentral Iran, surrounded in the west, north, and northeast by the active Alpine – Himalayan orogeny (i.e., the Alborz mountain range; Fig. 1, bottom). The ILPA network is so configured so that two stations IR1 and IR6 are located on an extrusive volcanic zone trending NW – SE. In the north, stations IR2 and IR7 are also located on the volcanic rocks on the border of the Qazvin depression. Station IR3 is located in the northeast of the ILPA network, surrounded in the west and north by the recent volcanic rocks, and in the east and south by the unfolded Quaternary cover alluviums and synorogenic con-
glomerate sediments. Stations IR4 and IR5 are placed on the Eocene – Oligocene Karaj formation in the south of the ILPA network (Fig. 1, bottom). Generally, the Cenozoic sediments on the southern side of the Alborz mountain range constitute a heterogeneous discontinuous package of polymicts siliciclastics with subordinate carbonate layers in the lower part with clear synorogenic characteristics. The lower part of the succession consists of lensoidal (up to 1500 m thick) elongated masses of thick-bedded calcarenites (Alavi, 1996). The Neogene– Quaternary sediments (including the Hezar-Darreh and Kahrizak Formations) show every aspect of sedimentation and include coarse to fine, northward-coarsening, polymict fanglomerates, and breccias with reworked sediments of the Neogene age that are mainly derived from the Alborz magmatic assemblage. The Alborz magmatic assemblage (including the Karaj Formation of Eocene age) consists of both submarine and subaerial, porphyritic and nonporphyritic, massive lava flows of andesitic and
120
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
andesitic – basaltic compositions; extensive well bedded; and calcareous and noncalcareous pyroclastics such as tuffs and agglomerates (andesitic and trachyandesitic in composition). Volcanic breccias (widely distributed volcaniclastic sandstones) and red, coarse, polymict conglomerates (which are intruded by diorites, monzodiorites, and diabase dikes and sills) are also important members of this magmatic assemblage (Amini and Emami, 1992). The major faults in the ILPA area trending NW – SE are as follows: (1) the Koshk-e-Nosrat fault in the south of stations IR4 and IR5; (2) the South Parandak fault in the north of IR5 and in the south of stations IR1 and IR6; (3) a major fault with a length of about 30 km, which crosses the extrusive volcanic zone at the north of stations IR6 and IR1 and terminates at synorogenic conglomerate sediments in the east; and (4) the WNW – ESE-trending Ipak reverse fault consisting of a complex of multiple segments. The 56-km-long eastern section of Ipak fault strikes E – W, N105jE and N70jE, and dips south (close to station IR7), while the 39-km-long western section strikes N110jE and N120jE and dips north (Berberian and Yeats, 2001). The time domain linearized inversion technique of Ammon et al. (1990) is used to invert radial receiver functions to wave velocity structure models. The following steps are taken to prepare velocity – depth functions from the radial receiver functions: – An initial velocity –depth model is defined based on previous geophysical investigations in the area; – A single inversion is performed to indicate the necessary number of iterations; – A smooth inversion is performed to identify a suitable smoothing parameter; – The initial models are perturbed using a pseudoMonte Carlo approach to generate a large number of different starting models; – The radial receiver function is inverted by minimising the difference between the observed receiver function and the synthetics computed for those models, while simultaneously constraining the model smoothness; – Those models that produced a good match to the data are plotted in one figure to classify the main features of the structure; – One of the best-fitted models is selected, and the main features of the structure are derived by group-
ing the thin layers with similar velocities into a single thicker layer; – The adjusted model is used as an initial model and inversion is repeated regarding observed radial receiver function; – Forward modelling is used to produce the synthetic radial receiver function for the adjusted model; and – The results of the forward modelling and inversion are compared to infer a final model. The possible resolution is limited by the frequency content of receiver functions of the data, which in this study are lowpass-filtered at f 0.48 Hz. The corresponding wavelength of the compressional waves is about 12.5 km at 6 km/s average crustal velocity. Based on recent receiver function studies (e.g., Darbyshire et al., 2000), various seismic tests have shown that interfaces separated by more than a quarter wavelength of the seismic waves are resolvable. Therefore, for the shortest compressional wavelength of 12.5 km presented here, the vertical resolution is about 3 km. In practice, our model is parameterised as a stack of horizontal layers to a depth of f 60 km with a layer thickness of 3 km. The inversion estimates the shear wave velocity as a free parameter. Initially, a constant Poisson’s ratio of 0.25 is assumed to calculate the compressional wave velocity. To perturb the initial model, a random perturbation of Poisson’s ratio is added onto layer velocities. The random perturbation is defined to apply broad changes to the initial model. The densities are approximated throughout the model using q = 0.32Vp + 0.77 (Berteussen, 1977), where q and Vp represent density and P-wave velocity, respectively. 3.1. Data analysis of station IR1 To derive the main features of the crustal structure beneath station IR1, 10 teleseismic events have been analysed. The results show strong amplitudes on the radial receiver function compared to the tangential component and the background noise (Fig. 3). In Fig. 3, the radial and tangential receiver functions are arranged based on increasing backazimuths. The radial receiver functions are plotted on the left panel and the tangential receiver functions on the right. The event number, which is written on the left of the each radial receiver function, is related to the event number in
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
121
Fig. 3. The radial and tangential receiver functions of 10 teleseismic earthquakes that were used for station IR1. The radial receiver functions are illustrated on the left side of figure and the tangential receiver functions on the right. First column numbers are related to the event numbers in Table 1. Second and third columns indicate the great circle path of the events and the backazimuths, respectively. The great circle path of events and backazimuths are in terms of degrees. The final column numbers indicate the normalised scaling value of each tangential receiver function to its corresponding radial receiver function.
Table 1. More information of teleseismic events can be found in Table 1. A 180-s time window was used in all processing steps, although only 100 s of receiver functions have been shown in Fig. 3. Three main similar patterns can be seen during the first f 12 s of the radial receiver functions (Fig. 3): (1) a strong direct P-wave arrival (Pp phase) with gradually increasing delay time versus backazimuth; (2) a prominent 5- to 6-s positive P-to-S converted arrival (Ps phase); and (3) a negative arrival between 9 and 11 s. Although the azimuthal coverage is small, a few minor differences in the patterns of events can be observed due to a laterally heterogeneous structure around the receiver. The first three events have the same epicentral distances and show a small positive arrival between 2 and 3 s on the radial receiver functions, while this arrival disappears in the rest of events. The negative arrivals between 9 and 11 s on the radial
receiver functions of the first five events are broader than the second five events. The maximum amplitudes of tangential receiver functions are less than half of their radial components. Relatively strong negative tangential receiver function first arrivals likely suggest a dipping Moho interface. However, deriving a dipping Moho through receiver function analysis does need a full azimuthal coverage. A large secondary phase immediately follows the direct P-wave on the tangential receiver functions, for those events that have backazimuths greater than 100j. Although the arrival time delay of the Pp phase on the radial receiver functions with a backazimuth of >100j indicates the presence of a thick layer of low-velocity sediments, a dipping sediment – basement interface can also be inferred from the large secondary phase on the tangential receiver functions. This result is well consistent with low-velocity heterogeneous sediments and other
122
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
Fig. 4. Inverse modelling result for the radial receiver function of the teleseismic earthquake of December 23, 1978 at station IR1. (a) Observed radial receiver function (thin solid line), synthetic radial receiver function (grey thick line), and observed tangential receiver function (dotted line). (b) Velocity – depth models from inverting the radial receiver function. Details of the process are explained in the text.
Fig. 5. Forward and inverse modelling results for one of the best-fitting models of Fig. 4. (a) One of the best-fitted models (dotted line) selected from Fig. 4 and the adjusted model (solid line). (b) Comparison between the observed radial receiver function (dotted line) and the synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model in comparison with observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of Fig. 4c. (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
123
Fig. 6. Inverse modelling result for stacked radial receiver functions of teleseismic earthquakes of group 1 at station IR1. (a) Stacked radial receiver function (thin solid line) and synthetic radial receiver functions (grey thick line). (b) Velocity – depth models from inverting the stacked radial receiver function.
Fig. 7. Forward and inverse modelling results for one of the best-fitting models of Fig. 6. (a) One of the best-fitting models (dotted line) selected from Fig. 6 and the adjusted model (solid line). (b) Comparison between observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver functions (grey solid line) from inversion of the adjusted model in comparison with observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
124
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
geological evidence of structure (Amini and Emami, 1992; Alavi, 1996) in the south – southeast of station IR1. The result of inverting the radial receiver function of the teleseismic earthquake of December 23, 1978 is illustrated in Fig. 4. Fig. 4a shows the observed radial receiver function (solid thin line), the corresponding synthetic (grey thick line), and the observed tangential receiver function (dotted line). Information about the event is written on the top right corner of the panel. The receiver function is identified by the origin time in the form yyjjjhhmmss, where y = year, j = Julian day, h = hour, m = minute, and s = second. Fig. 4b shows the corresponding velocity – depth models. Considering the lithology around station IR1 and keeping only those models that produced a good match to the data and reasonable velocity – depth models, the main features of the crust easily can be derived. One of the bestfitted models of Fig. 4b was selected and plotted in Fig. 5a (dotted line). The main features of the structure are clarified by grouping the thin layers with similar velocities into a single thicker layer (Fig. 5b, solid
line). This adjusted model is considered as initial model for reinversion of the radial receiver function, for which the results are presented in Fig. 5c and d. The adjusted model also is used to produce a synthetic radial receiver function through forward modelling (Fig. 5b). The results of modelling show that the observed radial receiver functions (dotted line) are consistent with the synthetic (solid line) radial receiver functions (Fig. 5b and c). Comparing the results of the reinversion with forward modelling of this event gives us the primary structural model for IR1: (1) positive velocity gradient from 4 km/s at the surface to 7.2 km/s at the depth of 25 km; (2) negative Ps-converted phase defining the low-velocity zone at the depth of f 30 km; (3) a single thick layer with a constant velocity of 6.4 km/s and a 20-km thickness at the depth of 32 km; and (4) increasing velocity with depth from depth of 52 –60 km in the uppermost mantle. To amplify signal-to-noise ratio and enhance the amplitude of main phases such as Pp, Ps, and PpPms, stacking of radial receiver functions is applied. Owens et al. (1984) suggested that the stacking of events that
Fig. 8. Inverse and forward modelling results for stacked radial receiver function of teleseismic event of group 2 (for details see text) at station IR1. (a) One of the best-fitting models (dotted line) selected from Fig. 6 and the adjusted model (solid line). (b) Comparison between observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver functions (grey solid line) from inversion of the adjusted model in comparison with observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
are separated in distance by 10 –15j may not be degraded amplitudes of multiples even from deep interfaces. However, inverting the stacked radial receiver function provides an average crustal structure. The area of study is surrounded in the west, north, and northeast by the Alborz mountain range, and in the south and southeast by the central Iranian plateau. Therefore, available radial receiver functions are clustered into two groups corresponding to two different tectonic settings. The first group contains a stack of radial receiver functions of five teleseismic earthquakes (i.e., events March 18, 1977; March 19, 1977 (event 7); September 19, 1977; December 23, 1977; and January 14, 1978). These events are coming to station IR1 from a backazimuth of 48 –100j and with a great circle path of 65– 79j. The second group contains a stack of radial receiver functions of four tele-
125
seismic earthquakes (i.e., events March 8, 1977; September 23, 1977; October 16, 1977; and November 18, 1977), which were recorded from backazimuths of 112 –118j and a great circle path of 57 –78j. The results of the inversion of the stacked radial receiver function of group 1 are shown in Fig. 6. A relatively good match between stacked radial receiver function (solid line) and synthetics (grey thick line) up to 30 s (Fig. 6a) indicates reliable velocity – depth structural models (Fig. 6b). To test the stability of this result and to compare with the forward modelling result, one of the best-fitted models of inverse modelling (Fig. 6) is selected to apply reinversion and forward modelling. The results of forward modelling and reinversion are shown in Fig. 7. In the adjusted model, the depth of the uppermost mantle is increased intentionally (Fig. 7a, solid line). By this means, we improved the goodness
Fig. 9. Velocity – depth models from inverse modelling for the seven stations of ILPA network. An arrow marks the approximate position of the Moho boundary.
126
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
of fit up to 35 s (Fig. 7c). The results of inverse modelling of group 1 are similar with the results of the teleseismic earthquake of December 23, 1978. The difference between the two crustal models (Figs. 5d and 7d) is in the P-wave velocities of the uppermost crust. Corresponding results of inverse modelling of the stacked data for group 2 are shown in Fig. 8. A sharp velocity contrast at the depth of 46 km is revealed from an analysis of group 2 (Fig. 8d). This feature, assumed to be the base of the crust, could not be seen from an inverse modelling of the teleseismic earthquake of December 23, 1978 (Fig. 5d) and the stacked data for group 1 (Fig. 7d). The P-wave velocities in the range of 5.8– 6.3 km/s at the depth of 30 km characterise a low-velocity zone in the middle crust that can be seen from inverse modelling results (Figs. 5d, 7d, and 8d). This low-velocity zone may be related to recent volcanism in the area of station IR1. This result can be seen from a receiver function analysis at stations IR4 and IR5, too. That is, the southern part of the area of study shows a common low-velocity zone in the middle crust, which might be related to the recent volcanic activities (Fig. 9). The results of modelling from an analysis of teleseismic earthquakes for the rest of the stations of ILPA network (stations IR2 – IR7) are presented in Appendices A – F.
4. Discussion and conclusions This study permits a thorough investigation of the crust and uppermost mantle structure beneath seven stations of the ILPA network in the southwest of Tehran, Iran. The results of inversion are compared with the forward modelling results at each station. To compare the results and to present a reasonable final velocity structure model for the Tehran region, the results of final inversion in all stations are summarised in Fig. 9. This figure is so configured to place each velocity –depth model resulting from inversion in the position of its own station. The velocity changes with depth at the upper 10 km show a clear dependency on the lithology around each station. The results of modelling show that the P-wave velocity in the middle crust has the same range for stations IR1, IR3, IR5, and IR7. The P-wave velocity in the middle crust also has the same range below stations IR2, IR4, and IR6.
These two groups of stations have P-wave velocities of f 6.3 and f 7.0 km/s in the middle crust, respectively. This wide range of P-wave velocity in the middle crust suggests a complex geological setting in the area of study. The basement of the crust in this area is found to be at the depth of 46 F 2 km. The approximate position of the Moho border in each station is marked with an arrow (Fig. 9). There is no deep refraction seismic profile in this area. Hence, an exact correlation between the velocity structure models and geological evidences remains unclear. Based on laboratory measurements from the northern part of Alborz Mountains (in southwestern Turkmenistan; Volarovich et al., 1980), three main crustal layers can be inferred by their P-wave velocities: (a) sediments, conglomerate, and consolidated sediment (Vp < 4.8 km/s); (b) granitic (Vp between 4.8 and 6.4 km/s); and (c) basaltic (6.4 < Vp < 7.4 km/s). The crust –mantle boundary has a similar velocity contrast at most stations. The upper mantle rocks can be indicated by Vp = 8.0 km/s. Strong tangential receiver function first arrivals may suggest a dipping Moho reflector below the area of study. Our dipping Moho suggestion is consistent with this hypothesis that the juxtaposition of the Alborz system and the central Iranian province is considered to be a result of the early Cenozoic arc – continent collisional processes due to the northeastward subduction of a narrow, elongated back arc basin beneath the western –southeastern part of the Alborz continental block (Alavi, 1996). In this study, we find that a large secondary phase immediately follows the tangential Pp phase for those events coming to station IR1 from the south –southeast. The presence of this large secondary phase is well consistent with the lowvelocity heterogeneous sediments at this area and suggests a sediment – basement-dipping interface. However, due to the small coverage of events, it is difficult define the direction of this sediment –basement-dipping interface. The presence of a low-velocity heterogeneous surface layer can be also inferred from our receiver function inversion results (i.e., Fig. 8). The uppermost mantle shows a positive velocity gradient from 7.5 km/s at a depth of 48 km to 7.9 km/s at a depth of f 60 km. This range of the P-wave velocity in the uppermost mantle indicates a low-velocity mantle beneath the area of study. Our receiver function modelling is the first attempt to derive the main features of the crust and uppermost
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
mantle structure down to 60 km depth below the Tehran region. Therefore, the only chance for discussion is to compare these results with the results of Asudeh (1982) and Sobouti and Arkani-hamed (1996). The results of our receiver function modelling confirm the previous result of Asudeh (1982). However, his results did not reveal a dipping reflector at the basement of the crust. Our estimated crustal thickness is also in good agreement with the results of Sobouti and Arkani-Hamed (1996).
Appendix A . Data analysis of station IR2 Three teleseismic earthquakes (i.e., events February 18, 1977; May 28, 1977; and October 16, 1977) have been analysed for station IR2. The observed radial receiver functions show a weak Ps-converted phase 67 s after the direct Pp phases. The radial receiver functions were stacked and the modelling results are presented in Fig. 10. Comparing velocity – depth models for IR2, the main crustal features can
127
be identified as: (1) the P-wave velocity of the surface layers changes between 5.7 and 6.0 km/s; (2) a weak positive P-wave velocity gradient from 6.2 km/s at a depth of 38 km to f 7.2 km/s at 46 km; and (3) low P-wave velocities in the uppermost mantle from f 7.5 km/s at 48 km to 7.9 km/s at about 60 km.
Appendix B . Data analysis of station IR3 Two teleseismic earthquakes have been analysed to estimate the main features of the structure beneath station IR3. The result of forward modelling and the inversion of radial receiver functions for IR3 show similar velocity –depth structure (Figs. 11 and 12). The observed and synthetic radial receiver functions show a good match for the events January 13, 1978 and June 14, 1978 up to 25 s. The P-wave velocity of the upper 4 km is about 3.5 km/s for one event and 4 km/s for the other. This may be explained by the difference in the backazimuth of
Fig. 10. Inverse and forward modelling results for the stacked radial receiver function of three teleseismic earthquakes at station IR2. (a) One of the best-fitting models (dotted line) selected from the primary inversion of the stacked radial receiver function and the adjusted model (solid line). (b) Comparison between the stacked radial receiver function (dotted line) and synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model in comparison with stacked radial receiver function (dotted line). (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
128
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
Fig. 11. Forward and inverse modelling results for one of the best-fitted model of the earthquake of January 13, 1978 at station IR3. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
Fig. 12. Forward and inverse modelling results for one of the best-fitting models of the earthquake of June 14, 1978 at station IR3. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
two events and lithology around the station. The positive P-wave velocity gradient from 5 to 6.5 km/ s down to a 20-km depth identifies another main crustal feature below station IR3. A P-wave velocity of f 6.7 km/s is the average lower crustal velocity down to 48 km.
Appendix C . Data analysis of station IR4 Only one reliable teleseismic event with three components has been recorded at station IR4. Due to technical problem in the north –south component, other events were not considered suitable for receiver function analysis. The teleseismic event May 12, 1977 occurred at an epicentral distance of 52j from IR4 with a backazimuth of 64.5j. The result of the receiver function analysis of this event is summarised in Fig. 13. Fig. 13a shows one of the best-fitting structural models calculated by inversion and the suggested model for forward modelling and reinversion. Fig. 13b shows the comparison between the
129
results of forward modelling and the observed radial receiver function. The comparison between reinversion results and observed receiver functions is shown in Fig. 13c. These two figures show that the synthetic radial receiver functions are in good agreement with the observed data especially for the main phases such as Pp, Ps, and PpPms. The amplitude of the tangential receiver function is moderate. The corresponding velocity – depth models of the reinversion are presented in the Fig. 13d. From Fig. 13, the main features of the crustal structure beneath station IR4 can be identified as: (1) the surface layers show a strong and positive P-wave velocity gradient with depth; (2) slowly increasing P-wave velocity with depth from 6.5 km/s at a depth of 14 km to 7.0 km/s at 20 km; (3) a thick layer with constant P-wave velocity of 7.2 km/ s down to 32 km; (4) a negative velocity gradient below 2 km and down to 40 km; (5) a P-wave velocity contrast at the depth of about 44 km at the basement of the crust; and (6) a low-velocity mantle with Pwave velocity of 7.2 km/s to 7.9 km/s down to about 60 km.
Fig. 13. Inverse and forward modelling results for the event of May 12, 1977 at station IR4. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
130
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
Appendix D . Data analysis of station IR5 Two events (March 2, 1977 and October 16, 1977) were analysed for station IR5. The same simple crustal structure is deduced from these two different events. Therefore, only the analysis of the second event (October 16, 1977) is illustrated in Fig. 14. These two events arrived at the IR5 station from different backazimuths (i.e., 94j and 112j, respectively) but through the same lithology. Fig. 14a shows one of the best-fitting models and the adjusted model. Fig. 14b illustrates the comparison between forward modelling result and observed radial receiver function. Fig. 14c shows the results of reinversion in comparison with the observed radial and tangential receiver functions. The corresponding velocity –depth model of the reinversion is shown in Fig. 14d. From Fig. 14a and d, three main features can be inferred for the structure beneath station IR5: (1) gradually increasing P-wave velocity with depth from 5.5 km/s in the surface to 6.4 km/s at about 8 km; (2) a thick layer with the
constant P-wave velocity of f 6.3 km/s down to 35 km depth; and (3) gradually increasing P-wave velocity with depth from 6.4 km/s at 36 km to 7.8 km/s at 46 km.
Appendix E . Data analysis of station IR6 Four teleseismic earthquakes (i.e., February 16, 1977; February 18, 1977; February 19, 1977; and March 2, 1977) are analysed for IR6. The information of these events is presented in Table 1. The radial receiver functions of these events have large amplitudes compared to the presignal noise and tangential receiver functions. In spite of different backazimuths, similar main crustal features can be seen from the analysis of each event. To amplify the main phases, radial receiver functions were stacked. The stability of results from inverse modelling of stacked radial receiver function is tested by taking one of the best-fitted models and by applying reinversion and forward modelling. The results of for-
Fig. 14. Inverse and forward modelling results for the teleseismic earthquake of October 16, 1977 at station IR5. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
131
Fig. 15. Inverse and forward modelling results for stacked radial receiver function at station IR6. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from forward modelling. (c) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from inversion. (d) The P-wave velocity – depth models determined from inverse modelling.
Fig. 16. Inverse and forward modelling results for stacked radial receiver function at station IR7. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from forward modelling. (c) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from inversion. (d) The P-wave velocity – depth models determined from inverse modelling.
132
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
ward and inverse modelling of one of the best-fitted models are presented in Fig. 15. In Fig. 15a, one of the best-fitted models (dotted line) and the adjusted model (solid line) are shown. Fig. 15b and c shows the results of the forward modelling (solid line) and inverse modelling (solid grey lines) in comparison with the stacked radial receiver function (dotted line), respectively. The corresponding velocity – depth models resulting from inverse modelling are shown in Fig. 15d. The main features of the crustal structure for IR6 can be indicated as (1) P-wave velocity gradient in the upper crust from 5.4 km/s at the surface to f 7.4 km/s at f 22 km; (2) a constant P-wave velocity of f 6.9 km/s in the middle crust down to f 40 km; and (3) increasing P-wave velocity with depth in the lower crust and crust –mantle boundary zone from f 7 to f 7.9 km/s down to a depth of f 60 km.
Appendix F. Data analysis of station IR7 To determine the main features of the crustal structure beneath station IR7, receiver functions of 10 teleseismic earthquakes are investigated. Then the radial receiver functions of six events (i.e., March 19, 1977, event 7; March 19, 1977, event 8; May 12, 1977; September 19, 1977; September 23, 1977; and 16 October 1977) were stacked. To show the reliability of stacking and receiver function analysis for IR7, one of the best-fitting models from the primary inversion result is selected (Fig. 16a). The results of forward and inverse modelling of the stacked radial receiver function based on the adjusted model are illustrated in Fig. 16b and c. The results of the forward and inverse modelling are similar and show a good match to the three main phases (i.e., Pp, Ps, and PpPms). Therefore, the corresponding velocity – depth models (Fig. 16d) were considered to reliably represent the main features of the crustal structure for IR7. The structural model for IR7 can be thus summarised as: (1) an upper crust with P-wave velocity of 5.6 km/s down to a 13-km depth; (2) a middle crust with P-wave velocity of f 6.3 km/s down to a 30-km depth; (3) a lower crust with P-wave velocity of f 6.7 km/s down to a 47-km depth; and (4) velocity gradient from f 6.8 to 8.05 km/s down to a depth of about 60 km.
Acknowledgements We thank S. Mangino and an anonymous reviewer for their valuable comments and suggestions that improved the paper. We are also grateful to Prof. M. Ghafory-Ashtiani, the president of IIEES, for his kindness and help during revision of the paper. The MSRT of Iran has supported JD.
References Alavi, M., 1996. Tectonostratigraphic synthesis and structural style of the Alborz mountain system in northern Iran. J. Geodyn. 21, 1 – 33. Amini, X., Emami, N.H., 1992. Geological map of the Tehran sheet: scale 1:100,000. Geol. Surv. Iran (Tehran, Iran). Ammon, J.C., 1991. The isolation of receiver function effects from teleseismic P waveforms. Bull. Seismol. Soc. Am. 81, 2504 – 2510. Ammon, J.C., Randal, G.E., Zandt, G., 1990. On the nonuniqueness of receiver function inversion. J. Geophys. Res. 95, 15303 – 15318. Asudeh, I., 1982. Seismic structure of Iran from surface and body wave data. Geophys. J. R. Astron. Soc. 71, 715 – 730. Berberian, M., Yeats, R.S., 2001. Contributions of archaeological data to studies of earthquakes history in the Iranian plateau. J. Struct. Geol. 23, 563 – 584. Berteussen, K.A., 1977. Moho depth determination based on spectral-ratio analysis NORSAR Long-Period P-waves. Phys. Earth Planet. Inter. 15, 13 – 27. Chen, C.Y., Chen, W.P., Molnar, P., 1980. The uppermost mantle P wave velocities beneath Turkey and Iran. Geophys. Res. Lett. 7, 77 – 80. Darbyshire, F.A., Priestley, K.F., White, R.S., Stefansson, R., Gudmundsson, G.B., Jackobsdottir, S., 2000. Crustal structure of central and northern Iceland from analysis of teleseismic receiver functions. Geophys. J. Int. 143, 163 – 184. Dehghani, G.A., Makris, J., 1984. The gravity field and crustal structure of Iran. Neues Jahrb. Geol. Palaeontol. 168, 215 – 229. De Martini, P.M., Hessami, K., Pantosti, D., D’Addezio, G., Alinaghi, H., Ghafory-Ashtiani, M., 1998. A geological contribution to the evaluation of the seismic potential of the Kahrizak fault (Tehran, Iran). Tectonophysics 287, 187 – 199. Giese, P., Makris, J., Akasheh, B., Roewer, P., Letz, H., Mostaanpour, M., 1984. The crustal structure in southern Iran derived from seismic explosion data. Neues Jahrb. Geol. Palaeontol. 168, 230 – 243. Hedayati, A., Brander, J.L., Berberian, M., 1976. Microearthquake survey of Tehran region, Iran. Bull. Seismol. Soc. Am. 66, 1713 – 1725. Langston, C.A., 1977. Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves. Bull. Seismol. Soc. Am. 67, 713 – 724. Langston, C.A., 1979. Structure under Mount Rainier, Washington,
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133 inferred from teleseismic body waves. J. Geophys. Res. 84, 476 – 4749. Mangino, S., Priestley, K., 1998. The crustal structure of the southern Caspian region. Geophys. J. Int. 133, 630 – 648. Owens, T.J., 1987. Crustal structure of the Adirondacks determined from broadband teleseismic waveform modelling. J. Geophys. Res. 92, 6391 – 6401. Owens, T.J., Zandt, G., Taylor, S.R., 1984. Seismic evidence from an ancient rift beneath the Cumberland Plateau, Tennessee: a detailed analysis of broadband teleseismic P waveforms. J. Geophys. Res. 89, 7783 – 7795.
133
Sobouti, F., Arkani-Hamed, J., 1996. Numerical modelling of the deformation of the Iranian plateau. Geophys. J. Int. 126, 805 – 818. Volarovich, M.P., Volynets, L.N., Levshin, A.L., 1980. Composition of the Earth’s crust and the upper mantle from DSS data in central Turkmenia and laboratory measurements for high thermodynamic parameters. Izv. Earth Phys. 16, 441 – 446. Zandt, G., Myers, S.C., Wallace, T.C., 1995. Crustal and mantle structure across the basin and range – Colorado Plateau boundary at 37jN latitude and implications for Cenozoic extensional mechanism. J. Geophys. R. 100, 10529 – 10548.
Crust and uppermost mantle structure of Tehran region from analysis of teleseismic P-waveform receiver functions Javan Doloei a,b,*, Roland Roberts a a b
Department of Earth Sciences and Geophysics, Uppsala University, Villava¨gen 16, Uppsala S-752 36, Sweden International Institute of Earthquake Engineering and Seismology (IIEES), P.O. Box 19395/3913, Tehran, Iran Received 7 August 2001; accepted 3 February 2003
Abstract Teleseismic P-waveform receiver function analysis has been used to estimate the main features of the crust and uppermost mantle structure of the Tehran region. In this study, 20 teleseismic earthquakes that were recorded at the seven stations of Iranian Long Period Array (ILPA) were used to calculate receiver functions. The radial receiver functions are analysed through inverse modelling and forward modelling. The results of modelling show a good match between observed radial receiver functions and synthetic radial receiver functions up to 30 s. Based on this study, the crustal structure of this area is characterised by three main layers: The upper crust has a P-wave velocity between 4 and 5.8 km/s and a 14-km thickness. The results representing the uppermost crustal structure indicate a strong dependency on the lithology around each station. The middle crust has a positive P-wave velocity gradient from 6 to 6.4 km/s down to f 30 km depth. A P-wave velocity gradient from 6.4 to 7.5 km/s characterises the lower crust with Moho at 46 F 2 km at this area. A low-velocity mantle, with velocity gradually increasing with depth down to f 60 km depth, has been revealed beneath the Tehran region. D 2003 Elsevier Science B.V. All rights reserved. Keywords: Crust; Iran; Receiver function; Uppermost mantle structure; Seismology; Tehran
1. Introduction The major crustal features and the seismic velocities of the uppermost mantle are important constraints in understanding the tectonic evolution of an area. Pwaveform teleseismic receiver function analysis is becoming a common method to estimate this information. In this study, the main features of the crust and
* Corresponding author. International Institute of Earthquake Engineering and Seismology (IIEES), P.O. Box 19395/3913, Tehran, Iran. E-mail address: [email protected] (J. Doloei).
uppermost mantle structure beneath the Iranian Long Period Array (ILPA) are investigated using the receiver function technique. The ILPA network consists of seven stations, which are located southwest of the metropolis of Tehran. This area is situated at the southern border of the central Alborz mountains, which are the branches of the Alpine – Himalayan orogeny in Iran (Fig. 1, top). Iran lies within the continental collision zone between Africa and Eurasia. Active faulting, recent volcanism, and high surface elevations are important characteristics of this area. Geomorphic and trenching investigations (De Martini et al., 1998) along the Kahrizak fault, which is
0040-1951/03/$ - see front matter D 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0040-1951(03)00049-0
116
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
117
Table 1 Information on the 20 teleseismic earthquakes recorded at the ILPA network Event number
Date
Origin time (hhmmss)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1977/2/16 1977/2/18 1977/2/19 1977/3/2 1977/3/8 1977/3/18 1977/3/19 1977/3/19 1977/5/12 1977/5/28 1977/9/19 1977/9/23 1977/10/16 1977/11/18 1977/12/23 1978/1/13 1978/1/14 1978/6/14 1978/7/23 1978/12/23
104020.90 205129.80 223404.10 095323.20 231728 214352.40 105625.10 193508 111753.10 055148.10 051309.20 055755.60 210917.70 101741 210207.50 200303.80 032439 113420 144236.90 112312
Latitude (j) 0.50 33.07 53.57 6.77 0.45 16.77 44.20 16.81 39.27 1.73 1.98 11.21 9.73 4.35 39.13 44.69 34.81 38.32 22.28 23.25
Longitude (j)
Depth (km)
Magnitude
GCARC (j)
Backazimuth (j)
125.98 140.82 170.03 123.74 100.02 122.33 148.20 122.35 117.71 120.52 126.63 118.22 117.12 102.02 143.16 149.70 139.26 142.36 121.51 122.07
33 42 33 52 22 37 70 39 22 54 33 33 33 33 19 52 14 40 17 33
6.1 6.0 6.7 6.1 6.0 7.0 6.0 5.8 5.8 6.0 6.2 6.0 6.0 5.9 6.0 6.0 6.7 6.3 7.4 7.2
78 72 77 72 57 65 71 65 52 74 79 78 77 62 70 71 70 70 62 62
98 62 32 94 116 86 48 86 64 103 100 112 112 118 55 48 61 56 81 80
The GCARC and backazimuths are related to station IR1.
located south of Tehran and west of the ILPA network, provided some new constraints on the activity of this region. Previous studies of seismic structure of Iran give us a general view of the crust and uppermost mantle across the Iranian plateau. The earliest study of Hedayati et al. (1976) utilised a microearthquake survey to show a 31-km crustal thickness below the Tehran region. A general study of seismic P-wave and S-wave velocities by using arrival time data of regional earthquakes (Chen et al., 1980) provided seismic structure information beneath Iran with average Pn velocity of 8.0 F 0.1 km/s. Surface wave analysis of a few events (Asudeh, 1982) suggested a one-layer model crust with 45 km thickness beneath the eastern part of Alborz mountain range. Asudeh (1982) showed that the crust and upper mantle in the eastern part of Iran are characterised by low body wave velocities. The gravity field and crustal structure of Iran, in the study by Dehghani and Makris (1984),
showed that the Bouguer gravity along the Alborz mountain range varies between 100 and 120 mgal and that the crustal thickness is less than 35 km. That is, these mountains have no root. The first and only deep seismic soundings in Iran were conducted in 1978 between the Lut Block in the east of the country and the Zagros mountain range in the west. The results of this study identified a 40-km crustal thickness for the central Iranian plateau (Giese et al., 1984). Numerical modelling of the deformation of the Iranian plateau (Sobouti and Arkani-Hamed, 1996) identified that crustal thickening and shear deformation tend to localise near the southwest of the country. Their results specified a 45-km crustal thickness along the Alborz mountain range. Based on receiver function analysis of teleseismic data, Mangino and Priestley (1998) provided more information on the crustal structure of the southern Caspian region in the north of the Alborz mountain range. To give an
Fig. 1. (Top) Map of major structural divisions of Iran. Rectangle = research area (bottom). (Bottom) Location map of stations (triangles) for ILPA network in the Tehran area. Q = unfolded Quaternary cover alluvium; Qpl = Paleocene – Eocene, deposited on pre-Paleocene folded substratum, folded by middle and late Alpine events, conglomerate, limestone, and volcanic rocks; A3 OM = sediments deposited on pre-Oligocene substratum, folded during late Alpine event; Ev = extrusive volcanic rocks, mainly andesites and pyroclastics, acidic and intermediate basic.
118
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
improved picture of the Iranian part of the tectonically active Alpine – Himalayan orogenic belt, a receiver function study of available teleseismic ILPA data is necessary and useful. The purpose of this paper is to determine the crustal and uppermost mantle seismic velocity structures of the Tehran area based on a receiver function analysis of 20 teleseismic earthquakes, which occurred in the range of 52– 79j from centre of the ILPA network with magnitude greater than 5.8. Table 1 gives information on the teleseismic earthquakes used. The information of teleseismic earthquakes of Table 1 is extracted from PDE catalogue of US Geological Survey.
2. Receiver function analysis The time domain receiver function technique has been described in detail (Langston, 1977, 1979; Owens, 1987; Owens et al., 1984; Ammon et al., 1990) for use with single-station three-component data. Deconvolving the vertical seismogram from the radial component gives us the radial receiver function. A similar deconvolution of the vertical component from the transverse component provides the tangential receiver function, which is useful to test the data regarding the assumptions of a laterally homogeneous, isotropic Earth structure. It is well known that a teleseismic P-wave incident to a velocity discontinuity (e.g., the Moho) will generally produce P-waves and S-waves [both transmitted (Pp, Ps) and Moho reverberation phases such as PpPms and PpSms + PsPms]. Because the S-wave will travel more slowly than the Pwave, the time difference in the arrival of the two phases Pp and Ps is a direct measure of the depth of the original interface. The time interval between the Ps converted phase and the Moho PpPms reverberation phase gives the ratio of P-wave and S-wave velocities (Zandt et al., 1995). The seismic crustal and uppermost mantle structure of the Tehran area is investigated using teleseismic Pwaveforms recorded at the ILPA network for 2 years (1977 – 1978). The continuous waveform database was obtained from the Data Management Centre (DMC) of the Incorporated Research Institution for Seismology (IRIS). The ILPA network consists of seven stations IR1 –IR7. Six stations (IR2 –IR7) are placed around a
circle with a radius of about 30 km, and one station (i.e., IR1) is located in the centre (Fig. 1, bottom). This configuration allows investigation of the nonuniqueness of the deduced velocity models through comparing and combining the results of receiver function inversions from all stations. Due to the low sampling rate of Long Period data (1 Hz), a time window with a length of 180 s of Pwaveforms is used to calculate true amplitude receiver functions (Ammon, 1991). To separate actual teleseismic P-waveforms from noise and tapering effects, all time windows of three-component seismograms are selected 60 s before and 120 s after P-wave onset. The maximum frequency content of this type of data is 0.5 Hz, which is sufficient for resolving the main features of the crust and uppermost mantle structure around the three-component single station. The mean and a trend are removed from the selected time window of threecomponent observations and a 20% cosine taper is applied to the signal. The gains of horizontal and vertical components of borehole sensor are the same at each station of ILPA network. Therefore, it was not necessary to correct instrument gains. The true amplitude receiver function technique of Ammon (1991) is used to calculate the radial and tangential P-waveform receiver functions. To stabilise deconvolution, a water level parameter of c = 0.000001 is used. This value was selected after various trials in the range c = 1.0 10 8 to 1.0 100. The value chosen produces acceptable noise levels in most of our receiver functions. The averaging function is defined by a deconvolution of the vertical component from itself at each water level c to find out the narrowest possible Gaussian curve. Due to the lowlevel noise in our observations, a small value of water level gives minimum distortion effects. Applying a Gaussian filter (Ammon, 1991) controls the frequency content of the receiver functions. A Gaussian filter width parameter a = 1.015 is used, which lowpasses the waveforms at f 0.48 Hz. As an example, Fig. 2 shows a three-component seismogram of a teleseismic earthquake (December 23, 1978), with a blow-up of the P-waveform and corresponding receiver functions. This event was recorded at station IR1 with a surface wave magnitude of 7.2 and an epicentral distance of 62j. Although the amplitude of PpPms phases is weak due to surface topography, fortunately, a few teleseismic earthquakes in this study
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
119
Fig. 2. (a) Three-component seismogram of the earthquake of December 23, 1978 that was recorded at station IR1 in BAZ = 80j. Origin time is fixed at zero time so the time of different phases shows the travel time of each phase from source to receiver. (b) Windowing P-waveform; a lag window of the P-waveform with a length of 180 s is selected by taking 60 s of the signal before and 120 s after P-onset. (c) The corresponding radial and tangential receiver functions. The Pp phase is plotted with a 10-s delay.
clearly show PpPms phases on radial receiver functions at the other stations.
3. Velocity structure of Tehran area from teleseismic receiver functions The ILPA network is situated in a complex tectonic setting in northcentral Iran, surrounded in the west, north, and northeast by the active Alpine – Himalayan orogeny (i.e., the Alborz mountain range; Fig. 1, bottom). The ILPA network is so configured so that two stations IR1 and IR6 are located on an extrusive volcanic zone trending NW – SE. In the north, stations IR2 and IR7 are also located on the volcanic rocks on the border of the Qazvin depression. Station IR3 is located in the northeast of the ILPA network, surrounded in the west and north by the recent volcanic rocks, and in the east and south by the unfolded Quaternary cover alluviums and synorogenic con-
glomerate sediments. Stations IR4 and IR5 are placed on the Eocene – Oligocene Karaj formation in the south of the ILPA network (Fig. 1, bottom). Generally, the Cenozoic sediments on the southern side of the Alborz mountain range constitute a heterogeneous discontinuous package of polymicts siliciclastics with subordinate carbonate layers in the lower part with clear synorogenic characteristics. The lower part of the succession consists of lensoidal (up to 1500 m thick) elongated masses of thick-bedded calcarenites (Alavi, 1996). The Neogene– Quaternary sediments (including the Hezar-Darreh and Kahrizak Formations) show every aspect of sedimentation and include coarse to fine, northward-coarsening, polymict fanglomerates, and breccias with reworked sediments of the Neogene age that are mainly derived from the Alborz magmatic assemblage. The Alborz magmatic assemblage (including the Karaj Formation of Eocene age) consists of both submarine and subaerial, porphyritic and nonporphyritic, massive lava flows of andesitic and
120
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
andesitic – basaltic compositions; extensive well bedded; and calcareous and noncalcareous pyroclastics such as tuffs and agglomerates (andesitic and trachyandesitic in composition). Volcanic breccias (widely distributed volcaniclastic sandstones) and red, coarse, polymict conglomerates (which are intruded by diorites, monzodiorites, and diabase dikes and sills) are also important members of this magmatic assemblage (Amini and Emami, 1992). The major faults in the ILPA area trending NW – SE are as follows: (1) the Koshk-e-Nosrat fault in the south of stations IR4 and IR5; (2) the South Parandak fault in the north of IR5 and in the south of stations IR1 and IR6; (3) a major fault with a length of about 30 km, which crosses the extrusive volcanic zone at the north of stations IR6 and IR1 and terminates at synorogenic conglomerate sediments in the east; and (4) the WNW – ESE-trending Ipak reverse fault consisting of a complex of multiple segments. The 56-km-long eastern section of Ipak fault strikes E – W, N105jE and N70jE, and dips south (close to station IR7), while the 39-km-long western section strikes N110jE and N120jE and dips north (Berberian and Yeats, 2001). The time domain linearized inversion technique of Ammon et al. (1990) is used to invert radial receiver functions to wave velocity structure models. The following steps are taken to prepare velocity – depth functions from the radial receiver functions: – An initial velocity –depth model is defined based on previous geophysical investigations in the area; – A single inversion is performed to indicate the necessary number of iterations; – A smooth inversion is performed to identify a suitable smoothing parameter; – The initial models are perturbed using a pseudoMonte Carlo approach to generate a large number of different starting models; – The radial receiver function is inverted by minimising the difference between the observed receiver function and the synthetics computed for those models, while simultaneously constraining the model smoothness; – Those models that produced a good match to the data are plotted in one figure to classify the main features of the structure; – One of the best-fitted models is selected, and the main features of the structure are derived by group-
ing the thin layers with similar velocities into a single thicker layer; – The adjusted model is used as an initial model and inversion is repeated regarding observed radial receiver function; – Forward modelling is used to produce the synthetic radial receiver function for the adjusted model; and – The results of the forward modelling and inversion are compared to infer a final model. The possible resolution is limited by the frequency content of receiver functions of the data, which in this study are lowpass-filtered at f 0.48 Hz. The corresponding wavelength of the compressional waves is about 12.5 km at 6 km/s average crustal velocity. Based on recent receiver function studies (e.g., Darbyshire et al., 2000), various seismic tests have shown that interfaces separated by more than a quarter wavelength of the seismic waves are resolvable. Therefore, for the shortest compressional wavelength of 12.5 km presented here, the vertical resolution is about 3 km. In practice, our model is parameterised as a stack of horizontal layers to a depth of f 60 km with a layer thickness of 3 km. The inversion estimates the shear wave velocity as a free parameter. Initially, a constant Poisson’s ratio of 0.25 is assumed to calculate the compressional wave velocity. To perturb the initial model, a random perturbation of Poisson’s ratio is added onto layer velocities. The random perturbation is defined to apply broad changes to the initial model. The densities are approximated throughout the model using q = 0.32Vp + 0.77 (Berteussen, 1977), where q and Vp represent density and P-wave velocity, respectively. 3.1. Data analysis of station IR1 To derive the main features of the crustal structure beneath station IR1, 10 teleseismic events have been analysed. The results show strong amplitudes on the radial receiver function compared to the tangential component and the background noise (Fig. 3). In Fig. 3, the radial and tangential receiver functions are arranged based on increasing backazimuths. The radial receiver functions are plotted on the left panel and the tangential receiver functions on the right. The event number, which is written on the left of the each radial receiver function, is related to the event number in
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
121
Fig. 3. The radial and tangential receiver functions of 10 teleseismic earthquakes that were used for station IR1. The radial receiver functions are illustrated on the left side of figure and the tangential receiver functions on the right. First column numbers are related to the event numbers in Table 1. Second and third columns indicate the great circle path of the events and the backazimuths, respectively. The great circle path of events and backazimuths are in terms of degrees. The final column numbers indicate the normalised scaling value of each tangential receiver function to its corresponding radial receiver function.
Table 1. More information of teleseismic events can be found in Table 1. A 180-s time window was used in all processing steps, although only 100 s of receiver functions have been shown in Fig. 3. Three main similar patterns can be seen during the first f 12 s of the radial receiver functions (Fig. 3): (1) a strong direct P-wave arrival (Pp phase) with gradually increasing delay time versus backazimuth; (2) a prominent 5- to 6-s positive P-to-S converted arrival (Ps phase); and (3) a negative arrival between 9 and 11 s. Although the azimuthal coverage is small, a few minor differences in the patterns of events can be observed due to a laterally heterogeneous structure around the receiver. The first three events have the same epicentral distances and show a small positive arrival between 2 and 3 s on the radial receiver functions, while this arrival disappears in the rest of events. The negative arrivals between 9 and 11 s on the radial
receiver functions of the first five events are broader than the second five events. The maximum amplitudes of tangential receiver functions are less than half of their radial components. Relatively strong negative tangential receiver function first arrivals likely suggest a dipping Moho interface. However, deriving a dipping Moho through receiver function analysis does need a full azimuthal coverage. A large secondary phase immediately follows the direct P-wave on the tangential receiver functions, for those events that have backazimuths greater than 100j. Although the arrival time delay of the Pp phase on the radial receiver functions with a backazimuth of >100j indicates the presence of a thick layer of low-velocity sediments, a dipping sediment – basement interface can also be inferred from the large secondary phase on the tangential receiver functions. This result is well consistent with low-velocity heterogeneous sediments and other
122
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
Fig. 4. Inverse modelling result for the radial receiver function of the teleseismic earthquake of December 23, 1978 at station IR1. (a) Observed radial receiver function (thin solid line), synthetic radial receiver function (grey thick line), and observed tangential receiver function (dotted line). (b) Velocity – depth models from inverting the radial receiver function. Details of the process are explained in the text.
Fig. 5. Forward and inverse modelling results for one of the best-fitting models of Fig. 4. (a) One of the best-fitted models (dotted line) selected from Fig. 4 and the adjusted model (solid line). (b) Comparison between the observed radial receiver function (dotted line) and the synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model in comparison with observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of Fig. 4c. (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
123
Fig. 6. Inverse modelling result for stacked radial receiver functions of teleseismic earthquakes of group 1 at station IR1. (a) Stacked radial receiver function (thin solid line) and synthetic radial receiver functions (grey thick line). (b) Velocity – depth models from inverting the stacked radial receiver function.
Fig. 7. Forward and inverse modelling results for one of the best-fitting models of Fig. 6. (a) One of the best-fitting models (dotted line) selected from Fig. 6 and the adjusted model (solid line). (b) Comparison between observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver functions (grey solid line) from inversion of the adjusted model in comparison with observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
124
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
geological evidence of structure (Amini and Emami, 1992; Alavi, 1996) in the south – southeast of station IR1. The result of inverting the radial receiver function of the teleseismic earthquake of December 23, 1978 is illustrated in Fig. 4. Fig. 4a shows the observed radial receiver function (solid thin line), the corresponding synthetic (grey thick line), and the observed tangential receiver function (dotted line). Information about the event is written on the top right corner of the panel. The receiver function is identified by the origin time in the form yyjjjhhmmss, where y = year, j = Julian day, h = hour, m = minute, and s = second. Fig. 4b shows the corresponding velocity – depth models. Considering the lithology around station IR1 and keeping only those models that produced a good match to the data and reasonable velocity – depth models, the main features of the crust easily can be derived. One of the bestfitted models of Fig. 4b was selected and plotted in Fig. 5a (dotted line). The main features of the structure are clarified by grouping the thin layers with similar velocities into a single thicker layer (Fig. 5b, solid
line). This adjusted model is considered as initial model for reinversion of the radial receiver function, for which the results are presented in Fig. 5c and d. The adjusted model also is used to produce a synthetic radial receiver function through forward modelling (Fig. 5b). The results of modelling show that the observed radial receiver functions (dotted line) are consistent with the synthetic (solid line) radial receiver functions (Fig. 5b and c). Comparing the results of the reinversion with forward modelling of this event gives us the primary structural model for IR1: (1) positive velocity gradient from 4 km/s at the surface to 7.2 km/s at the depth of 25 km; (2) negative Ps-converted phase defining the low-velocity zone at the depth of f 30 km; (3) a single thick layer with a constant velocity of 6.4 km/s and a 20-km thickness at the depth of 32 km; and (4) increasing velocity with depth from depth of 52 –60 km in the uppermost mantle. To amplify signal-to-noise ratio and enhance the amplitude of main phases such as Pp, Ps, and PpPms, stacking of radial receiver functions is applied. Owens et al. (1984) suggested that the stacking of events that
Fig. 8. Inverse and forward modelling results for stacked radial receiver function of teleseismic event of group 2 (for details see text) at station IR1. (a) One of the best-fitting models (dotted line) selected from Fig. 6 and the adjusted model (solid line). (b) Comparison between observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver functions (grey solid line) from inversion of the adjusted model in comparison with observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
are separated in distance by 10 –15j may not be degraded amplitudes of multiples even from deep interfaces. However, inverting the stacked radial receiver function provides an average crustal structure. The area of study is surrounded in the west, north, and northeast by the Alborz mountain range, and in the south and southeast by the central Iranian plateau. Therefore, available radial receiver functions are clustered into two groups corresponding to two different tectonic settings. The first group contains a stack of radial receiver functions of five teleseismic earthquakes (i.e., events March 18, 1977; March 19, 1977 (event 7); September 19, 1977; December 23, 1977; and January 14, 1978). These events are coming to station IR1 from a backazimuth of 48 –100j and with a great circle path of 65– 79j. The second group contains a stack of radial receiver functions of four tele-
125
seismic earthquakes (i.e., events March 8, 1977; September 23, 1977; October 16, 1977; and November 18, 1977), which were recorded from backazimuths of 112 –118j and a great circle path of 57 –78j. The results of the inversion of the stacked radial receiver function of group 1 are shown in Fig. 6. A relatively good match between stacked radial receiver function (solid line) and synthetics (grey thick line) up to 30 s (Fig. 6a) indicates reliable velocity – depth structural models (Fig. 6b). To test the stability of this result and to compare with the forward modelling result, one of the best-fitted models of inverse modelling (Fig. 6) is selected to apply reinversion and forward modelling. The results of forward modelling and reinversion are shown in Fig. 7. In the adjusted model, the depth of the uppermost mantle is increased intentionally (Fig. 7a, solid line). By this means, we improved the goodness
Fig. 9. Velocity – depth models from inverse modelling for the seven stations of ILPA network. An arrow marks the approximate position of the Moho boundary.
126
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
of fit up to 35 s (Fig. 7c). The results of inverse modelling of group 1 are similar with the results of the teleseismic earthquake of December 23, 1978. The difference between the two crustal models (Figs. 5d and 7d) is in the P-wave velocities of the uppermost crust. Corresponding results of inverse modelling of the stacked data for group 2 are shown in Fig. 8. A sharp velocity contrast at the depth of 46 km is revealed from an analysis of group 2 (Fig. 8d). This feature, assumed to be the base of the crust, could not be seen from an inverse modelling of the teleseismic earthquake of December 23, 1978 (Fig. 5d) and the stacked data for group 1 (Fig. 7d). The P-wave velocities in the range of 5.8– 6.3 km/s at the depth of 30 km characterise a low-velocity zone in the middle crust that can be seen from inverse modelling results (Figs. 5d, 7d, and 8d). This low-velocity zone may be related to recent volcanism in the area of station IR1. This result can be seen from a receiver function analysis at stations IR4 and IR5, too. That is, the southern part of the area of study shows a common low-velocity zone in the middle crust, which might be related to the recent volcanic activities (Fig. 9). The results of modelling from an analysis of teleseismic earthquakes for the rest of the stations of ILPA network (stations IR2 – IR7) are presented in Appendices A – F.
4. Discussion and conclusions This study permits a thorough investigation of the crust and uppermost mantle structure beneath seven stations of the ILPA network in the southwest of Tehran, Iran. The results of inversion are compared with the forward modelling results at each station. To compare the results and to present a reasonable final velocity structure model for the Tehran region, the results of final inversion in all stations are summarised in Fig. 9. This figure is so configured to place each velocity –depth model resulting from inversion in the position of its own station. The velocity changes with depth at the upper 10 km show a clear dependency on the lithology around each station. The results of modelling show that the P-wave velocity in the middle crust has the same range for stations IR1, IR3, IR5, and IR7. The P-wave velocity in the middle crust also has the same range below stations IR2, IR4, and IR6.
These two groups of stations have P-wave velocities of f 6.3 and f 7.0 km/s in the middle crust, respectively. This wide range of P-wave velocity in the middle crust suggests a complex geological setting in the area of study. The basement of the crust in this area is found to be at the depth of 46 F 2 km. The approximate position of the Moho border in each station is marked with an arrow (Fig. 9). There is no deep refraction seismic profile in this area. Hence, an exact correlation between the velocity structure models and geological evidences remains unclear. Based on laboratory measurements from the northern part of Alborz Mountains (in southwestern Turkmenistan; Volarovich et al., 1980), three main crustal layers can be inferred by their P-wave velocities: (a) sediments, conglomerate, and consolidated sediment (Vp < 4.8 km/s); (b) granitic (Vp between 4.8 and 6.4 km/s); and (c) basaltic (6.4 < Vp < 7.4 km/s). The crust –mantle boundary has a similar velocity contrast at most stations. The upper mantle rocks can be indicated by Vp = 8.0 km/s. Strong tangential receiver function first arrivals may suggest a dipping Moho reflector below the area of study. Our dipping Moho suggestion is consistent with this hypothesis that the juxtaposition of the Alborz system and the central Iranian province is considered to be a result of the early Cenozoic arc – continent collisional processes due to the northeastward subduction of a narrow, elongated back arc basin beneath the western –southeastern part of the Alborz continental block (Alavi, 1996). In this study, we find that a large secondary phase immediately follows the tangential Pp phase for those events coming to station IR1 from the south –southeast. The presence of this large secondary phase is well consistent with the lowvelocity heterogeneous sediments at this area and suggests a sediment – basement-dipping interface. However, due to the small coverage of events, it is difficult define the direction of this sediment –basement-dipping interface. The presence of a low-velocity heterogeneous surface layer can be also inferred from our receiver function inversion results (i.e., Fig. 8). The uppermost mantle shows a positive velocity gradient from 7.5 km/s at a depth of 48 km to 7.9 km/s at a depth of f 60 km. This range of the P-wave velocity in the uppermost mantle indicates a low-velocity mantle beneath the area of study. Our receiver function modelling is the first attempt to derive the main features of the crust and uppermost
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
mantle structure down to 60 km depth below the Tehran region. Therefore, the only chance for discussion is to compare these results with the results of Asudeh (1982) and Sobouti and Arkani-hamed (1996). The results of our receiver function modelling confirm the previous result of Asudeh (1982). However, his results did not reveal a dipping reflector at the basement of the crust. Our estimated crustal thickness is also in good agreement with the results of Sobouti and Arkani-Hamed (1996).
Appendix A . Data analysis of station IR2 Three teleseismic earthquakes (i.e., events February 18, 1977; May 28, 1977; and October 16, 1977) have been analysed for station IR2. The observed radial receiver functions show a weak Ps-converted phase 67 s after the direct Pp phases. The radial receiver functions were stacked and the modelling results are presented in Fig. 10. Comparing velocity – depth models for IR2, the main crustal features can
127
be identified as: (1) the P-wave velocity of the surface layers changes between 5.7 and 6.0 km/s; (2) a weak positive P-wave velocity gradient from 6.2 km/s at a depth of 38 km to f 7.2 km/s at 46 km; and (3) low P-wave velocities in the uppermost mantle from f 7.5 km/s at 48 km to 7.9 km/s at about 60 km.
Appendix B . Data analysis of station IR3 Two teleseismic earthquakes have been analysed to estimate the main features of the structure beneath station IR3. The result of forward modelling and the inversion of radial receiver functions for IR3 show similar velocity –depth structure (Figs. 11 and 12). The observed and synthetic radial receiver functions show a good match for the events January 13, 1978 and June 14, 1978 up to 25 s. The P-wave velocity of the upper 4 km is about 3.5 km/s for one event and 4 km/s for the other. This may be explained by the difference in the backazimuth of
Fig. 10. Inverse and forward modelling results for the stacked radial receiver function of three teleseismic earthquakes at station IR2. (a) One of the best-fitting models (dotted line) selected from the primary inversion of the stacked radial receiver function and the adjusted model (solid line). (b) Comparison between the stacked radial receiver function (dotted line) and synthetic radial receiver function (solid line) through forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model in comparison with stacked radial receiver function (dotted line). (d) The P-wave velocity – depth models calculated from inversion of the adjusted model.
128
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
Fig. 11. Forward and inverse modelling results for one of the best-fitted model of the earthquake of January 13, 1978 at station IR3. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
Fig. 12. Forward and inverse modelling results for one of the best-fitting models of the earthquake of June 14, 1978 at station IR3. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
two events and lithology around the station. The positive P-wave velocity gradient from 5 to 6.5 km/ s down to a 20-km depth identifies another main crustal feature below station IR3. A P-wave velocity of f 6.7 km/s is the average lower crustal velocity down to 48 km.
Appendix C . Data analysis of station IR4 Only one reliable teleseismic event with three components has been recorded at station IR4. Due to technical problem in the north –south component, other events were not considered suitable for receiver function analysis. The teleseismic event May 12, 1977 occurred at an epicentral distance of 52j from IR4 with a backazimuth of 64.5j. The result of the receiver function analysis of this event is summarised in Fig. 13. Fig. 13a shows one of the best-fitting structural models calculated by inversion and the suggested model for forward modelling and reinversion. Fig. 13b shows the comparison between the
129
results of forward modelling and the observed radial receiver function. The comparison between reinversion results and observed receiver functions is shown in Fig. 13c. These two figures show that the synthetic radial receiver functions are in good agreement with the observed data especially for the main phases such as Pp, Ps, and PpPms. The amplitude of the tangential receiver function is moderate. The corresponding velocity – depth models of the reinversion are presented in the Fig. 13d. From Fig. 13, the main features of the crustal structure beneath station IR4 can be identified as: (1) the surface layers show a strong and positive P-wave velocity gradient with depth; (2) slowly increasing P-wave velocity with depth from 6.5 km/s at a depth of 14 km to 7.0 km/s at 20 km; (3) a thick layer with constant P-wave velocity of 7.2 km/ s down to 32 km; (4) a negative velocity gradient below 2 km and down to 40 km; (5) a P-wave velocity contrast at the depth of about 44 km at the basement of the crust; and (6) a low-velocity mantle with Pwave velocity of 7.2 km/s to 7.9 km/s down to about 60 km.
Fig. 13. Inverse and forward modelling results for the event of May 12, 1977 at station IR4. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
130
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
Appendix D . Data analysis of station IR5 Two events (March 2, 1977 and October 16, 1977) were analysed for station IR5. The same simple crustal structure is deduced from these two different events. Therefore, only the analysis of the second event (October 16, 1977) is illustrated in Fig. 14. These two events arrived at the IR5 station from different backazimuths (i.e., 94j and 112j, respectively) but through the same lithology. Fig. 14a shows one of the best-fitting models and the adjusted model. Fig. 14b illustrates the comparison between forward modelling result and observed radial receiver function. Fig. 14c shows the results of reinversion in comparison with the observed radial and tangential receiver functions. The corresponding velocity –depth model of the reinversion is shown in Fig. 14d. From Fig. 14a and d, three main features can be inferred for the structure beneath station IR5: (1) gradually increasing P-wave velocity with depth from 5.5 km/s in the surface to 6.4 km/s at about 8 km; (2) a thick layer with the
constant P-wave velocity of f 6.3 km/s down to 35 km depth; and (3) gradually increasing P-wave velocity with depth from 6.4 km/s at 36 km to 7.8 km/s at 46 km.
Appendix E . Data analysis of station IR6 Four teleseismic earthquakes (i.e., February 16, 1977; February 18, 1977; February 19, 1977; and March 2, 1977) are analysed for IR6. The information of these events is presented in Table 1. The radial receiver functions of these events have large amplitudes compared to the presignal noise and tangential receiver functions. In spite of different backazimuths, similar main crustal features can be seen from the analysis of each event. To amplify the main phases, radial receiver functions were stacked. The stability of results from inverse modelling of stacked radial receiver function is tested by taking one of the best-fitted models and by applying reinversion and forward modelling. The results of for-
Fig. 14. Inverse and forward modelling results for the teleseismic earthquake of October 16, 1977 at station IR5. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) The observed radial receiver function (dotted line) and synthetic radial receiver function (solid line) resulting from forward modelling. (c) The synthetic radial receiver function (grey solid line) from inversion of the adjusted model and the observed radial receiver function (dotted line). The observed tangential receiver function (dashed line) is plotted in the lower part of (c). (d) The P-wave velocity – depth models calculated from the inversion of the adjusted model.
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
131
Fig. 15. Inverse and forward modelling results for stacked radial receiver function at station IR6. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from forward modelling. (c) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from inversion. (d) The P-wave velocity – depth models determined from inverse modelling.
Fig. 16. Inverse and forward modelling results for stacked radial receiver function at station IR7. (a) One of the best-fitting models (dotted line) and the adjusted model (solid line). (b) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from forward modelling. (c) Stacked radial receiver function (dotted line) and the synthetic radial receiver function (solid line) resulting from inversion. (d) The P-wave velocity – depth models determined from inverse modelling.
132
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133
ward and inverse modelling of one of the best-fitted models are presented in Fig. 15. In Fig. 15a, one of the best-fitted models (dotted line) and the adjusted model (solid line) are shown. Fig. 15b and c shows the results of the forward modelling (solid line) and inverse modelling (solid grey lines) in comparison with the stacked radial receiver function (dotted line), respectively. The corresponding velocity – depth models resulting from inverse modelling are shown in Fig. 15d. The main features of the crustal structure for IR6 can be indicated as (1) P-wave velocity gradient in the upper crust from 5.4 km/s at the surface to f 7.4 km/s at f 22 km; (2) a constant P-wave velocity of f 6.9 km/s in the middle crust down to f 40 km; and (3) increasing P-wave velocity with depth in the lower crust and crust –mantle boundary zone from f 7 to f 7.9 km/s down to a depth of f 60 km.
Appendix F. Data analysis of station IR7 To determine the main features of the crustal structure beneath station IR7, receiver functions of 10 teleseismic earthquakes are investigated. Then the radial receiver functions of six events (i.e., March 19, 1977, event 7; March 19, 1977, event 8; May 12, 1977; September 19, 1977; September 23, 1977; and 16 October 1977) were stacked. To show the reliability of stacking and receiver function analysis for IR7, one of the best-fitting models from the primary inversion result is selected (Fig. 16a). The results of forward and inverse modelling of the stacked radial receiver function based on the adjusted model are illustrated in Fig. 16b and c. The results of the forward and inverse modelling are similar and show a good match to the three main phases (i.e., Pp, Ps, and PpPms). Therefore, the corresponding velocity – depth models (Fig. 16d) were considered to reliably represent the main features of the crustal structure for IR7. The structural model for IR7 can be thus summarised as: (1) an upper crust with P-wave velocity of 5.6 km/s down to a 13-km depth; (2) a middle crust with P-wave velocity of f 6.3 km/s down to a 30-km depth; (3) a lower crust with P-wave velocity of f 6.7 km/s down to a 47-km depth; and (4) velocity gradient from f 6.8 to 8.05 km/s down to a depth of about 60 km.
Acknowledgements We thank S. Mangino and an anonymous reviewer for their valuable comments and suggestions that improved the paper. We are also grateful to Prof. M. Ghafory-Ashtiani, the president of IIEES, for his kindness and help during revision of the paper. The MSRT of Iran has supported JD.
References Alavi, M., 1996. Tectonostratigraphic synthesis and structural style of the Alborz mountain system in northern Iran. J. Geodyn. 21, 1 – 33. Amini, X., Emami, N.H., 1992. Geological map of the Tehran sheet: scale 1:100,000. Geol. Surv. Iran (Tehran, Iran). Ammon, J.C., 1991. The isolation of receiver function effects from teleseismic P waveforms. Bull. Seismol. Soc. Am. 81, 2504 – 2510. Ammon, J.C., Randal, G.E., Zandt, G., 1990. On the nonuniqueness of receiver function inversion. J. Geophys. Res. 95, 15303 – 15318. Asudeh, I., 1982. Seismic structure of Iran from surface and body wave data. Geophys. J. R. Astron. Soc. 71, 715 – 730. Berberian, M., Yeats, R.S., 2001. Contributions of archaeological data to studies of earthquakes history in the Iranian plateau. J. Struct. Geol. 23, 563 – 584. Berteussen, K.A., 1977. Moho depth determination based on spectral-ratio analysis NORSAR Long-Period P-waves. Phys. Earth Planet. Inter. 15, 13 – 27. Chen, C.Y., Chen, W.P., Molnar, P., 1980. The uppermost mantle P wave velocities beneath Turkey and Iran. Geophys. Res. Lett. 7, 77 – 80. Darbyshire, F.A., Priestley, K.F., White, R.S., Stefansson, R., Gudmundsson, G.B., Jackobsdottir, S., 2000. Crustal structure of central and northern Iceland from analysis of teleseismic receiver functions. Geophys. J. Int. 143, 163 – 184. Dehghani, G.A., Makris, J., 1984. The gravity field and crustal structure of Iran. Neues Jahrb. Geol. Palaeontol. 168, 215 – 229. De Martini, P.M., Hessami, K., Pantosti, D., D’Addezio, G., Alinaghi, H., Ghafory-Ashtiani, M., 1998. A geological contribution to the evaluation of the seismic potential of the Kahrizak fault (Tehran, Iran). Tectonophysics 287, 187 – 199. Giese, P., Makris, J., Akasheh, B., Roewer, P., Letz, H., Mostaanpour, M., 1984. The crustal structure in southern Iran derived from seismic explosion data. Neues Jahrb. Geol. Palaeontol. 168, 230 – 243. Hedayati, A., Brander, J.L., Berberian, M., 1976. Microearthquake survey of Tehran region, Iran. Bull. Seismol. Soc. Am. 66, 1713 – 1725. Langston, C.A., 1977. Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves. Bull. Seismol. Soc. Am. 67, 713 – 724. Langston, C.A., 1979. Structure under Mount Rainier, Washington,
J. Doloei, R. Roberts / Tectonophysics 364 (2003) 115–133 inferred from teleseismic body waves. J. Geophys. Res. 84, 476 – 4749. Mangino, S., Priestley, K., 1998. The crustal structure of the southern Caspian region. Geophys. J. Int. 133, 630 – 648. Owens, T.J., 1987. Crustal structure of the Adirondacks determined from broadband teleseismic waveform modelling. J. Geophys. Res. 92, 6391 – 6401. Owens, T.J., Zandt, G., Taylor, S.R., 1984. Seismic evidence from an ancient rift beneath the Cumberland Plateau, Tennessee: a detailed analysis of broadband teleseismic P waveforms. J. Geophys. Res. 89, 7783 – 7795.
133
Sobouti, F., Arkani-Hamed, J., 1996. Numerical modelling of the deformation of the Iranian plateau. Geophys. J. Int. 126, 805 – 818. Volarovich, M.P., Volynets, L.N., Levshin, A.L., 1980. Composition of the Earth’s crust and the upper mantle from DSS data in central Turkmenia and laboratory measurements for high thermodynamic parameters. Izv. Earth Phys. 16, 441 – 446. Zandt, G., Myers, S.C., Wallace, T.C., 1995. Crustal and mantle structure across the basin and range – Colorado Plateau boundary at 37jN latitude and implications for Cenozoic extensional mechanism. J. Geophys. R. 100, 10529 – 10548.