Thickness of the lithosphere beneath Turkey and surroundings from S-receiver functions

We analyze S-receiver functions to investigate variations of lithospheric thickness below the entire region of Turkey and surrounding areas. The teleseismic data used here have been compiled combining all permanent seismic stations which are open to public access. We obtained almost 12 000 S-receiver function traces characterizing the seismic discontinuities between the Moho and the discontinuity at 410 km depth. Common-conversion-point stacks yield wellconstrained images of the Moho and of the lithosphere– asthenosphere boundary (LAB). Results from previous studies suggesting shallow LAB depths between 80 and 100 km are confirmed in the entire region outside the subduction zones. We did not observe changes in LAB depths across the North and East Anatolian faults. To the east of Cyprus, we see indications of the Arabian LAB. The African plate is observed down to about 150 km depth subducting to the north and east between the Aegean and Cyprus with a tear at Cyprus. We also observed the discontinuity at 410 km depth and a negative discontinuity above the 410, which might indicate a zone of partial melt above this discontinuity.


Introduction
The lower boundary of the lithospheric plates is a very important parameter for understanding plate tectonics, although it is still one of the less known quantities. Our current knowledge of the lithospheric thickness beneath Turkey relies on studies that examined the data from several temporary and permanent seismic networks (e.g. Angus et al.,20 2006; Sodoudi et al., 2006Sodoudi et al., , 2015Gök et al., 2007Gök et al., , 2015Vanacore et al., 2013;Vinnik et al., 2014). Interpretations from these studies are either confined to a limited region or to a limited depth extent, i.e., to crustal depths only. Thus, the variations of lithospheric thickness have not yet been homogeneously characterized in Turkey and surroundings. A robust estimate of the lithospheric thickness is an important constraint for our under-

Tectonic setting
The Cenozoic closure of the Tethys Ocean and the following continental collision of the northward moving African and Arabian plates with the Eurasian plate are considered as the main driving processes of the present tectonic setting in the eastern Mediterranean and Anatolia (Şengör and Kidd, 1979;Dewey and Şengör, 1979;Taymaz et al., 10 1990Taymaz et al., 10 , 1991aTaymaz et al., 10 , b, 2004Jackson et al., 1992;Armijo et al., 1999;Faccenna et al., 2014;Schildgen et al., 2014). Figure 1 summarizes the tectonic setting of the target area. GPS measurements have highlighted that the retreating Hellenic trench and associated slab roll-back of the subducting African lithosphere is currently the major driving force for the western movement of the Anatolian plate (e.g. McKenzie, 1978;Reilinger 15 et al., 2006). Current deformation is mainly accommodated through strain localization along the North Anatolian Fault (NAF) and the East Anatolian Fault (EAF). The deformation along the EAF is purely left lateral with no compressional component except localized thrust ridges associated with strike-slip tectonics (Bulut et al., 2012). Body wave tomography images on global and regional scale indicate that the Hellenic sub-Introduction   Özacar et al., 2008;Bécel et al., 2009Bécel et al., , 2010Vanacore et al., 2013;Karabulut et al., 2013;Sodoudi et al., 2015). The lithosphere-asthenosphere boundary beneath the East Anatolian Accretionary Complex (EAAC) has been observed at 60-80 km depth based on S-receiver functions (Angus et al., 2006). They consider this lithospheric thickness to be anomalously thin and interpret this structure to be the 5 remnant of the detachment of an oceanic slab. Further evidence for a thin lithosphere overlying hot and partially molten asthenospheric material beneath Eastern Anatolia has been found in other P-and S-receiver function images (Zor et al., 2003;Özacar et al., 2008;Vanacore et al., 2013). Other indications of a shallow LAB are strong attenuation of Lg and Sn phases (Gök et al., 2003) and relatively low Pn-and uppermost mantle S-velocities (Al-Lazki et al., 2004;Maggi and Priestley, 2005). The slow velocity anomaly was attributed to the ascending asthenosphere that is emplaced beneath the plateau following the detachment of the northward subducting Arabian oceanic lithosphere (Keskin, 2003;Faccenna, 2003;Şengör et al., 2003;Biryol et al., 2011;Fichtner et al., 2013a, b). Recently, tomography studies by Biryol et al. (2011) and Fichtner 15 et al. (2013a) have confirmed that a hot and buoyant asthenospheric body supports the ∼ 2 km elevation of the Eastern Anatolian Plateau in the presence of an about 45 km thick crust (Şengör et al., 2003;Zor et al., 2003). At about 350 km depth, a fast anomaly has been found by tomographic studies and interpreted to represent a slab detachment (Lei and Zhao, 2007;Zor, 2008).

