APPLICATION OF RECEIVER FUNCTIONTECHNIQUE TO WESTERN TURKEY
Timur TEZEL∗
MEE07169

Supervisor:

Takuo SHIBUTANI∗∗

ABSTRACTIn this study I tried to determine the shear wave velocity structure in the crust and uppermost upper
mantle using receiver function technique. For this purpose I selected teleseismic earthquakes in theepicentral distances between 30 and 90 degrees, magnitudes of which are greater than 6.3. I collectedwaveform data recorded at Turkey’s broadband seismic stations between 2006 and 2008. I used radialreceiver functions which were calculated using extended time multitaper method to determine thecrustal and uppermost upper mantle structure. My study consists of three steps: receiver function (RF)calculation, RF image and RF inversion. Because of the data availability, I studied data of 33 stations
in western Turkey which belong to the General Directorate of Disaster Affairs Earthquake Research
Department and Bogazici University Kandilli Observatory. First, I calculated radial receiver functions
for each station. In the next step, I made two RF images using the RFs from stations along two lines:
the one is running from the northeast of Bozcaada to Konya city and the other is running from thenorth of Izmir to the northeast corner of Rhodes island. In the third step, I applied genetic algorithm
inversion method to determine the shear wave velocity structure beneath six seismic stations.
The results show that the Moho depth changes from region to region; we observed that theMoho depth is around 35 km and 25 km beneath the former and latter profiles, respectively. The
results of receiver function inversion are consistent with these RF images. Shear wave velocities were
estimated to be in the range between 3.5 and 3.9 km/s for the lower crust and between 4.0 and 4.6
km/s for the uppermost upper mantle, respectively.

Keywords: Receiver Function, Crust, Teleseismic Wave, Turkey.
INTRODUCTION
Turkey has experienced many natural disasters which have caused serious casualty, collapse ofbuildings, economical losses up to now. Among these disasters, earthquakes are the first to beconcerned rather than land slides, floods, rock falls, avalanches, and droughts. In Turkey, earthquakes
occur in the upper crust generally, and one of the important issues is to determine discontinuities andvelocity changes in the crust and also to determine the transition to the Moho. There are someseismological studies to find a velocity model using travel time tomography, surface wave groupvelocity inversion, reflection profiling, and receiver function method etc. Determination of the
discontinuities and especially shear wave velocity structure has an important role in the planning of theurban areas and cities. In this study, we applied receiver function method to determination of shear
wave velocity structures beneath the seismic stations in Turkey.

∗General Directorate of Disaster Affairs Earthquake Research Department, Ankara, Turkey.
Assoc. Professor, Disaster Prevention Center, Research Center of Earthquake Prediction, Kyoto University,Japan.

∗∗1

