Variations in the crustal structure beneath western Turkey

Summary

We use teleseismic receiver functions to investigate the crustal structure at two locations in western Turkey using seismic data recorded on small arrays of temporary broad-band seismographs. The results from these analyses are compared with receiver function results from the GDSN station ANTO on the Anatolian Plateau in central Turkey. The crust is ∼ 30km thick in the region of western Turkey where active normal faulting reveals present-day extension in the upper crust and alkali-basaltic volcanism reveals recent extension within the subcrustal lithosphere The crust is ∼ 34km thick further east where crustal extension is still evident but less pronounced. In the Anatolian Plateau, which is not currently extending, the crust is ∼ 38km thick. The level of extension estimated from these measurements of crustal thickness implies a β -factor of ∼ 1.2. This value agrees with the amount of extension estimated in the upper crust from the integrated seismic strain rate (β -factor of ∼ 1.3), from surface faulting(β -factor of ∼ 1.25) and from the amount of extension in the subcrustal lithosphere estimated from the volcanism (β -factor <2), all indicating that the extension is approximately uniformly distributed vertically throughout the lithosphere. The Moho transition in this region appears to thin slightly as the degree of extension increases westwards.

1 Introduction

The Aegean Sea and the surrounding coastal areas of Greece and western Turkey form one of the most seismically active and rapidly deforming continental regions in the world. In this study we investigate the variation in crustal thickness across western Turkey and relate the crustal structure to extension rates as indicated by normal faulting in the upper crust and by extension-related volcanism in the subcrustal lithosphere.

Crustal deformation in western Turkey occurs as a consequence of the westward motion of the Anatolian Plateau and the southwestward motion of the southern Aegean Sea relative to Eurasia (McKenzie 1972, 1978a). The collision of Arabia and Eurasia in the Caucasus (Şengör & Kidd 1979; Dewey et al. 1986) has thickened the crust of the Anatolian Plateau, and the gravitational potential energy stored in the elevated topography, coupled with the continuing northwards motion of Arabia, forces the Anatolian Plateau westwards (McKenzie & Yilmaz 1991). This westward motion of the plateau is accommodated by right-lateral strike-slip motion on the North Anatolian Fault system (Ketin 1948; McKenzie 1972; Şengör, Görür & Şaroğlu 1985) and left-lateral strike-slip motion on the East Anatolian Fault zone (Taymaz, Eyidoğan & Jackson 1991). Seismicity within the Anatolian Plateau is low to the east of ∼ 31 °E, but high to the west. The geomorphology of western Turkey is dominated by a series of E–W-trending, normal-fault-bounded horst and graben structures; the N–S extension inferred from these structures is consistent with regional earthquake focal mechanisms (McKenzie 1972, 1978a; Jackson & McKenzie 1984; Eyidoğan & Jackson 1985). Eyidoğan (1988) used the seismic moment tensors of M_S ≥ 5.5 earthquakes occurring between 1943 and 1983 to estimate ∼ 13.5 mm yr⁻¹ N–S extension across western Turkey, implying an extensional β-factor (McKenzie 1978b) of ∼1.3 in this area.

Paton (1992a,b) studied the relationship between the young (<12Ma) extension-related volcanics and the surface expression of the active faulting in western Turkey. From the tilting of Neogene sediments in the footwall of the faults and from 2-D gravity modelling, he estimated 6–10km of extension across each graben, implying average β -factors of ∼ 1.25 for the upper crust. The extension-related magmas are thought to arise from the stretching of the upper mantle but because these cannot be produced by direct melting of the sublithospheric mantle, their volume does not con-strain the amount of extension. The Kula basalts are the most voluminous (∼ 2.3km³, Bunbury 1992), and their high potassium content makes it unlikely that much melt has underplated the lower crust, yet the high Mg content rules out significant high-level fractionation. The basalts are generated by lithospheric melting and only constrain the amount of extension to be less than β≈ 2 (McKenzie & Bickle 1988). Whether the amount of extension in the upper crust (indicated by the integrated seismic strain or surface faulting) and in the mantle lithosphere (indicated by the volume of volcanism) is matched by thinning of the lower crust is not clear. In the Basin and Range province of western North America, where similar or greater crustal extension has occurred, the lower crust has, at least on a local scale, flowed to maintain its thickness (Gans 1987), rather than thinned in response to the deformation.

The few prior estimates of crustal thickness in Turkey are too imprecise to determine whether the crust is thicker beneath the unstretched regions than beneath the stretched regions. Makris & Vees (1977) and Makris (1978) used topography, gravity and limited seismic refraction data in the Aegean Sea to estimate that the crust thickens from about 22km in the central Aegean Sea to over 40km in the Anatolian Plateau. This change in crustal thickness coincides with the decrease in the amount of extension inferred from the eastwards decrease in seismicity and normal faulting in western Turkey, but is not precise enough to reveal how the lower crust has responded to this extension. In this study we analyse teleseismic receiver functions to determine the crustal structure at two sites in western Turkey (Fig. 1): at Kula, where large active normal faults nearby indicate significant present-day extension in the upper crust, and recent alkali-basaltic volcanism suggests recent extension in the lithospheric mantle; and at Us,/e1>14;ak, where fewer faults suggest that some upper-crustal extension is still occurring but much reduced compared to Kula. We compare the crustal thickness at these two sites with the crustal thickness at Ankara in the Anatolian Plateau, where there is no surface evidence for present-day extension.

2 Seismological Constraints On Crustal Structure

We investigate the crustal structure of western Turkey primarily using teleseismic receiver functions. This technique is useful in determining the crustal shear-wave structure beneath a three-component broad-band seismograph. However, a receiver function is sensitive to the relative arrival times of converted phases from discontinuities, but contains little absolute velocity information. Thus, while a receiver function is sensitive to the crustal impedance contrasts, a large range of velocity models may fit the observations equally well. The range of acceptable earth models can be greatly reduced by including in the analysis other seismological data which further constrain the velocities. The apparent velocities of crust and upper-mantle phases from local and regional earthquakes provide rough estimates of the crust and uppermost mantle velocities. However, the main constraint we use is surface-wave phase and group velocities, which are sensitive to the average shear-wave velocity of the crust but weakly constrain the depths to interfaces. Thus surface-wave phase-velocity data complement receiver function analyses.