Data and method
Recent extensions of seismic and geodetic networks in Turkey have considerably contributed to our understanding of seismicity patterns and crustal and lithospheric structures beneath Anatolia and the surrounding regions. These developments have been triggered in particular by the destructive 1999 Izmit (Mw 7.4) and Düzce (Mw 7.2) earthquakes. Our data set consists of teleseismic waveforms that were extracted from 1028 earthquakes recorded at 195 broadband stations (Fig. 2 GEOFON (Hanka and Kind, 1994;GEOFON Data Centre, 1993). We apply the S-receiver function method to this dataset in order to image seismic discontinuities in the crust and upper mantle. A description of the S-receiver function 5 analysis scheme and examples are found in Yuan et al. (2006) and Kind et al. (2012Kind et al. ( , 2014. The method employs teleseismic S phases, which are converted to P waves at a discontinuity below the station site.The main advantage of this approach over the Preceiver functions techniques is that there is practically no interference from multiples, which can otherwise significantly interfere with the imaging of deeper convertors, in 10 particular for the depth range of 80-200 km in continental regions. S-receiver functions provide a broad sampling of the upper mantle comparable to the steep angle reflection images of the crust. This technique requires a sufficient number of high quality observations of teleseismic S phases, typically at least about 60 records per station with signal-to-noise ratio above two. Permanent broadband stations or densely deployed 15 temporary broadband networks can usually provide this amount of data. The limit of resolution is determined by the frequency content of teleseismic S phases, which usually only give convincing results for periods longer than a few seconds. Since the essential part of S-receiver functions are weakly converted P-precursors of the SV phase, the signal-to-noise ratio needs to be improved by stacking large number of records at 20 the same station, or within a confined region of conversion points at a fixed depth, or by migrating the rays into the depth domain based on a 1-D reference model and the assumption that conversions occur at horizontal boundaries. The latter approach is known as a common-conversion-point stack (CCP). To reduce source side effects, the distribution of sources should have a well-balanced azimuthal and epicentral coverage. In We have restricted the dataset by visual inspection to reasonably simple and clear waveforms. Deconvolution results that fail to transform the SV signal (Q component) 5 into a sharp spike, have also been excluded based on a further visual inspection. This resulted in 11 660 S-receiver functions. For initial evaluation, all traces are combined within non-overlapping 0.5 • epicentral distance bins irrespective of the station and source locations (Fig. 3). This figure shows all observable precursors of the S phases, including phases that might disturb the expected S-to-P conversion from upper mantle 10 discontinuities below the stations. A positive conversion (red) at −3 s, corresponding to the Moho, and a negative (blue) signal at −10 s, interpreted as the LAB, is clearly visible. We can also identify the SKS660p and ScS660p phases caused by S-to-P conversions of SKS and ScS at the 660 km discontinuity. These signals have significantly different slowness compared to the S-to-P conversions of the S phase at the LAB or 15 Moho, meaning that they are cutting through all other phases before the S signal. It means that, if we apply proper delays before summation according to the slowness of the converted S phase, all other phases traveling with different slowness are suppressed. The strong signal following S is the SPn phase, which is caused by an S-to-P conversion traveling horizontally as a P wave through the uppermost mantle layer at 20 the receiver side.

Synthetic tests for imaging inclined structures with S-receiver functions
Since we expect to image an inclined slab in the south of the study region, we first use synthetic seismograms to quantify the magnitude of migration artifacts and infer the maximum range of inclination at which the LAB of a slab can be imaged using Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | receiver functions for inclined layers. In these figures, the converters used for the forward modeling are displayed by black lines. Positive and negative amplitudes of the synthetic S-RFs (S-receiver functions) are shown by red and blue, respectively. Based on available reference tomographic model (e.g. Biryol et al., 2011;Bakirci et al., 2012), an inclined 45 km thick high velocity zone (slab) dipping to the north is used as input  Figure 4g shows the emergent angle of the Sp converted wave at the lower boundary of the slab relative to the slab normal for the corresponding inclination angles of the slab as a function of the back azimuth. For zero and small incidence 15 angles of the high velocity zone, the discontinuities are properly imaged at expected positions (see Fig. 4a and b). A signal is generated by events from nearly all azimuths. For a slab inclination greater than 20 • , the lower boundary of the slab is imaged only by events from southern directions since for back azimuths of 90 • the emergent angle of the converted Sp wave reaches its critical value (see green line in Fig. 4g). Thus, 20 for events from northern azimuths (0-90 • ), no transmitted Sp conversion exists. For larger inclinations of the slab, the back azimuthal range for which Sp converted energy is transmitted becomes even smaller (see Fig. 4d-f and red, yellow and violet lines in Fig. 4g). Moreover, the amplitudes from the inclined structure show a reversal for slab inclinations larger than 30 • and events from southern directions (160-180 • ). This am-25 plitude reversal occurs similarly for the Q-component of P-RFs from inclined structures (Schneider et al., 2013) and might cause destructive interference in CCP stacks. 7,2015 Thickness of the lithosphere beneath Turkey and surroundings from S-receiver functions R. Kind et al.

Observations and interpretation
We extracted several north-south and east-west profiles through Anatolia. The results from time domain stacks and depth-domain migration are shown in Figs. 5-10. The receiver functions are filtered with an 8 s low-pass filter in those figures. In the time domain, a move-out correction is applied using a reference slowness of 6.4 s •−1 and 5 time-corrected traces are stacked for S-to-P piercing points at 200 km depth within a certain window of latitude and longitude along the selected profiles. We chose the piercing point depth of 200 km in the time domain images in order to focus on potentially subducting structures. The summation windows have been placed along east-west and north-south profiles. The size of the windows is 0.3 • in the north-south (or east-10 west) direction and the width of the section in the east-west (or north-south) direction given in each figure caption (e.g. Fig. 5). Neighboring boxes are not overlapping, i.e., no lateral smoothing is applied along the profiles. The relatively large profile widths (3-9 • ) were required to collect a sufficient number of traces to ensure a good signalto-noise ratio of the originally weakly converted phases in the summation trace. In 15 order to obtain an image directly in the depth-distance domain, we also calculated CCP stacks, which, unlike the time domain summation, does not require the arbitrary choice of a piercing point depth. We assumed the IASP91 global reference model for the time-to-depth conversion. The receiver function migration technique does not depend strongly on the given velocity model.

Moho
A positive converter associated with the Moho is very obvious in all profiles in Figs. 5-10. The Moho is only briefly discussed here as it is not the main target of this study. Basically, shorter period P-receiver functions are more suitable for a detailed image of the 5 Moho or intracrustal discontinuities. However, we point out a few crustal observations obtained from S-receiver functions. The crustal thickness increases from the Aegean to Eastern Anatolia from ∼ 25 to ∼ 40 km (see Fig. 9). As the seismic traces are averaged with piercing points in relatively large grids (1 • longitude, 4 • latitude in Fig. 9), smallscale variations of the Moho are not well-resolved. The observed range of the Moho 10 depths is consistent with results from previous receiver function studies (e.g. Saunders et al., 1998;Zor et al., 2006;Özacar et al., 2008;Vanacore et al., 2013;Karabulut et al., 2013). However, in Eastern Anatolia, Şengör et al. (2003) and Zor et al. (2003) reported that the Moho depth increases from less than 40 km at the Bitlis Suture to about 50 km near the Black Sea. This deep Moho is not seen in our data (Fig. 7), im- 15 plying that this change in Moho depth only occurs in a limited region, such that it is not visible in our relatively large-scale north-south profile in Fig. 7. Vanacore et al. (2013) presented a detailed Moho map from P-receiver functions. They also observed a general increase in Moho depth from west to east, which is roughly consistent with our S-RF results. The subducting African Moho is visible to a depth of nearly 100 km in 20 the southern part of western Anatolia and the Aegean (Fig. 5a). No significant change in Moho depth is visible across the North Anatolian Fault in the south-north profiles (see e.g. Fig. 6). However, our data are filtered with an 8 s low-pass filter and therefore not optimized for higher resolution determination of the Moho depth across the NAF. P-receiver function are more useful for this purpose.

Lithosphere-asthenosphere boundary
The stacked S-RF images in Figs. 5-10 show a clear observation of a relatively homogeneous negative discontinuity at 80-100 km depth beneath the entire region, which we interpret to represent the LAB. The smooth appearance of the LAB is due to the chosen filter (8 s low-pass). Using shorter wavelengths, the LAB appears to be less 5 well defined (see Fig. 11). This indicates that the LAB may consist of several neighboring discontinuities and due to the superposition of neighboring signals, the depth of the LAB may depend on the period. We have chosen this relatively long period filter to emphasize the dominant structures in the mantle as clearly as possible.
Our observations of a shallow LAB beneath Anatolia confirms earlier results. Previous findings mainly from body and surface wave tomography studies (Biryol et al., 2011;Bakirci et al., 2012;Salaün et al., 2012;Fichtner et al., 2013a, b) indicate a similarly thin lithosphere beneath entire Anatolia. Additional constraints confirming a thin lithosphere are provided by the much smaller-scale P-RF and S-RF images in Sodoudi et al. (2006Sodoudi et al. ( , 2015, Angus et al. (2006), Özacar et al. (2008) and Gök et al. (2011).

Subduction of African lithosphere
In our stacked S-RFs and corresponding depth migrated images, the African LAB dips flatly to the north with an angle of 20-30 • and reaches down to almost 150 km depth 25 beneath Anatolia and the Aegean (Figs. 5 and 6). This is similar to the observed slab geometry found by a recent Rayleigh wave tomography by Bakirci et al. (2012) also similar to receiver function results from Li et al. (2003) and Sodoudi et al. (2006Sodoudi et al. ( , 2015. In contrast to that, the African slab observed with tomographic methods reaches into the lower mantle with a dip of 40-50 • (e.g. Wortel and Spakman, 2000). What is the reason for the differences in the observations of the two types of data? Comparing our data with the synthetics in Fig. 4, we realize that for a slab dip > 40 • only events with 5 a backazimuth between 160-200 • would contribute to the slab images. However, we have only very few events from the south (see Fig. 2). A steep part of the slab is thus not expected to show on our images, and its absence on our image does not rule out the presence of a deeper slab. Figure 4 also shows that slabs with a dip of 20-30 • can be observed relatively well with S-receiver functions. All our images of the subducting 10 African plate are only obtained from the shallow and flat dipping part of the slab. Southward retreat of the Aegean slab has long been considered as the major mechanism to explain the N-S oriented extensional regime in western Anatolia and the Aegean Sea, as well as the rapid deformation and westward extrusion of Anatolia along the North and East Anatolian Faults (NAF and EAF) (McKenzie, 1978;Okay and Satır, 15 2000; Reilinger et al., 2006). The slab retreat cannot be imaged directly in this study. However, our observations of the thin lithosphere provide seismological support for the removal of lithosphere beneath entire Anatolia. Indications of lithospheric removal in the region have been given in earlier regional tomography studies (i.e. Wortel and Spakman, 2000;Piromallo and Morelli, 2003). Şengör et al. (2003) proposed the de-20 tachment of the Arabian plate and subsequent upwelling of asthenospheric mantle to explain the high elevation and the lack of a thick mantle lithosphere beneath the Eastern Anatolian Plateau. Features of slab break-off were inferred later from teleseismic tomography images (Lei and Zhao, 2007;Zor, 2008). Without requiring slab break-off models, recent geodynamic models (Göğüş and Pysklywec, 2008;Komut et al., 2012;25 Göğüş, 2015) associate present-day high topography, lithospheric thinning and hot surface heat flow with vertical flow and very high thermal conditions (e.g. 1300-1400 • C) of the uppermost mantle that resulted in the delamination of lithosphere, which was originally proposed by Bird (1979) Biryol et al. (2011) observed a wide gap between the high velocity anomalies of the subducting Aegean and Cyprus slabs, which they interpreted as slab tear. Our northsouth profile at the transition of the Aegean and Cyprus slabs (Fig. 6) is too broad to resolve such a slab tear. The dipping features of the African Moho along the Hellenic and Cyprus trenches differ, implying that the subduction is only shallow toward the 5 south of Cyprus (Figs. 5 and 6). The absence of trench-perpendicular extension on the over-riding central Anatolia is different from the western Anatolia and the Aegean and was previously considered to be the consequence of a lower subduction angle of the African slab along the Cyprus trench (Barka and Reilinger, 1997). Changing slab geometry with decreasing dip angle from west to east was previously suggested to result from the Eratosthenes Sea Mountain beneath Cyprus resisting subduction (Kempler and Ben-Avraham, 1987; Barka and Reilinger, 1997). The east-west profile along the southern boundary of Anatolia (Fig. 8) provides more detailed information.
Here, we see that the African (or Aegean) slab is apparently dipping to the east towards Cyprus where it reaches about 150 km depth; this geometry can be interpreted as 15 a northward-directed dip of the African Plate, resulting in a structural boundaries that are inclined to the north and east. To the east of Cyprus, we see the flat Arabian LAB near 100 km depth (Fig. 8), clearly discontinuous with the African slab LAB. This transition occurs significantly west of the Dead Sea Transform Fault (which is at 35-36 • E, see Fig. 1), which might indicate that the Dead Sea Transform Fault is inclined 20 to the west. we also do not observe any considerable contrast in the delay times of the LAB signals across the NAF (Fig. 6). We conclude therefore that lithospheric depths and average velocities are very similar on both side of the fault. We also tested narrower northsouth profiles and obtained similar results, although the data quality is worsening when fewer traces are summed. The oblique pattern of the fast polarization direction inferred 5 from the SKS measurements to the actual orientation of the NAF is another indication that the NAF is a relatively shallow feature not linked to a deep deformation zone (e.g. Paul et al., 2014). These observations are different from those found in northeast Tibet (Leon Soto et al., 2012;Eken et al., 2013) where fast anisotropy directions within the entire lithosphere are nearly parallel to the strike of the North and South Kunlun faults.

Depth of the NAF
In contrast to northern Anatolia, significant differences in crustal and lithospheric thickness are observed across the San Andreas Fault in southern California based on P-and S-receiver functions (Lekic et al., 2011;Miller et al., 2014).

Mantle transition zone discontinuities
We observed the 410 discontinuity beneath most of Anatolia. This discontinuity is shal-15 lowest beneath Western Anatolia at ∼ 390 km depth and deepens smoothly towards central Anatolia to 420 km. The depth is approximately constant below central Anatolia and deepens again across eastern Anatolia from 410 to 430 km (Fig. 9). We leave this observation to a further analysis using P-receiver functions that are more suitable for this purpose. In the context of this study, these observations are mainly used to verify 20 the reliability of our data. Another observation is a negative, but more scattered discontinuity about 50 km above the 410 in the entire area. Such a discontinuity is known globally (Tauzin et al., 2010) and considered to be caused by a low velocity layer containing liquids or partial melts. We do not observe any other prominent discontinuity between the LAB and the low velocity zone on top of the 410 on larger scale, especially 25 no clear indications for the Lehmann discontinuity, which has been observed below the north-western United States (Kind et al., 2015) Geology, 36, 723-726, doi:10.1130/G24982A.1, 2008. Gök, R., Sandvol, E., Törkelli, N., Seber, D., and: Sn attenuation in the Anatolian and Iranian plateau and surrounding regions, Geophys. Res. Lett., 30, 8042, doi:10.1029/2003GL018020, 2003 Gök, R., Pasyanos, M., and Zor, E.: Lithospheric structure of the continent-continent collision zone: eastern Turkey, Geophys. J. Int., 169, 1079-1088, doi:10.1111/j.1365-246X.2006.03288.x, 2007 Teoman, U., Turkelli, N., Godoladze, T., and Javakishvirli, Z.: Lithospheric velocity struc-  Lithos, 120, 14-29, 2010. Kalafat, D., Öütcü, Z., Suvarikli, M., Güneş, Y., Yilmazer, M., Kekovali, K., Kara, M., Kiliç, K., Bekler, F. N., Küsmezer, A., Deniz, P., Çomolu, M., Berberolu, M., Berberolu, A., Poyraz, S. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Rychert, C. A. and Shearer, P. M.: A global view of the lithosphere-asthenosphere boundary, Science, 324, 495-498, doi:10.1126/science.1169754, 2009. Rychert, C. A., Shearer, P. M., andFischer, K. M.: Scattered wave imaging of the lithosphereasthenosphere boundary, Lithos, 120, 173-185, 2010. Salaün, G., Pedersen, H. A., Paul, A., Farra, V., Karabulut, H., Hatzfeld, D., Papazachos, C., 5 Childs, D. M., Pequegnat, C., and SIMBAAD Team: High-resolution surface wave tomography beneath the Aegean-Anatolia region: constraints on upper-mantle structure, Geophys. J. Int., 190, 406-420, doi:10.1111/j.1365-246X.2012.05483.x, 2012 Tan Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Şengör, A. M. C. and Kidd, W. S. F.: The post-collisional tectonics of the Turkish-Iranian Plateau and a comparison with Tibet, Tectonophysics, 55, 361-376, 1979. Şengör, A. M. C., Özeren, S., Genç, T., and Zor, E.: East Anatolian high plateau as a mantle supported, north-south shortened domal structure, Geophys. Res. Lett., 30, 8045, doi:10.1029/2003GL017858, 2003.      , 2000) for different inclination angles of a high velocity zone (slab). The converters of the input models are depicted as black lines. The velocities used for the modeling are given in a). The locations of the converters as reconstructed by CCP stacking of the synthetic S-receiver function amplitudes are shown in red (positive) and blue (negative). As for the data migrations, a 1-D velocity model is used for calculating the conversion points, resulting in the offset between input and recovered anomalies for inclined discontinuities. g) The emergent angles (relative to the slab surface) of the Sp converted and transmitted waves are calculated in dependence of the back azimuth of the incident S-wave for the velocity contrast at the lower boundary of the slab. The results for the different considered inclination angles (dip) are shown in different colors. The range of back azimuths for which Sp converted energy is transmitted is getting smaller for increasing dipping angle. 7,2015 Thickness of the lithosphere beneath Turkey and surroundings from S-receiver functions At shorter periods (a and b) the LAB (blue signal around 10 s) may consist of several smaller discontinuities. Only for the 6 and 8 s low-pass filters it appears as single discontinuity. In contrast, the Moho and the discontinuity at 410 km depth (red signals around 5 and 45 s, respectively) remain as sharp single discontinuities at all periods, suggesting that this behavior is due to structure rather than to the higher noise levels at shorter periods.