DATA
I used teleseismic waveform data for 50 earthquakes recorded by General Directorate of DisasterAffairs Earthquake Research Department (here after ERD) and Boğaziçi University Kandilli
Observatory (KOERI) broadband stations. I selected 33 broadband stations which are located inwestern Turkey. Selection criteria were as follows: magnitude of earthquakes should be greater than
6.3 and epicentral distance should be between 30˚ and 90˚.
Data were obtained from KOERI via internet (http://barbar.koeri.boun.edu.tr) and from ERD
via CD-R.

THEORY AND METHODOLOGYThe Receiver Function technique is used for determination of crustal structure beneath seismic stations.
Teleseismic records (30˚ ≤ ∆ ≤ 90˚) are used in this method. Teleseismic P waveforms contain
information related to source, propagation path and local structure beneath the recording station. The
method uses the coda part of the P wave which includes converted phases and reverberations
generated at discontinuities beneath each station and convolved with source function and instrumental
impulse responses. If we eliminate the source and instrumental effects from waveforms, they provide
information about the local velocity structure under the seismic station. In Figure 1, we can see thesimple ray diagram for the incident wave and its converted and reverberated phases.

Three componentstation
Free Surface

P waveS waveInterface h
PpPhs PpPhs+PpShs

Pp Ps

Figure 1. Simplified ray diagram identifying the major P to S converted phases which comprise thereceiver function for a single layer over a half-space under the seismic station (modified fromAmmon(1990)).
A teleseismic P wave arrives at the recording station with both a constant and relatively large
horizontal phase velocity (15-25 km/s). This phase velocity justifies a plane wave and simplifies thestudy of the resulting ground motion (Ammon, 1990).
Receiver Function Calculation
Determination of the velocity structure of the crust and upper mantle beneath a single seismic station
can be done by using teleseismic receiver functions. The receiver function analysis uses the converted
phases and multiples recorded on the horizontal seismograms (e.g. Burdick and Langston, 1977;Langston 1977, 1979; Owens et al., 1984; Ammon, 1991).

2DV (t ) = I (t ) * S (t ) * EV (t )DR (t ) = I (t ) * S (t ) * E R (t )(1)

DT (t ) = I (t ) * S (t ) * ET (t )Here, V, R, T shows vertical, radial and tangential components respectively. Also, I(t) isimpulse response of the recording instrument, S(t) is the seismic source function, EV(t), ER(t), ET(t) are
the vertical, radial and tangential impulse response of the earth structure. And asterisk (*) shows theconvolution operator.
As earth structure beneath a station will produce phase conversions of the P to S type,horizontal components of ground motion will be different from the vertical component. DV(t) contains
the factors which we wish to remove from observed seismograms, so isolating ER(t) and ET(t) can beaccomplished by deconvolving DV(t) from DR(t) and DT(t). In receiver function method, there areseveral stabilization methods; they are the water level method, the multitaper method, the extended
time multitaper method. In this study I used the extended time multitaper method which was improvedby Shibutani et al. (2008) based on Park and Levin (2000) and Helffrich (2006). In this method, theyused the three lowest-order 4π prolate eigentapers of 50 sec duration. They summed the multitapers
with 75 % window overlap. The multitapers smoothly are connected at the junctions, and the resultant
taper has a flat level. If we use the flat part for windowing the P onset and the P coda, we can estimate
the relative amplitudes of receiver functions (Shibutani et al., 2008). As a result, the frequency-domain
receiver function is defined byK −1

E R (ω ) =∑ H

k =0(k ) (ω )U~ k (ω )2K −1

⎛ J ⎞ K −1 (k ) ~ k(k ) ~ (k )∑ U (ω )U (ω ) + ⎜⎜ S ⎟⎟ ∑ N (ω )N (ω )⎝ J N ⎠ k =0
k =0

G (ω )(2)

where H (k ) (ω ) and U (k ) (ω ) denotes the Fourier spectrum of the radial or transverse and verticalcomponents of the waveform data with the kth prolate eigentaper, and N (k ) (ω ) is the Fourier spectrumof presignal noise of the vertical component waveform. JS and JN are the numbers of the multitapers
used for the signal and presignal. The second term in the right-hand side of equation (2) is a Gaussianhigh cut filter in which a controls the the corner frequency. In this study, I set a to 2, and then thecorner frequency becomes 0.3 Hz. I showed calculated radial receiver functions with this method fortwo stations (GDZ and KDHN) in Fig. 2. In these figures, we can see the difference between
waveforms of receiver functions related with backazimuth. It reflects the complexity beneath thestations and the lateral change of earth structure.

Figure 2. Calculated receiver functions for the GDZ and KDHN stations are shown on the left andright panel, respectively. Receiver functions are ordered by backazimuth and positive amplitudes arecolored by red.

3

Receiver Function Image
The peaks and troughs of receiver functions correspond to boundaries of the S-wave velocity structure,and the time axis of the receiver functions can be converted to the depth axis using a 1-D velocity
model. Delay time between converted Ps phase and direct Pp phase is given by the below equation,z=Z⎛⎞
11T ps ( Z ) = ∫ ⎜− p2 −− p 2 ⎟dz
2⎜⎟α (z )2z = 0 ⎝ β (z )⎠(3)where α and β are P and S wave velocities as a function of depth z respectively, Z indicates the depthof the interfaces, and p is the ray parameter. The converted receiver functions can be represented by abending ray with the ray parameter and a backazimuth. I projected the rays into 1 km by 1 km cells. Iused four 1-D velocity models: the three models W1, W2, and W3 are based on Tezel et al. (2007) andthe fourth model is AK135 (Kennett et al. 1995). Figure 3 shows the 1-D velocity models (W1, W2,W3 and AK135) and two lines: the one is running from the northeast of Bozcaada to Konya city andthe other is running from the north of Izmir to the northeast corner of Rhodes island along which RF
images are processed. Figure 4 shows the images obtained by this method. In the figures, the positive
amplitudes of the receiver functions are indicated by the red color and indicate velocity discontinuities
at the top of high velocity layers. In Figure 4, 2-D RF images reflect the relation between earth
structure and 1-D velocity models W1, W2, W3 and AK135 respectively.

Figure 3. Input 1-D velocity models based on Tezel et al. (2007) and Kennett et al. (1995) and selected
image lines ( A-AA and B-BB).
Receiver Function Inversion
In previous studies, a linearization procedure was used to invert the receiver function which requires
the initial model to be close to the true velocity structure. Ammon (1990) showed that the final models
were dependent on the initial models. For this reason in this study I used the genetic algorithm (GA)(Shibutani et al. 1996).
In this study, I applied this inversion technique to radial receiver functions of the six
broadband stations. For this purpose radial receiver functions were stacked according to thebackazimuths and RF waveform similarities for each station. I tried to model the crust and uppermost
mantle down to 50 km with six major layers: a sediment layer, basement layer, upper crust, middle
crust, lower crust and uppermost mantle. The model parameters in each layer are the thickness, the Swave velocity. The velocity ratio between P and S waves (Vp/Vs) and the density in the sediment andbasement layers are also model parameters. For each model parameter, upper and lower bounds and 2n4

possible values are specified. The size of the model space to be searched is 246 ≈ 7.04 x 1013. Figure 5
shows the Moho depth estimated for these stations.

DISCUSSIONS AND CONCLUSIONSIn this study I applied RF (receiver function) technique to determine the velocity discontinuities suchas Moho etc. Moreover, I applied genetic algorithm using receiver functions to determine the shear
wave velocity structure beneath six seismic stations. I used only radial receiver functions (I calculated
the tangential receiver functions but not used them in the interpretation). After calculation of receiver
functions, I employed 2-D image technique which converts the time axis to depth axis using 1-D
velocity models and shows the amplitudes with colors. Red (blue) color shows positive (negative)
amplitude which indicates the increase (decrease) in the velocity. Along A-AA and B-BB lines I
applied this method and some velocity discontinuities were detected and indicated by red dotted lines(Fig. 4).
AAA

BBBFigure 4. RF images for the A-AA and B-BB lines using W3 1-D velocity model.
We interpreted that Moho should lie between 30 and 40 km for line A-AA and probably around 35 km.
For some stations we can see another signature around the 15 km, which was interpreted as Conrad
discontinuity. In line B-BB Moho discontinuity seems placed between 20 and 30 km depth andprobably around 25 km.
After RF image process I employed the genetic algorithm to determine the shear wave velocity
model beneath seismic stations. I selected six broadband stations (BALB, GDZ, ALT, KDHN, LADK,KONT) that have more data and show good RF results. I stacked the RFs according to thebackazimuths and then applied inversion process for each station.
Synthetic receiver functions generally show good fitting with the observed waveforms. Shearwave velocity models which were derived from inversion indicate that general character of crust –mantle boundary in the region. The depth of the boundary for each station was determined. Shearwave velocities near the surface varied between 1.5 and 3.5 km/s and at the mid – crustal depth range
between 3 and 4 km/s. Sub Moho velocities are between 4 and 4.6 km/s. The results for some groups(made based on backazimuths) for the stations KONT, ALT and GDZ indicate low velocity zone
between 25 and 35 km.
Figure 5 show the estimated Moho depths for six stations. These velocity values are generally
similar to those of previous studies, while we can not observe low velocity layer around 10 km thatwas mentioned by previous studies.

5Figure 5. Moho depths calculated from inverison
Results of this study with RF image and RF inversion show the Moho thinner in the west of study area
than in the east part. These results similar to the previous studies which were mentioned the Moho
depth changes from 30 km to 38 km from the west to the east of Turkey. Moreover, other stations
which were used for the receiver function show complex receiver waveforms, which reflect thecomplex tectonic structure around the study area. Also, difference between A-AA and B-BB profiles
suggest that the study area is not a homogenous area and that deformation is different in each area. In
these RF images Moho discontinuity has been seen as a wide band caused by the Gaussian filter used
the mean value for the Moho depth considering the error range (i.e. ± 5 km). This study is apreliminary attempt to determine the Moho discontinuity based on RF image and will contribute toconstruction of a reliable reference velocity model in Turkey.

AKNOWLEDGEMENT
I would like to express my sincere gratitude to my adviser Dr. T. Hara for his suggestions andcontributions to my study.

REFERENCES
Ammon, C.J., 1990, J. Geophy. Res., 95, 15303-15318.
Ammon, C.J., 1991, Bull. Seism. Soc. Am., 81, 2504-2510.
Burdick, L.J., and Langston, A.C., 1977, Bull. Seism. Soc. Am., 67, 677-691.
Helffrich, G., 2006, Bull. Seism. Soc. Am., 96, 344-347.
Kennett B.L.N., Engdahl E.R., and Buland, R., 1995, Geophys. J. Int., 122, 108-124.
Langston, C.A., 1977, Bull. Seism. Soc. Am., 67, 1029-1050.
Langston, C. A., 1979, J.Geophys.Res., 84, 4749-4762.
Owens, T. J., Zandt, G., and Taylor, S.R., 1984, J.Geophys. Res., 89, 7783-7795.
Park, J., and Levin, V., 2000, Bull. Seismo. Soc. Am., 90, 1517-1520.
Shibutani, T., Sambridge, M., and Kennett, B., 1996, Geophy. Res. Let., 23, 1829-1832.
Shibutani, T., Ueno, T., and Hirahara, K., 2008, Bull. Seismo. Soc. Am., 98, 812-816.
Tezel, T., Erduran, M., and Alptekin, Ö., 2007, Annals of Geophysics, Vol.50, No.2, pp. 177-190.

6