The seismograms used to determine crustal structure in western Turkey were recorded on small three-station arrays of broad-band seismographs which we installed and operated in the vicinity of Kula and Us,/e1>14;ak (Fig. 1). Each seismograph in the array consisted of a Guralp CMG-3T broad-band sensor recorded on a Refraction Technology 72a–02 16-bit data logger at 10sampless^{− 1}. Absolute timing was provided by an Omega receiver and site locations (¹) were determined by GPS. The seismographs were calibrated and data were collected from the sites at approximately six-week intervals. Once during the deployment the seismographs were calibrated with a random binary input (Bergeret al. 1979). We had hoped that the small three-station array data could be used to enhance the teleseismic signals, but restrictions imposed by site locations and cultural noise made this impossible. Data for the station ANTO at Ankara (Fig. 1) were obtained from the Incorporated Research Institution for Seismology (IRIS) Data Management Center. ANTO consists of a Geotech KS-36000-i borehole seismometer recorded on a model Q-680 24-bit Quanterra data logger.

2.1 Apparent velocities of local and regional phases

The crustal structure of western Turkey is not well known. Several studies have used local and regional earthquake traveltimes to estimate crust and upper-mantle velocities, but with varying results. Chen, Chen & Molnar (1980) measured the P_n velocity in northwest Turkey from arrival-time data at Istanbul (IST) and found a P_n velocity of 7.70 ± 0.04 kms^1-1: ; Kadinsky-Cade (1981) used data from the same station and found a P_n velocity of 8.3 ± 0.3kms^1-1, Necioglu, Maddison & Turkelli (1981) analysed travel-times from 43 earthquakes occurring in western Turkey recorded on a number of regional seismographs and found the P_n velocities to be 7.83 ± 0.17kms^1-1 in the vicinity of Ankara, 8.10 ± 0.26kms^1-1 in the vicinity of Us,/e1>14;ak and 8.04 ± 0.12kms^1-1 in the region north of Kula. Hearn & Ni (1994)) inverted arrival times taken from the catalogue of the International Seismological Centre (ISC) for local and regional earthquakes and found the P_n velocity to be ∼ 7.8kms^1-1 in the Ankara region and ∼ 8.1kms^1-1 in the Kula and Us,/e1>14;ak region. The refraction data of Gürbüz & Evans (1991) from the Tuz Gölü basin 150km SSE of Ankara primarily constrain the velocity of the sedimentary layers, but they also show a 6.15kms^1-1: basement at depths between 6 and 10km. These studies say little about the crustal thickness and allow a large range of P_n velocities.

Saunders (1996) constructed composite record sections from local and regional earthquake seismograms recorded at Kula and Us,/e1>14;ak (Fig. 1). From these, the P_n and L_g group velocities were measured as 8.1 ± 0.1kms^1-1: and 3.5 ± 0.2kms^1-1 respectively at both Kula and Us,/e1>14;ak, but the P_g velocity was 6.4 ± 0.1kms^1-1 at Kula and 6.1 ± 0.1kms^1-1 at Us,/e1>14;ak. ISC traveltime data for all well-located local and regional events with impulsive P_g, S_g and P_n arrivals at station ALT (∼ 60km NE of Us,/e1>14;ak) give 6.1 ± 0.2kms^1-1 for P_g, 3.3 ± 0.2kms^1-1 for S_g and 8.0 ± 0.2kms^1-1 for the P_n phase.

2.2 Surface-wave group and phase velocities

Mindevalli & Mitchell (1989) measured fundamental-mode Rayleigh- and Love-wave group velocities for eastern and western Turkey in the 8–47s period range from ANTO seismograms (Fig. 2). Even with the large uncertainties in their group-velocity data, Mindevalli & Mitchell (1989) feel there are resolvable differences in the crust and upper-mantle structure of western and eastern Turkey.

We have measured Rayleigh-wave phase velocities in the5–100s period range between Kula and Us,/e1>14;ak (Fig. 2). We also intended to measure the dispersion between ANTO and Us,/e1>14;ak but found that for the entire period during which our stations were deployed at Kula and Us,/e1>14;ak the clock at ANTO was in error. The Kula and Us,/e1>14;ak interstation phase-velocity measurements were made using the constrainedleast-squares method of Gomberg et al. (1988)) in which the phase-velocity measurement is considered as a linear filter estimation problem. If S₁(ω) is the Fourier transform of the Rayleigh-wave seismogram at station 1, then S₂(ω), the transformed seismogram at station 2 after propagating the distance Δr along a great-circle path between the stations, is given by

where ω is the angular frequency. The Earth filter F(ω) may be written

where k(ω) is the wavenumber, A(ω) describes the amplitude decay, and the phase velocity, C(ω), is determined from the phase term

A number of seismogram pairs can be used simultaneously to estimate the phase velocity by taking the least-squares solution. The phase-velocity curve is smoothed by limiting the final group velocity to be within some specified range of the starting model group velocity. The phase term is computed by solving for a correction vector to an initial starting phase-velocity curve. The data are weighted in the least-squares solution by a weight matrix consisting of a subjective weight factor, which depends on the relative quality of each pair of recordings and on the squared coherence of the data.

Gomberg . (1988) expanded e^iδkΔr and neglected second- and higher-order terms. This introduces non-linear errors in the estimation of both the amplitude and the wavenumber correction vectors and requires iteration to obtain the true correction vectors. This is unnecessary in the linear filter estimation and we solve for the wavenumber correction vector directly without making this approximation and find that this is both more stable and more accurate when noise is present in the data. It has the added advantage that the amplitude and phase information are kept separate in the calculations, enabling the phase velocity to be evaluated independently of the amplitude calibration of the waveforms.

All of the earthquakes used in phase-velocity determinations (Table 2) lie within ∼ 5° of the great-circle path between the stations, so the effects of difference in source phase at the two stations and event mislocation are negligible. The accuracy of the dispersion curve could be affected by the short path length between Kula and Us,/e1>14;ak (90km); however, timing errors should be negligible in the digital data. Phase velocities were measured for only those portions of the waveform which showed simple dispersion. The squared coherency plot indicates the portions of the dispersion measurement which are well constrained by the data. The waveform fits (Fig. 3) between the far seismogram and the filtered near seismogram also give some indication of the accuracy of the phase term in eq.(1). Mindevalli & Mitchell (1989)) inverted the Rayleigh- and Love-wave group-velocity data to obtain crustal models for eastern and western Turkey. Our reason for measuring the Rayleigh-wave phase velocity was to provide additional velocity constraints in the receiver function analysis discussed below. Instead of inverting the phase-velocity data, we computed theoretical dispersion curves for simple models with the crustal thickness varying from 25 to 50 km, the crustal velocity varying from 3.1 to 3.9 kms^1-1 and the mantle velocity varying from 4.3 to 4.9kms^1-1 in an attempt to match the observed dispersion. Fitting the short-period dispersion required crustal velocities of 3.4 ± 0.1kms^1-1, whereas matching the long-period portion of the curves required 4.5 ± 0.2kms^1-1 half-space velocities. With these velocities, a 32km thick single-layer crust fits the Kula–Us,/e1>14;ak dispersion curve. The phase-velocity dispersion curve can be better matched by a model consisting of a two-layer crust, each layer 15km thick with a 3.3kms^1-1 upper layer, a 3.5kms^1-1 lower layer and a mantle consisting of a 4.4kms^1-1 layer 50km thick overlying a 4.5kms^1-1 half-space. We tested more complex models with a greater νmber of layers, but we found that the fits could not be significantly improved over the models discussed here. These models are less complex than those of Mindevalli & Mitchell (1989)) but have about the same average velocity structure. The models in Fig. 2 show that lower-crustal shear-wave velocities are low. The average velocity in the upper 15km of the crust is ∼ 3.3kms^1-1 ; the lower-crustal velocity is ∼ 3.5kms^1-1, These dispersion-derived estimates of the crustal shear-wave velocities are compatible with the S_g and L_g apparent velocity measurements ((Kadinsky-Cade et al. 1981; Saunders 1996). The upper-mantle structure consists of either a 4.5kms^1-1 half-space or a 50km thick low-velocity (∼ 4.4kms^1-1) layer overlying a 4.5kms^1-1 half-space. These upper-mantle shear-wave velocities are low for the observed P_n velocity, assuming a Poisson’s ratio of 0.25.

2.3 Teleseismic receiver functions

The use of receiver functions to determine crust andupper-mantle structure beneath three-component broad-band seismograph stations is now a well-established technique. We have determined true-amplitude radial and tangential receiver functions by computing the spectral ratio of the radial or tangential component spectrum and the vertical component spectrum (Langston 1979; Ammon 1991) for the events noted in Table 2. (The events in the following discussion are referred to by a νmber consisting of the year, month and day.) The deconvolution is stabilized by specifying a minimum spectral trough fill, c, for the vertical spectrum. The resulting receiver function is smoothed by a Gaussian filter whose width factor, a, limits the high-frequency content of the final waveform. Values for c were selected using a qualitative assessment of deconvolution stability and noise levels (Saunders 1996) and by considering the form of the averaging functions (the vertical component deconvolved from itself, with given c values). Values for a were chosen by considering noise levels in the data and a polarization analysis as a function of frequency to assess the bandwidth of scattered energy in the signal (Saunders 1996). All deconvolutions in this study were computed for c=0.01–0.001 and a=1.0. This Gaussian low-pass filter passes frequencies up to ∼ 0.4Hz, which is sufficient for resolving the main features of the crustal structure

When available, several high signal-to-noise ratio (SNR) receiver functions from a small range of backazimuths and epicentral distances were stacked to enhance the receiver function and reduce the noise level (Fig. 4). However, because of the short duration of the seismograph deployment at Kula and Us,/e1>14;ak, we did not have a large νmber of teleseisms for stacking. In several cases where there was a cluster of events, one event was clearly superior in terms of the SNR or averaging function. In such cases, we compared the receiver functions from all of the events in the cluster to assess the stable features of the receiver function waveform (Fig. 4) but then used the event with the higher SNR or superior averaging function for analysis.

We determined the crustal velocity structure beneath the stations using a combination of inversion and forward modelling. 1-D velocity models were determined using the time-domain linearized inversion procedure of Ammon, Randall & Zandt (1990). The starting model was parametrized as a stack of thin, horizontal layers to a depth of 60km. The S-wave velocity was the free parameter in the inversion, the layer thicknesses were fixed, the P-wave velocity was set assuming a Poisson’s ratio of 0.25, and the density was set using the relationship ρ=0.32V_p+ 0.77 (Berteussen 1977). The starting model was randomly perturbed into a large number of different starting models, and the radial receiver function was inverted by minimizing the difference between the observed receiver function and the synthetics computed for these models, while simultaneously constraining the model smoothness. This pseudo-Monte Carlo approach reduces the dependence of the inversion convergence on the form of the initial starting model. The inversion produces a range of solutions that fit the observed receiver function to differing °s. Those models which produced a good match to the data, given the noise levels present and waveform coherence, were selected for further study. As the layering became more apparent during the inversions we grouped adjacent thin layers with similar velocities into a single thicker layer and reinverted the receiver function.

We then used forward modelling to reduce the model complexity and assess how well individual model features were constrained by arrivals in the observed receiver functions. In the forward modelling we required that these simplified earth models contain the general features of the inversion models but with a small number of model parameters. We also required them to be consistent with the average crust and upper-mantle velocities measured from regional earthquakes and with the surface-wave dispersion. Large tangential motion in the receiver functions indicates laterally heterogeneous structures. Where large tangential motions are observed we attempt to match only the largest-amplitude radial arrivals in the forward modelling so as to minimize the effect of varying lateral structure on the 1-D earth model.

An important consideration in the interpretation of receiver function data is the lateral extent of the regions sampled by individual arrivals in the waveform, that is the degree of structural ‘averaging’ present in the data. This may be quantified, to an extent, by considering the size of the first Fresnel zone associated with a given arrival, illustrated schematically in Fig. 5. Paths from the edge of the first Fresnel zone differ in length by a half-wavelength from paths through the zone centre. Contributions from the wavefield propagating along all paths within the first Fresnel zone will interfere constructively to produce the observed signal at the receiver. The energy travelling from the zone centre follows a ‘geometric’ path, that is it obeys Snell’s law and Fermat’s principle. The radius of the Fresnel zone is defined as_R=√ (h+λ/2)2 −h² and is taken to represent approximately the region sampled by energy of wavelength λ at the interface at depth h. For the case of the Moho at 30km depth and receiver functions with a dominant period of 3–4s corresponding to a wavelength of λ≈ 14km, the radius of the first Fresnel zone is ∼ 32km. Therefore, stations separated by ∼ 10km laterally at the surface, as at the Kula and Us,ak stations, record P_s converted phases which sample much of the same structure at the Moho. This region can be defined for any point along the propagation path and for all phases which contribute to the receiver function. Thus, similar considerations apply to the Moho multiples PpP_ms and PpS_ms+PsP_ms except that these arrivals sample at even greater distances from the seismograph location. Therefore, the Moho depth determined from a teleseismic receiver function is not a point measurement but is an average over a considerable region (30–50km) surrounding the seismograph.

2.3.1 Analysis of the Kula data

KEN was the quietest of the Kula array stations, and the following analysis concentrates on data from it. Two SAL receiver functions were forward modelled but no data from SIH were deemed suitable for analysis. Several single-event KEN receiver functions were inverted, and two examples are shown in Fig. 6. The 930515 receiver function (Figs 6a–b) consists of the P_s (3.8s), PpP_ms (13s) and PpS_ms+PsP_ms (17.5s) Moho phases. Inversion of this receiver function resulted in two groups of models (A and B) whose synthetic receiver functions fit the observed receiver function equally well. These and all of the inversion receiver functions discussed below are the final inversion results after testing the effects of various layer parametrizations and smoothing. Group A has a ∼ 3.6kms^1-1 upper crust overlying a ∼ 3.3kms^1-1 lower crust and a 32–34km deep Moho; group B has a ∼ 3.8kms^1-1 upper crust overlying a ∼ 3.4kms^1-1 lower crust and a 33–35km deep Moho. The 940314 KEN receiver function (Fig.6c) has a larger tangential component than 930515, but the radial is similar except that the PpP_ms phase is advanced by about 1s. Inversion of the 940314 receiver function results in a set of crustal models similar to model group A of 930515.

Inversion model group A of 930515 and the inversion results of 940314 are similar, and for this reason inversion model group A of 930515 was taken as the starting model for the KEN forward modelling test. The 930515 receiver function is well matched by the synthetic receiver function for a simple crustal model (Fig.7a) consisting of a 2km thick 2.45kms^1-1 surface layer, a 30km thick 3.5kms^1-1 crust and a 4.5kms^1-1 mantle (Model K; Fig.7a). Inclusion of the low-velocity surface layer produces the observed delay in the direct arrival and enhancement of the negative motion at 2s. The inclusion of the thin layer above the Moho improves the fit of the delay and amplitude of the positive and negative reverberations at 13 and 17.5s, respectively. The 4.5kms^1-1 sub-Moho velocity is low considering the observed 8.1 ± 0.1kms^1-1P_n velocity. If the sub-Moho velocity is increased to 4.7kms^1-1 (consistent with a Poisson’s ratio of 0.25 and a P_n velocity of ∼ 8.1kms^1-1), the amplitude of the P_s phase is fitted better but the PpP_ms and PpS_ms+PsP_ms amplitudes are then too large (Fig.7a). A 4.5kms^1-1 half-space velocity matches the surface-wave dispersion (see below), whereas a 4.7kms^1-1 sub-Moho velocity requires an upper-mantle velocity reversal to match the surface-wave dispersion. Inversion model group A of 930515 (Fig.6a) suggests a velocity reversal in the mid-crust, and including a thin mid-crustal low-velocity zone (LVZ) (Fig.7b) in the forward model improves the fit of the phase at ∼ 8s. The inclusion of a mid-crustal discontinuity rather than the LVZ (Fig.7c) also improves the fit of the ∼ 8s phase but further degrades the fit of the P_s phase. In this instance and in most of the following cases we favour the model with amid-crustal LVZ, although we must acknowledge that there is only weak evidence for choosing this over the model with a mid-crustal discontinuity.

In modelling the other KEN receiver functions (Figs 8 and 9), Model K was taken as the starting model, then minor adjustments were made in the layer thicknesses and velocities to maximize the fit of the synthetic waveform to the observed waveform. It is possible to fit the majority of the KEN receiver functions with only small changes to Model K. The 930525 noise level (Fig.8a) is higher than either 930515 or 940314, but the same main arrivals are clearly visible and are matched by the synthetic arrivals from Model K. However, the inclusion of a thin LVZ results in a better match of the waveform at ∼ 8s. It is difficult to match the timing of both the PpP_ms and the PpS_ms+PsP_ms phases arriving at 13 and 17.5s respectively, which may indicate that there is more structure to the Moho or that lateral inhomogeneities exist within the crust. No low-velocity surface layer is necessary in matching the 930809b receiver function (Fig.8b), but the prominent arrival at ∼ 8s suggests a 3.25kms^1-1 LVZ between 15 and 19km depth. The timing of the arrival at ∼ 13s is not well matched, and no simple structures were found which simultaneously matched this and the subsequent ≠gative arrival at ∼ 16s. The amplitudes of the P_s phase at 3.5s and the following negative arrival are underestimated by this model. There is also substantial unmatched energy on the radial component at times greater than ∼ 20s. The 930807 receiver function (Fig.9a) contains noise and large tangential motion during the first 10s; however, the main conversions and reverberation from the Moho can still be identified. A mid-crustal LVZ is not required by these data since there is no evidence of significant energy arriving between the P_s (4 s) and PpP_ms (14s) phases. Fig. 9(b) shows similar results for 931113. The speed of the surface layer is slightly faster than for Model K, and a deeper mid-crustal LVZ (15.5–19.5km depth) is required to fit the timing of the positive arrival at ∼ 9s.

The SAL 930918 receiver function (Figs 10a and b) is similar to those at KEN except that the tangential amplitudes are greater, the positive phase at ∼ 8s is stronger and the PpS_ms+PsP_msphase is less prominent. There is also a small delay of the direct and P_s arrivals relative to KEN, indicating a thicker low-velocity surface layer. The 930809b SAL receiver function is almost identical to 930918 (Fig. 4), but we have modelled only the 930918 receiver function because it has lower noise at long periods. A simple crustal section (Fig.10a) similar to Model K, with constant velocities between depths of 2 and 29.5km, can match the Moho converted phases and reverberations but is unable to match the ∼ 8s arrival. This phase is better matched by including a mid-crustal LVZ at 13–17km depth. This phase could instead result from a velocity step in the mid-crust (Fig.10b). Although a small step contrast from 3.2 to 3.5kms^1-1 at a depth of 17km generates a positive arrival at ∼ 8s, the amplitude of this arrival and the negative arrivals throughout are then poorly matched and the fit to the direct arrival is degraded. By increasing the contrast across the mid-crustal step (Fig.10b), the amplitude of the 8s arrival is fitted better, but the negative arrivals and the P_s phase amplitude are still not well matched. It therefore seems that both the positive and the negative contrasts associated with an LVZ are necessary in fitting the data. None of the models produces a good fit to the signal beyond ∼ 17s.

The SAL 930929 receiver function (Fig.10c) is similar to the 930918 receiver function except that it contains significant energy on both the radial and the tangential components at times greater than 20s. The model shown is the same as that for the mid-crustal LVZ for 930918 (Fig.10a) except that the lower-crustal section is 0.5km thinner. The waveform is well matched except for some minor misfit in the amplitude of the direct P and P_s phases and the mismatch of the arrivals at times greater than 15s. The amplitude of the arrival at 8s, which is modelled as a reverberation from the mid-crustal LVZ, is slightly overestimated. This may indicate that the velocity of the material between depths of 13 and 18km is slightly higher than in the model, but given the level of tangential energy arriving at about this time, this is not necessarily the case.

The majority of the Kula array receiver functions can be matched with simple crustal models (see Fig. 17 below) consisting of a 2km thick low-velocity surface layer, a ∼ 15km thick upper crust with shear-wave velocity in the range3.0–3.5kms^1-1, a ∼ 12km thick lower crust with shear-wave velocity in the range 3.5–3.7kms^1-1, and a crustal thickness value in the range 29–31 km. Most waveforms considered are better matched by structures containing a ∼ 4–5km thick LVZ at depths between 13 and 20km in which the shear velocities drop to values of ∼ 3.25kms^1-1, Converted phases and reverberations from the LVZ fit the positive arrival consistently observed between the P_s and PpP_ms arrivals and the negative arrivals on either side of P_s. Receiver functions provide little constraint on the upper-mantle shear velocity, except that it should lie in the range 4.5 ± 0.2kms^1-1, This upper-mantle shear-wave velocity also agrees with the surface-wave dispersion data. The increased delay in the direct arrival at SAL compared with KEN can be modelled by introducing a slightly thicker or slightly lower-velocity surface layer. This apparent delay of the direct arrival results from interference of the immediately following P_s conversion from the sediment–basement interface. It seems unlikely that the relative P_sdelay between KEN and SAL is related to different crustal thicknesses beneath the stations given the significant overlap of the sampling regions of this phase at the Moho of the two stations (Fig. 5).

2.3.2 Analysis of the Kula data

The discussion in this section concentrates on data from station DAN, the quietest of the Us,/e1>14;ak array stations. The DAN results are supported by the result from the one GUS receiver function modelled. DAN receiver functions for the stack of the three Hindu Kush events (930809a, 930809b and 930918; Fig.11a) and a single Indian event (Fig.11b) were inverted. The Hindu Kush receiver function at DAN consists of the P_s (∼ 4s), PpP_ms (∼ 14s) and PpS_ms+PsP_ms (∼ 17s) Moho phases. The inversion models (Fig.11a) have low average crustal velocities (3.3–3.5kms^1-1) and a gradational Moho between 30 and 35km depth. The 930929 receiver function (Fig.11b and c) is similar to the Hindu Kush receiver function except for the presence of the positive arrival at ∼ 8s. Its inversion produces two groups of models whose synthetic receiver functions satis-factorily fit the data. The velocity profile of the first model group (Fig.11b) is similar to the Hindu Kush model, while the second set (Fig.11c) has a much higher average crustal velocity (3.6–3.7kms^1-1), a zone of high-velocity material (∼ 3.8kms^1-1) from 8 to 16km depth, an LVZ in the range28–32km and a deeper Moho (32–35 km).

The Hindu Kush inversion model was taken as the starting velocity model in the forward modelling. Fig. 12(a) shows a simplified velocity model consisting of a three-layered crust with an average velocity ∼ 3.4kms^1-1 and a gradational Moho between depths of 29.5 and 33.5km (Model U; Fig.12a) whose receiver functions match the DAN Hindu Kush receiver function. The other Us,/e1>14;ak crustal models are all modifications of Model U. This model does not match the arrival time of the negative phase at 17s, but there is some tangential motion associated with this arrival which may be distorting the wave shape. No simple changes in the crustal structure were found which caused the PpS_ms+PsP_msphase to arrive any earlier while leaving the times of the P_sand PpP_msunchanged. The model (Fig.12a) which fits the DAN–Hindu Kush receiver function better has the same crustal structure as Model U, but has a 4.6kms^1-1 layer between 33.5 and 41.5km depth, a 4.2kms^1-1 LVZ from 41.5 to 46.5km depth and a 4.5kms^1-1 half-space.

The 930918 GUS receiver function (Fig.12b) is matched by the synthetic from a model identical to Model U except that the crust is 0.5km thinner and the speed of the surface layer is 0.3kms^1-1 slower. The major phases are matched, including the ≠gative PpS_ms+PsP_ms arrivals at ∼ 18s. The amplitude of the direct arrival is underestimated and the time of the P_s phase is early, but there is some tangential motion in phase with P_s and it is difficult to assess the stability of these features with a single waveform. However, it is clear that the major features in the waveform can be matched without requiring any significant variations in upper-mantle structure. Given the overlap of the sampled regions beneath DAN and GUS, it is likely that the differences between the PpS_ms+PsP_msarrivals result from structural differences in the shallow crust.

Fig. 13 shows the forward-modelling solution for events 940205 and 930929 recorded at DAN. The 940205 data are noisy, but the major features of the receiver function are matched by a simple model having an upper-crustal velocity slightly lower than Model U and a lower-crustal velocity slightly higher. The small positive arrival at ∼ 8s is matched by a reverberation from the mid-crustal discontinuity. The fit to the negative arrivals in the first 10s of this waveform is not good, but this is acceptable given the extremely large side-lobes (17 per cent of the direct arrival amplitude) present in the averaging function for this event. The 930929 receiver function is matched by synthetic receiver functions from an almost identical structure except for a lower-velocity surface layer. The 930929 receiver function is fitted equally well by changing the mid-crustal discontinuity into a mid-crustal LVZ. The negative arrival at ∼ 17s is difficult to match with the PpS_ms+PsP_msreverberation, but given the results of the forward modelling of the Hindu Kush receiver functions, the misfit is not significant.

The majority of the Us,/e1>14;ak receiver function data can be matched with simple crustal models consisting of a 2km thick low-velocity surface layer, a ∼ 17km thick upper crust with a shear-wave velocity in the 3.3–3.4kms^1-1 range and a ∼ 14km thick lower crust with a shear-wave velocity in the 3.5–3.6kms^1-1 range (see Fig. 17 below). There is less variation between the Us,/e1>14;ak crustal models than among the Kula crustal models. Most receiver functions at Us,/e1>14;ak do not require a mid-crustal velocity discontinuity. The Moho is more gradational and deeper beneath the Us,/e1>14;ak array than beneath the Kula array. An upper-mantle velocity of ∼ 4.5kms^1-1 matches both the majority of the receiver function data and the surface-wave dispersion data (see below).

2.3.3 Analysis of the ANTO data

The ANTO receiver functions were derived from events larger than those recorded at the Kula and Us,/e1>14;ak seismographs, and this generally resulted in a higher SNR. Small apparent delays in the direct radial arrival of the receiver functions at ANTO indicate the presence of low-velocity, near-surface material. Fig. 14(a) shows the inversion results for a stack of waveforms from events 930904 and 941025 from the Hindu Kush. The majority of features in the waveform can be matched with a structure which has a low-velocity (∼ 2.5–3kms^1-1) 4km thick surface layer, a nearly uniform (3.5–3.6kms^1-1) crustal section between 4 and 34km depth, and a gradational Moho between 34 and 38km depth. Below ∼ 40km, there is a LVZ where velocities drop to ∼ 4.2kms^1-1, The P_s phase arrives at ∼ 5s, and the upward trend in the radial data at around 15s may represent the PpP_ms reverberation from the Moho. Beyond times of ∼ 16s, however, the nature of the radial waveform becomes more complex, and it is difficult to identify clearly the PpS_ms+PsP_ms reverberations. The negative arrival at 7.5s is poorly matched.

Fig. 14(b) shows inversion results for a stack of Japanese Islands and Kuril Islands events (930115, 940828, 941004, 941009, 941016 and 941228). The inversion model consists of a ∼ 2.6–3.1kms^1-1, 2km thick surface layer and an oscillatory crustal section with average velocity 3.5–3.6kms^1-1, The spread in mid-crustal velocities is fairly large. It was not possible to simplify further the inversion models by combining layers and then reinverting the data as we did in the other inversions. There is a prominent velocity step from about 3.5 to 3.7kms^1-1 at 12km depth, a gradational Moho between 34 and 38km depth, a LVZ (4.2–4.3kms^1-1) from 46 to 54km depth and a 4.5–4.7kms^1-1 half-space.

The 940802 Kamchatka receiver function (Figs 14c and d) is similar to the Japan/Kuril receiver function except that the large positive arrival seen at 9s in the Japan–Kuril receiver function is absent. Synthetics from two groups of inversion models match the 940802 receiver function. The first model group (Fig.14c) fits the main features of the receiver function with a simple model consisting of a 2km thick ∼ 2.5kms^1-1 surface layer, an 8km thick 3.45kms^1-1 upper crustal layer, a 24km thick ∼ 3.65kms^1-1 lower crustal layer, a gradational Moho between 34 and 38km depth and a 4.5–4.6kms^1-1 upper mantle. This simple model produces a good match to the P_s (5s), PpP_ms (15s) and PpS_ms+PsP_ms (19s) arrivals. The amplitudes of the P_sphase and the immediately following negative arrival are slightly underestimated. However, this later negative arrival is in phase with large-amplitude tangential motion. Crustal velocities of the second model group (Fig.14d) are in general ∼ 0.2kms^1-1 faster than those of the former model group, resulting in a deeper Moho (36–40km depth).

The inversion models for the Hindu Kush (Fig.14a), Japan/Kuril (Fig.14b) and the first model group for 940802 (Fig.14c) are similar, so we have used these as the starting model in the forward modelling. The Japan/Kuril receiver function (Fig.15a) is fitted by synthetic receiver functions from a crustal model consisting of a 2km thick 2.4kms^1-1 surface layer, an 8km thick 3.2–3.3kms^1-1 upper-crustal layer, a 23km thick 3.6kms^1-1 lower-crustal layer, a transitional Moho between 34 and 37km depth and a 4.5kms^1-1 upper-mantle half-space. However, increasing the crustal thickness by 0.5km and including a mantle LVZ between 45.5 and 53.5km depth results in a better fit of the receiver function between 18 and 22s (Model A; Fig.15a). The other ANTO forward models are derived from Model A. The forward model obtained for the Hindu Kush receiver function consists of a 2.4kms^1-1 2km thick surface layer, a 10.5km thick 3.2–3.3kms^1-1 layer, a 20km thick 3.6kms^1-1 layer and a 4.5kms^1-1 upper-mantle half-space. Including a 4.2kms^1-1 LVZ between 40.5 and 44.5km depth slightly improves the fit between 18 and 22s in the receiver function, although the improvement is minor.

The receiver function for event 940818 (Fig.16a) constrains the structure to the south of ANTO. The fit achieved by both models shown is good given the level of noise and tangential motion. The model with a slightly thinner crustal section and deeper mid-crustal velocity step, however, produces a better match to the negative arrivals at 8 and 19s. The 940802 and 930515 receiver functions are also fit by the same model as that derived from the Japan–Kuril receiver function (Figs 16b and c).

The majority of the ANTO receiver functions can be matched with a simple crustal model (see Fig. 17), despite the large tangential amplitudes observed at this station. The average crustal structures consist of a 2km thick 2.4kms^1-1 surface layer, an 8 to 10km thick ∼ 3.3kms^1-1 upper crust and a 20–22km thick ∼ 3.6kms^1-1 lower crust. There is little variation between the ANTO crustal models. The P_s amplitude and dispersion data can both be matched with a 4.5kms^1-1 upper-mantle half-space. The receiver function arrivals constraining the upper-mantle LVZ arrive late in the waveform (>20s). This part of the waveform is also likely to be more affected by scattered energy and for this reason we do not feel that the upper-mantle LVZ is required by the data. The velocity structure at ANTO is therefore similar to that at the Kula and Us,/e1>14;ak arrays, except for the thicker crustal section and the slightly higher average crustal velocities.

3 Discussion And Conclusions

In this study we have investigated the crustal structure beneath two sites in the extending region of western Turkey and compared these to the crustal structure beneath a site in the central Anatolian Plateau. The crustal models were determined primarily by analysing teleseismic receiver functions, but apparent velocities of local and regional seismic phases and surface-wave phase- and group-velocity dispersion pro-vided additional constraints on the crust and upper-mantle velocities. We determined 1-D velocity models from the receiver function data, although lateral variations in structure are clearly present. We first inverted the receiver function data at each station for those azimuths least affected by off-azimuth arrivals. We then simplified the models and tested their features by forward modelling to reveal the main crustal features while minimizing the effects of the lateral variations in structure Although the tangential amplitudes are relatively large in the three sets of receiver functions, there is no indication from these of systematic tangential arrivals diagnostic of dipping interfaces or crustal anisotropy as discussed by Peng & Humphreys (1997).

The results of our receiver function analysis are summarized in Fig. 17. The best-constrained feature of the crustal models is the depth to the Moho, which clearly increases from beneath Kula (∼ 30km), in the west within the extending crust, to Us,/e1>14;ak (∼ 34km), approximately on the edge of the extending zone, to Ankara (∼ 37.5km), in the unextended Anatolian Plateau. This is also apparent in the receiver function data themselves (Fig. 18), especially for the first Moho conversion P_s, but also to a lesser extent for the Moho multiples PpP_ms and PpS_ms+PsP_ms. For all of the receiver function data the differential time between direct P and P_s increases from 3.67 ± 0.12s at Kula, to 4.07 ± 0.12s at Us/e1>14;ak, to 4.68 ± 0.12s at Ankara. The differential time between direct P and the Moho multiples, although having more variation, shows the same pattern [PpP_ms: Ankara (14.87 ± 0.24), Us,/e1>14;ak (13.70 ± 0.11) and Kula (13.07 ± 0.14); PpS_ms+PsP_ms: Ankara (19.34 ± 0.18), Us,/e1>14;ak (16.95 ± 0.60) and Kula (16.91 ± 0.57)].

The average crustal shear-wave velocities beneath the three sites are similar. The lower crust at Kula and Us,/e1>14;ak exhibits low velocities (∼ 3.5kms^1-1) for continental crust. Moreover, velocities much higher than this result in a poor fit of the short-period (<30s) surface-wave phase velocity (Fig. 19). These lower velocities suggest the presence of silica-rich material in the lower crust beneath western Turkey. If this is the case, this also indicates that significant underplating of the Moho by basalts has not occurred in western Turkey. This range of velocities is consistent with low-velocity quartz schists and granites, although typical values for such rock types are slightly higher at lower crustal pressures (Carmichael 1982). Since both the crustal stretching and the presence of volcanics suggest that moderately high heat flow is likely in western Turkey, it is possible that elevated lower-crustal temperatures may have reduced the velocity of material in the lower crust. The variation in gravity across central and western Turkey is small (Lemoine et al. 1997). Converting each of the crustal sections in Fig. 17 into a density profile assuming a Poisson’s ratio of 0.25 and the V_p–density relationship of Berteussen (1977), and balancing these density profiles isostatically against the structure of a mid-ocean ridge (Cochran 1982) predicts an excess elevation due to the thickened crust of 1630m at Ankara, 1290m at Us,/e1>14;ak and 780m at Kula, all of which are close to the observed elevations.

If the crustal section beneath Ankara represents crust unstretched by the current phase of extension and the crust beneath Kula and Us,/e1>14;ak differs only in that it has been stretched, then the β -factors (=original thickness/new thickness) are ∼ 1.27 in the vicinity of Kula and ∼ 1.12 around Us,/e1>14;ak. These values agree with the geologically derived β -factor of ∼ 1.25 estimated by Paton (1992a, b) from surface faulting observations in western Turkey, with the β -factor estimated from integrating seismic strain rates of ∼ 1.3 (Eyidoğan 1988) in western Turkey, and also support the general observation of increasing extension westwards. The deformation within the crust resulting from the extension is approximately uniformly distributed with depth, and the discrete faulting observed at the surface in western Turkey provides a good estimate of the amount of extension throughout the crust. These data say nothing about the level of extension in the lithospheric mantle.

Dispersion curves for the Us,/e1>14;ak and Ankara velocity models match the Kula–Us,/e1>14;ak dispersion data better than does the dispersion curve for the Kula model (Fig. 19). This may indicate that the structure beneath Us,/e1>14;ak is more representative of the path-averaged structure between Kula and Us,/e1>14;ak, and that the greater crustal thinning beneath Kula compared to Us,/e1>14;ak does not extend very far east of Kula. This is also indicated by the faulting, seismicity and topography (Fig. 20). Most of the large normal faults occur from the region of Kula westwards. This is also apparent in the mean elevations estimated to the nearest 100m from the ETOPO5 topo-graphic data set. Between Ankara (∼ 1400m elevation) and Us/e1>14;ak (∼ 1100m elevation), about 300km distance, the mean elevation decreases by ∼ 300m; between Us,/e1>14;ak and Kula (∼ 700m elevation), a distance of about 75km, the mean elevation decreases by ∼ 400m. In the 125km to the west of Kula the elevation drops to sea level, and there is an additional drop of ∼ 500m elevation between the coast and the central Aegean about 150km further west. This variation in elevation suggests that most of the extension has occurred west of Kula. ) and ) find that the crustal thickness beneath the central Aegean sea is about 22 km, which also suggests that most of the extension-related crustal thinning occurs west of Kula.

The sub-Moho shear velocity is ∼ 4.5kms^1-1 beneath both central and western Turkey and allows the receiver function and long-period surface-wave dispersion data to be simultaneously matched with simple mantle models. The apparent P_n velocity is 8.1 ± 0.01kms^1-1 in western Turkey and 7.8 ± 0.01kms^1-1 in central Turkey, requiring an upper-mantle Poisson’s ratio of ∼ 0.25 beneath Ankara and ∼ 0.28 beneath Kula and Us,/e1>14;ak. The presence of strong P_n arrivals on the regional seismograms indicates that a positive P-wave velocity gradient exists immediately beneath the Moho. The absence of a high-frequency regional S_n phase indicates either a negative sub-Moho shear-velocity gradient, and hence an increasing Poisson’s ratio with depth in the upper mantle, or a zone of significant shear-wave attenuation immediately below the Moho. An upper-mantle shear-wave velocity of 4.5kms^1-1 in western Turkey is similar to that observed in the Basin and Range province of western North America, a tectonically similar region. However, it is not clear from the phase-velocity data whether the upper mantle of the two regions is similar. The short-period (< 30s) phase velocity in western Turkey is identical to that in the Basin and Range (); however, the long-period (>60s) phase velocity is ∼ 0.2kms^1-1 faster in western Turkey than in the Basin and Range, which has a substantial upper-mantle low-shear-wave-velocity zone below the upper-mantle lid. No such LVZ is required to fit the western Turkey long-period phase-velocity curve. However, if the sub-Moho shear-wave velocity is taken to be 4.7kms^1-1, then an upper-mantle shear-wave LVZ is required to fit the observed dispersion. ) found the sub-Moho shear-wave velocity east of Ankara to be about 4.30 ± 0.05kms^1-1 with a positive upper-mantle gradient. In western Turkey they found the sub-Moho shear velocity to be as low as 4.20kms^1-1 for the inversion of their Rayleigh-wave dispersion data, and as high as ∼ 4.75kms^1-1 for the inversion of their Love-wave dispersion data. The joint Rayleigh- and Love-wave dispersion inversion indicates a 4.60kms^1-1 sub-Moho velocity. In each case there is a positive shear-wave gradient in the upper mantle. They suggest that the difference indicated by the Rayleigh- and Love-wave data for western Turkey may result from crustal aniso-tropy. Our surface-wave dispersion data do not address this question; however, crustal anisotropy is not required by the receiver function data. The lack of an S_n phase calls into question Mindevalli & Mitchell’s (1989) positive upper-mantle shear-wave-velocity gradients.

Two other features of Fig. 17 are worth noting even though they are less well constrained by the receiver function analysis. First, the Moho appears to be fairly sharp beneath Kula in the extending region, and gradational over about 6km depth beneath the unextending region ≠ar Ankara. This is also suggested by the relative amplitudes of the P_sphases in Fig. 18. The second point is that the spread of crustal velocities between the models at each site becomes less in going from Kula to Ankara. This may be the result of more pronounced crustal deformation in the region of greater extension, leading to more scattered energy in the P-wave coda at Kula. Ankara is in the Cretaceous suture zone in central Turkey where the crust is severely deformed (Şengör & Yilmaz 1981). However, the reduced spread in the Ankara models may result from higher SNR due to the larger events used to determine the ANTO receiver function. This explanation does not, however, apply to the Kula and Us,/e1>14;ak data.

One intriguing feature of these receiver function data is the presence of the arrival consistently seen between the P_sand PpP_msphases on many of the radial receiver functions, primarily at the Kula array. This arrival is small compared to the main Moho arrivals and may result from scattered energy in the P-wave coda resulting from lateral variations in the crustal structure in the vicinity of the seismograph site. However, the consistency of this arrival at Kula, and to a lesser extent at Us,/e1>14;ak, has led us to consider the possibility that it may result from a mid-crustal discontinuity or a mid-crustal LVZ. This arrival is best matched by including a mid-crustal low-velocity layer with a ∼ 7 per cent velocity reduction(3.24–3.35kms^1-1) in the model.

A crustal LVZ may result from compositional changes, high temperatures, fluids at high pore pressure, or the presence of partial melt. The heat flow in western Turkey is unknown, but based on the volume of the volcanics (Bunbury 1992) and the amount of stretching (Eyidoğan 1988; Paton 1992a, b; this study), it is unlikely to be high enough for ∂ melt in the mid-crust (McKenzie 1978b). The abundant alkali basalt extrusives in the vicinity of Kula suggest that any crustal LVZ may be related to the extensional volcanism, but the volume of extruded material is small (∼ 2.3 km³) and the volcanics show little evidence of crustal contamination (Paton 1992b). It is therefore unlikely that sufficient quantities of this material exist in the crustal section to be detectable at the wavelengths considered in the receiver function analysis. The simplest explanation for the possible crustal LVZ is the existence of a compositional change in the mid-crust. Mueller (1980) observe a LVZ (∼ 3.12kms^1-1) at depths between ∼ 13 and 18km along a traverse between the Rhine graben and the Po Plain, which they equate with a granitic intrusion. It is also possible that the mid-crustal LVZ, if it exists, may be related to an earlier phase of andesitic volcanism in western Turkey (Şengör & Yilmaz 1981).

Acknowledgements

We would like to acknowledge the invaluable logistical assistance of Celal Şengör for getting the Kula and Us/e1>14;ak seismographs sited and operating. Judith Bunbury helped in the site selection, James McKenzie built the random binary calibration and seismometer self-centring systems, Cansun Güralp helped with various aspects of the instrumentation, Brek Betts helped to fix Fig. 5 and Dave Sandwell helped to produce Fig. 20. We would also like to thank Dan McKenzie and James Jackson for numerous enlightening discussions and for critically reading the manuscript. This research was supported by NERC grant GR3/7843. This is Cambridge University Department of Earth Sciences contribution 5171.

References

Author notes