# Constraints on absolute S velocities beneath the Aegean Sea from surface wave analysis.

E. Bourova, I. Kassaras, Helle A. Pedersen, T. Yanovskaya, Denis Hatzfeld, A. Kiratzi

► To cite this version:

E. Bourova, I. Kassaras, Helle A. Pedersen, T. Yanovskaya, Denis Hatzfeld, et al. Constraints on absolute S velocities beneath the Aegean Sea from surface wave analysis. Geophysical Journal International, 2005, 160 (3), pp.1006-1019. 10.1111/j.1365-246X.2005.02565.x. hal-00109942

HAL Id: hal-00109942

https://hal.science/hal-00109942v1

Submitted on 22 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Constraints on absolute S velocities beneath the Aegean Sea from surface wave analysis

E. Bourova,<fnref n="1" />,<fnref n="2" /> I. Kassaras,<fnref n="3" /> H. A. Pedersen,<fnref n="1" /> T. Yanovskaya,<fnref n="2" /> D. Hatzfeld<fnref n="1" /> and A. Kiratzi<fnref n="4" />

Accepted 2004 December 15. Received 2004 December 13; in original form 2003 September 17

## SUMMARY

The 3-D structure of the lithosphere beneath the Aegean Sea is investigated through surface wave dispersion analysis. Rayleigh and Love waves recorded by 12 broad-band stations installed for a duration of 6 months in the Aegean region are processed through array analysis and Wiener filtering. Data from three GEOFON stations in the area of Crete were also used. The resulting two-station phase velocities are used to determine lateral variations of Rayleigh wave phase velocities between periods of 20 and 100 s by a 2-D ray tomography method. The obtained phase velocities are inverted to calculate variation of S-wave velocity with depth using a combination of linearized inversion and a Monte Carlo based non-linear inversion.

The absolute S-wave velocity is resolved to a depth of approximately 200 km. A high-velocity anomaly of 3 per cent is observed in the southern Aegean attributed to the Hellenic subduction. In the northern part of the Aegean, in the prolongation of the North Anatolian Fault which is influenced by strong extensional movements, we found low absolute S-wave velocities at 50–100 km depth. This supports a model of a distributed deformation of the upper mantle in the area. Separate Rayleigh and Love wave phase velocity inversions along common profiles reveal a strong Love–Rayleigh discrepancy in the northern Aegean down to at least 150 km depth, i.e. most probably including the top of the asthenosphere.

**Key words:** Aegean Sea, anisotropy, phase velocity, surface wave dispersion.

## 1 INTRODUCTION

It is still largely unknown whether continental lithospheric deformation is localized along major faults which extend through the whole lithosphere or is distributed over large areas (England & Jackson 1989; Thatcher 1995). In the latter hypothesis major faults are the surface expression of a distributed underlying deformation. Further complexity arises from our poor knowledge of the interaction between the lithosphere and asthenosphere, and the implication of the latter in plate tectonic driving forces (e.g. Bourne et al. 1998; Bokelmann 2002).

These questions are particularly difficult to answer because most of the available information related to the deformation is located at or close to the surface (GPS geodesy, tectonics, seismicity) and very little is related to parameters associated with the displacement or the deformation of the uppermost mantle.

The Aegean region, which is continental and lies between the African and Eurasian lithospheric plates, is particularly interesting in this context. First, the kinematics are well documented by geodetic measurements (McClusky et al. 2000), seismic focal mechanisms

(Papazachos & Kiratzi 1996; Jackson & McKenzie 1988), and tectonic observations for the Miocene time (Gautier & Brun 1994) and the Quaternary (Mercier et al. 1987; Armijo et al. 1996). Secondly, the lateral movements are very large. Even though the two plates converge at a rate of 1 cm yr⁻<fnref n="1" />, the convergence across the Hellenic Trench is approximately 3.5 cm yr⁻<fnref n="1" />. This high rate is due both to the westward motion of Anatolia relative to Africa and to the rapid and intense deformation that affects the region. The northern part of the Aegean, into which the North Anatolian Fault extends, is submitted to strong extensional strain.

The crustal structure and thickness of the Aegean is reasonably well known from a wide variety of studies using deep seismic soundings (Makris 1978; Makris & Stobbe 1984), gravity data (Tsokaz & Hansen 1997), surface wave data (Papazachos 1969; Kalogeras & Burton 1996; Karagianni et al. 2002) and receiver function inversion (Sodoudi et al. 2004), but we have a relatively poor idea of the mantle structure, which is a key to the understanding of the dynamics of the area. The mantle structure has been studied mostly by P-wave traveltime inversion using local and teleseismic events (Spakman et al. 1988; Papazachos et al. 1995; Papazachos & Nolet

<footnote n="1">¹Laboratoire de Geophysique Interne et Tectonophysique, Observatoire de Grenoble, France. E-mail: helle.pedersen@obs.ujf-grenoble.fr</footnote>

<footnote n="2">²Institute of Physics, Saint Petersburg State University, Russia</footnote>

<footnote n="3">³Seismological Laboratory, Department of Geophysics and Geothermics, University of Athens, Greece</footnote>

<footnote n="4">⁴Geophysical Laboratory, Aristotle University of Thessaloniki, Greece</footnote>
1997; Tiberi *et al.* 2000), the latter providing velocity perturbations relative to an input velocity model. Absolute velocities can be better estimated by traveltime inversion from local events, but the resolution of the data is poor for the deeper part of the lithosphere and the underlying asthenosphere due to the very limited intermediate depth seismicity (with a maximum depth of 180 km) in the Aegean region. This is particularly true for the northern Aegean where there is no resolution below 100 km depth. The present study supplements previous tomography results from the area (Spakman *et al.* 1988; Papazachos *et al.* 1995; Papazachos & Nolet 1997; Tiberi *et al.* 2000).

We here show results from a surface wave study using the best broadband station coverage of the region to date. The advantage of the surface waves, in spite of a poorer lateral resolution than in teleseismic body wave tomography, is the good constraints on absolute shear wave velocities, which are difficult to obtain with other tomography techniques. We use a teleseismic surface wave tomography for 3-D imaging of the lithosphere using a method proposed by Ditmar & Yanovskaya (1987) and Yanovskaya & Ditmar (1990). Phase velocities are calculated for periods of up to 100 s thereby giving information at depths down to approximately 200 km. Due to the short duration of the field experiment, we increased the number of useable events by combining two-station measurements with array analysis (e.g. Baumont *et al.* 2002).

## 2 DATA AND PHASE VELOCITY MEASUREMENTS

The data were obtained during a 6-month temporary seismological experiment in 1997 (Hatzfeld *et al.* 2001). The temporary network was composed of 30 digital three-component seismological stations covering the Aegean. For this study we used the 12 broad-band stations equipped with CMG3-T (100 s) and CMG40 (60 s) seismometers which recorded continuously at a sample frequency of 62.5 Hz. Data from three GEOFON stations equipped with STS2 (100 s) seismometers are also included. The distribution of the stations is presented in Fig. 1. Technical details are described in Hatzfeld *et al.* (2001).

Phase velocity dispersion curves for Rayleigh and Love waves can be determined between selected pairs of stations located on the same great circle as the epicentre, within a few degrees. The phase velocity analysis encounters three difficulties. First, the interstation distance at long periods is smaller than or approximately equal to the wavelength. In this case the phase differences between the two stations are small, and therefore even low noise can induce errors in large phase velocity. Second, the seismic rays are deflected from the great circle path due to variations in lateral velocity. Observed ray deflections are as high as 30° (Cotte *et al.* 2000) and result in incorrect phase velocity estimations as the propagation distance

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The main features of the region under study and the distribution of three-component broad-band seismological stations used in the surface wave tomography. The temporary network is composed of 12 broad-band stations equipped with CMG3-T (100 s) and CMG40 (60 s) and three GEOFON stations equipped with STS2 (100 s) sensors.</figcaption></figure>
between the two stations is different from the one inferred from the great circle path. This problem is aggravated with increasing angular difference between the station-station and the epicentre-station great circles. Third, the surface wave propagation velocity between two stations may not be the same as the velocity that corresponds to the phase velocity of the underlying medium (Wielandt 1993; Friederich et al. 1994). This problem can be at least partly resolved by avoiding the use of seismic events in the same epicentral area. However, the use of a variety of epicentral locations is in contradiction with the constraints due to the ray deflection, which imposes that the events for a given station pair are chosen within a small corridor around the great circle.

One solution to these conflicting priorities is to use an array tomographic approach allowing for non-plane wave fronts (Friederich 1998; Bruneton et al. 2002). However, such tomographic methods have a very large number of parameters to invert for, due to the need to describe the wave front of each event. We did not consider such an approach justified as only (at best) 12 seismic stations were operating simultaneously, and over a limited amount of time.

To accommodate these problems we applied a data processing method which is based on two-station measurements but which includes array measurements to make first-order corrections for frequency-dependent deviations from the great circle (Cotte et al. 2000; Baumont et al. 2002). This approach makes it possible to use events with angles between epicentre-station and station-station great circles of up to 25°, thereby increasing the amount of data and making it possible to use a wider variety of epicentre locations. This procedure increases the amount of data sufficiently that the average dispersion curve should also be relatively independent of random noise.

The array analysis is carried out by inverting the traveltime differences between the stations to obtain the propagation direction. The phase difference between the stations in the profile under analysis can then be converted into phase velocities by using the corrected propagation distance. The array analysis requires a minimum of three stations installed in a 'triangle-like' geometry, and in our case it was typically carried out using data from three or four stations. The obtained dispersion curves were subsequently transformed into a 2-D phase velocity map at nine selected periods T, equally sampled between T = 20 s and T = 100 s, using a simple tomographic approach which is described by a limited number of parameters (see next section).

The combination of array analysis and two-station measurements does not dispense with the need for very careful data selection and analysis. We selected 65 teleseismic events distributed with backazimuths mostly between -30° to 80° and 120° to 200°, with magnitudes larger than 5.5, in a distance range from 17° to 159° and a good signal-to-noise ratio. All records were corrected for the instrument response. Multiple filter analysis (Dziewonski et al. 1969) and subsequent frequency-time filtering (Levshin 1989) were applied to Rayleigh (vertical component) and Love (transverse component) waves to minimize contamination by noise, spectral holes, higher modes and multipathing. The frequency-time filtering was carried out with manual control at all steps. Records with a low signal-to-noise ratio and poorly constrained group velocities were eliminated. Approximately 40 per cent of the initial data were rejected. The careful data selection and processing make it possible in some cases to extend the dispersion measurements beyond the seismometer's cut-off frequency. The time differences used in the phase velocity analysis were calculated using Wiener filtering (e.g. Taylor & Toksoz 1982). Only phase differences estimated with a coherence higher than 0.95 were accepted.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The profiles for which good-quality phase velocity measurements were obtained for (a) Rayleigh waves (25 profiles) and (b) Love waves (9 profiles in total). Phase velocity dispersion curves are determined by a combination of two-station and array measurements.</figcaption></figure>

For each interstation profile we calculated the average phase velocity and the standard deviation using all the available phase velocity measurements for the profile. The profiles for which we obtained phase velocities for Rayleigh and Love waves are shown in Fig. 2. For Rayleigh waves, the profiles of phase velocity measurements cover the region rather uniformly and can be used for tomography, in spite of a rather poor east-west resolution. For Love waves, phase velocities are measured along some profiles only and cannot be used for separate Love wave tomography. Therefore, we combined Love and Rayleigh wave phase velocity measurements to detect areas of Love-Rayleigh discrepancy in the Aegean Sea (see Section 4.3).
It was not possible to obtain phase velocity measurements along the three profiles between the three GEOFON stations (GVD, KRIS, SANT), for three reasons: (1) the stations are located very close to each other; (2) the signal-to-noise ratio at long periods was generally low; and (3) of the events with a good signal-to-noise ratio, the majority suffered from strong multipathing effects in this highly heterogeneous area. Similar effects were observed by Shapiro *et al.* (1998) along the Mexican coast. We therefore calculated the average structure beneath these stations using the array technique of Pedersen *et al.* (2003). However, with only one high-quality event, the measurement may be biased. We therefore attributed high error bars to this particular dispersion curve, to decrease its influence on the subsequent tomography. We verified that the inclusion of these stations did not alter the tomography results in other areas which are covered by better measurements, including the area immediately north of Crete.

Fig. 2(a) shows that the number of Rayleigh wave phase velocity measurements used for the tomography is rather limited. The profiles of phase velocity measurements are mostly oriented east-west so we can hope to obtain good constraints on the shear wave velocity structure mainly in the north-south direction.

The phase velocity measurements for Rayleigh and Love waves in the period range 20–100 s are presented in Figs 3 and 4. Some phase velocity dispersion curves, especially in the southern part of the region, have large error bars. This is mostly due to the location of the seismological stations on islands and poor weather conditions during the experiment that result in a strong microseismic noise.

# 3 INVERSION METHODS

We applied a two-stage inversion to obtain constraints on the 3-D *S*-velocity structure from the dispersion curves of Fig. 3. In the first stage we used a 2-D tomographic method to obtain maps of phase velocities in the area for a series of periods. In the second stage, the locally averaged dispersion curves at different points were inverted to obtain *S*-wave velocity profiles with depth.

## 3.1 Inverting for phase velocity maps

To solve the 2-D tomography problem we used the method developed by Ditmar & Yanovskaya (1987). It is a generalization of the Backus–Gilbert method for the 1-D case (Backus & Gilbert 1968), when the data kernels are not square integrable. The resulting 2-D phase velocity maps are made to be smooth by including the norm of the phase velocity gradient in the cost function (Yanovskaya & Ditmar 1990). The simple tomographic method that we use is particularly well adapted to our limited data set because the number of parameters to invert for is not very high: this method estimates only the phase velocity distribution from traveltimes without taking into account other factors which can affect these traveltimes such as wave front distortion and diffraction. It can provide a stable result for the large-scale lateral velocity variations in the region, but due to the limited number of data and dominance of east-west-oriented profiles, small-scale lateral heterogeneities cannot be resolved.

The sphericity is taken into account by an Earth-flattening transformation (Yanovskaya 1982; Jobert & Jobert 1983). The region of the sphere covered by the rays is transformed into a plane with the saved invariance of traveltimes. The solution is constructed on the plane, and then the plane is transformed back into the sphere. This technique is acceptable for our region of study as it can be

used in areas up to approximately 2000 × 2000 km² (e.g. Barmin *et al.* 2001).

Measurement errors, multipathing and interference of Rayleigh and Love waves can all lead to incompatible traveltime observations in the region. To reduce the effect of individual phase velocities which for some reason (e.g. multipathing, noise) are incompatible with those from neighbouring profiles, we ran the tomographic inversion several times, at each run eliminating profiles where the traveltime residual between the predicted and observed phase velocity was more than 3 rms (Yanovskaya *et al.* 1998; Karagianni *et al.* 2002). In our case, the number of eliminated paths varies from 0 to 2 depending on period.

The output velocity maps correspond to velocities that are averaged over some area which depends on the profile density and orientation. Since in our case the paths have a predominantly east-west direction, the averaging areas should also be stretched in this direction. Following the procedure proposed by Yanovskaya (1997), we characterize the averaging areas by ellipses which are characterized by their size and the direction of the major axis (e.g. Yanovskaya *et al.* 1998; Karagianni *et al.* 2002).

## 3.2 Inversion for *S* velocity with depth

For the subsequent 1-D inversion of local dispersion curves we developed a variant of the two-step inversion technique of Shapiro *et al.* (1997). In the original method, a linearized inversion of the dispersion curve is carried out (Herrmann 1987) to obtain a reasonable layered starting model, after which random changes are applied to obtain a family of models which can explain the observed dispersion curve. The main drawback of this method is that mantle layers are uncorrelated so that highly oscillating and unrealistic models may be part of the final family of mathematically accepted models. The uncertainty of the earth model, which is obtained as the average of all the acceptable models, is consequently artificially high. The unrealistic models must then either be manually or automatically rejected. For larger areas, it is possible to reduce the area of search to models that are constrained by other data such as heat flow (Shapiro & Ritzwoller 2004).

To take into account the continuity of velocities in the mantle we modified the method of Shapiro *et al.* (1997). Practically, the model from which the Monte Carlo inversion starts is constructed as follows: the velocity profile which results from the linearized inversion is approximated in the mantle by a small number of layers. The velocity gradient is constant within each layer and no velocity jumps are allowed across layer interfaces. As the mantle velocity model is continuous, only two parameters (velocity at the top of the layer, depth to the bottom of the layer) are necessary to describe the velocity function in each layer. In the crust, the model is parametrized by the thickness and constant velocity of each layer. With this parametrization, the number of parameters is small, making the subsequent Monte Carlo inversion fast in terms of CPU time. The difference from the method of Shapiro *et al.* (1997) resides in this different description of the mantle part of the model.

At each iteration of the Monte Carlo inversion, random changes within given limits are added to all model parameters. If the new predicted dispersion curve falls within the error bars of the ‘observed’ dispersion curve, i.e. that obtained after the 2-D tomographic inversion curve, the model is accepted and used as the initial model for next iteration. In the opposite case, the new gradient model is rejected and the next iteration is operated from the last accepted model. To be certain of exploring different parts of the parameter
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Rayleigh wave mean phase velocity dispersion curves and their standard deviation for each individual profile (vertical axis: phase velocity in km s⁻¹, horizontal axis: period in s). The number in the lower right-hand corner of each figure indicates the number of phase velocity measurements for the profile.</figcaption></figure>

space, we run the Monte Carlo inversion several times, with different random series.

The main advantage of the alternative parametrization is that only a few unrealistic models are part of the final set of accepted models. Consequently, the error bars on the model, calculated as the standard deviation over all accepted models, become realistic without requiring a *posteriori* rejection of unrealistic solutions.

# 4 RESULTS

## 4.1 Tomographic images

We calculated the lateral variations of phase velocities beneath the Aegean by applying the 2-D tomographic method described in the previous section for the period range of 20–100 s.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Fig. 3 for Love wave dispersion curves.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of mean phase velocity dispersion curves of Rayleigh waves in the Aegean region. Phase velocity dispersion curves for Papazachos (1969), Papazachos *et al.* (1995) and IASP91 are calculated from corresponding S-wave velocity models. Our mean dispersion curve differs significantly from the curve of Papazachos *et al.* (1995) between 20 and 70 s with a maximum difference at 40 s.</figcaption></figure>

As a starting model for the 2-D tomography, we used the average of the observed dispersion curves of Fig. 3. Our average dispersion curve is shown in Fig. 5. For comparison, we added the dispersion curves that correspond to the S-wave models suggested by Papazachos (1969), Papazachos *et al.* (1995) and IASP91 (Kennett & Engdahl 1991). Compared with these regional and global reference dispersion curves, our dispersion curve is characterized by low velocities at intermediate periods. The shift is too big to be

explained by realistic differences in the average Moho depth in the area and must therefore reflect average differences in uppermost mantle structure between the investigated areas.

Due to the limited number of phase velocity measurements the tomographic inversion was damped towards a very smooth model. Rayleigh wave phase velocities obtained from 2-D tomography are shown in Fig. 6 for selected periods. The results for 30, 40, 60, 80 and 90 s are omitted, as they are similar to the presented maps of 20, 50, 70 and 100 s. We present here the results only where the mean size of the averaging area is less than 400 km. Note again that this value of the averaging area in most places covers big differences in north-south and east-west resolution.

The quality of the solution can be estimated by comparing initial and remaining rms traveltime residuals. The initial mean square traveltime residual for our data set is 3.1 s and the mean remaining residual is 1.5 s, i. e. approximately half of the observed residuals can be explained by the lateral heterogeneities of Fig. 6. This relatively modest residual reduction is mainly due to the strong smoothing.

The resolution of the inversion was estimated by the averaging ellipses which are displayed in Fig. 7 in some selected points for Rayleigh waves for the 70 s period. As expected, their long axis is parallel to the predominant east-west profile orientation. Average resolution varies from approximately 200 km, typically in the north-south direction to up to 800 km in the east-west direction.

Phase velocities at periods of 20–30 s correspond to a layer of approximately 40 km thickness and are strongly influenced by the crustal velocity structure and in particular by Moho depth. The crustal thickness varies considerably across the Aegean. In continental Greece, the Moho reaches a depth of 40 km (e.g. Kalogeras & Burton 1996), whereas in the Sea of Crete it is less than 25 km deep (e.g. Karagianni *et al.* 2002). In the subsequent inversions, we use the preliminary Moho depth map of D. Hatzfeld (personal communication) which is based on receiver function analysis and
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Phase velocities of Rayleigh wave for periods of 20, 50, 70 and 100 s. As a rough approximation, the tomographic maps for 20, 50, 70 and 100 s sample the Earth from the surface to depths of 30, 100, 140 and 200 km respectively. Phase velocities of a period of 20 s are strongly influenced by the variations in crustal thickness. The Moho depth variations (D. Hatzfeld, personal communication) are superimposed on the phase velocity map at *T* = 20 s. For periods of 40 s and longer a high-velocity anomaly is observed in the southern Aegean which we attribute to the Hellenic subduction zone.</figcaption></figure>

which has the advantage of covering the whole study region. These Moho depths are superimposed onto the phase velocity map of Fig. 6 (*T* = 20 s). Ongoing work on receiver functions (Sodoudi *et al.* 2004) should soon provide even better Moho maps of the area.

There is a good correlation between the phase velocity anomalies in Fig. 6 (*T* = 20 s) and the Moho depths. Low phase velocities are generally associated with crustal thickening. These anomalies (that reach a value of 1.5–2 per cent) are consistent with the results of the tomography by Papazachos & Nolet (1997) even though their depth slices show more details than our velocity maps. These differences arise from the fact that only large-scale crustal thickness variations are identified in our tomography and the phase velocity map for the 20 s period corresponds to an integration over the velocity structure in the crust and, to a lesser extent, the uppermost mantle.

In the southern Aegean we observe a high-velocity anomaly of approximately 3 per cent. This anomaly, present for all periods

from 40 to 100 s and trending east–west, can be attributed to the Hellenic subduction zone. The lowest velocity zone on 40–100 s tomographic maps is found beneath the northern Aegean Sea and shows a trend parallel to the North Aegean Trough. This large tectonic structure is the continuation of the North Anatolian Fault within the northern Aegean Sea. It appears on the maps (Fig. 6, *T* = 50, 70 and 100 s) that the Aegean Trough is associated with southwest–northeast-trending low-velocity anomaly. However, the shape of the anomaly may be influenced by the orientation of the profiles.

## 4.2 S-wave velocity structure

To ensure that we used the mean 1-D S-wave velocity profile beneath the study area, we recalculated the average dispersion curve from
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Averaging areas of elliptic shape at selected points (+) for Rayleigh waves for a period of 70 s. The azimuth-dependent resolution at any point (x₀, y₀) is approximated by an ellipse centred at (x₀, y₀) with axes oriented along the directions of the best and the worst resolution. The long axes of the ellipses are parallel to the predominant profile orientation. The average resolution varies from 200 km mostly in a north-south direction to 800 km in an east-west direction.</figcaption></figure>

the phase velocity maps of Fig. 6. It only differs very slightly (and within the error bars) from the dispersion curve in Fig. 5. We then inverted the updated average dispersion curve to obtain the mean S-wave velocities as a function of depth. This mean model is calculated by the two-step inversion described in Section 3. IASP91 (Kennett & Engdahl 1991) was used as starting model in the linearized inversion. In the subsequent Monte Carlo inversion, the random perturbations that were applied to the model parameters at each iteration were up to 1 per cent for S-wave velocities and up to 5 per cent for interface depths. By interface depths we mean here the depths where the velocity gradient changes. In this and all other inversions we repeated the Monte Carlo inversion with 50 series of random numbers to be sure that we explored different parts of the parameter space. Each Monte Carlo inversion had 50 iterations, so a total of 2500 models were tested in the inversion.

At each depth we calculated the average S-wave velocity of all accepted models shown in Fig. 8. In Fig. 8 we also plot in grey all acceptable models obtained by the stochastic inversion. The error bars are approximately 2 per cent in the mantle. Our mean S-wave model is compared with the one suggested by Papazachos *et al.* (1995) for the Aegean region and IASP91 (Kennett & Engdahl 1991). As compared with Papazachos *et al.* (1995), our mean S-wave velocities are lower in the uppermost mantle between 20 and 80 km. As the S-wave model of Papazachos *et al.* (1995) includes information from continental Greece which is not well sampled by

our data, we conclude that the S-wave model of Papazachos *et al.* (1995) is not well adapted for the central Aegean region between depths of 20 and 80 km. The low uppermost mantle S velocities are rather surprising, as one would expect the influence of the slab to increase the average velocities for the area. This is a first indication of the presence of very low velocities in other parts of the area.

To visualize the 3-D S-wave velocity structure beneath the Aegean we applied the same inversion procedure to the local dispersion curves obtained from the 2-D tomography in a 0.6° × 0.6° grid. These models were smoothed horizontally and assembled into S-wave velocity maps for selected depths.

Fig. 9 shows the absolute S-wave velocity distribution in the Aegean for depths of 50, 70, 100 and 200 km. The maximum wavelength of our data set is approximately 450 km so the S-velocity distribution is resolved down to approximately 200 km. As we used dispersion curves in the range of periods between 20 and 100 s, details of crustal structure are not resolved. Instead, we focus on lateral velocity variations that are associated with differences in the mantle structure. The resolution of lateral velocity variations varies from approximately 200 km to 800 km depending on azimuth, with the best resolution in the north-south direction.

As the east-west resolution is poor, a north-south vertical cross-section is a better representation of our results. Fig. 10 shows a cross-section along the profile AA'.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Mean S-wave velocity model for the Aegean region. The mean S-wave velocity model of this study corresponds to the mean dispersion curve in the Aegean, inverted for S-wave velocity by the two-step Monte Carlo inversion technique (Section 3.2). All acceptable models obtained using the stochastic Monte Carlo inversion are presented in grey (in total 274 models). The error bars of our mean S-wave velocity model are calculated as the standard deviation for all accepted models, and are up to 2 per cent. Compared with Papazachos *et al.* (1995), our mean S-wave velocities are lower in the uppermost mantle, between 20 and 80 km.</figcaption></figure>

The main feature on Figs 9 and 10 is the high-velocity anomaly associated with the slab. We have, however, a limited resolution of the slab geometry, which is better shown by deep body wave tomography (Spakman *et al.* 1988). In Fig. 10 the slab appears subhorizontal with an apparent ‘break-off’ at a depth of 100–150 km. These features are artefacts created by the strong smoothing imposed in the tomographic inversion. A significantly better path coverage is required to obtain constraints on the slab geometry. The S-wave velocity contrast observed in the slab is up to approximately +4 per cent compared with IASP91 (Kennett & Engdahl 1991) and +3 per cent compared with the area average. These values may be somewhat underestimated because of the lateral smoothing in the 2-D tomographic inversion.

The amplitude of the velocity variations associated with the slab and the absolute velocities are in agreement with previous studies based on local and teleseismic tomography (Spakman *et al.* 1988; Papazachos & Nolet 1997; Tiberi *et al.* 2000). The absolute shear velocities in the slab at a depth of 100 km are up to 4.68 km s⁻¹ (Fig. 9) corresponding to a *P* velocity of 8.09 km s⁻¹ (assuming a Poisson ratio of 0.25) which is in good agreement with the *P* velocity of 8.1 km s⁻¹ obtained in the same location by Papazachos & Nolet (1997). S-wave velocities of 4.7–4.8 km s⁻¹ in the slab (Meier *et al.* 2004) in Crete and the Sea of Crete are also very close to our results.

In the northern part of the Aegean region, at a distance of approximately 500 km from the Hellenic Trench, we observe a significant low-velocity anomaly situated between 50 and 130 km

depth (Fig. 10). The anomaly resolved by the data is 250 km wide; however, the resolution is insufficient to determine its shape. The S-wave velocity contrast compared with the average of the whole area is approximately –3 per cent at 100 km depth. However, absolute velocities are only low down to approximately 100 km, where they meet IASP91. The apparent low-velocity zone visible at 200 km depth in Fig. 9 in the northern Aegean actually corresponds to a +2 per cent contrast as compared with IASP91. Fig. 9 shows ‘low’ velocities in this zone only by contrast with a strong high-velocity anomaly in the southern part associated with the slab.

Within the low-velocity layer, the maximum velocity difference from IASP91 is up to –4 per cent, at 70 km depth. Even though any error in Moho depth would translate into erroneous mantle velocities, such problems cannot explain our observed low velocities: the most recent results from that particular area (F. Sodoudi, personal communication) show that the Moho may even be a little bit shallower, so the low-velocity layer in the upper mantle would become even stronger. This anomaly is also detected in a recent tomography of the whole of the Mediterranean by Marone *et al.* (2004). The location of the low-velocity anomaly coincides with the northern part of the Aegean, which is an area of rapid deformation and intense extension.

## 4.3 Love–Rayleigh discrepancy

A detailed analysis of anisotropy in the Aegean region is not possible as there are too few Love wave phase velocity measurements to perform a tomographic analysis (see Fig. 2b). However, it is still possible to compare Rayleigh and Love wave dispersion curves to identify zones of Love–Rayleigh discrepancy, which would be an indication of anisotropy. We therefore independently inverted Rayleigh and Love dispersion curves measured along the same profiles.

In a tomographic analysis it is possible to include profiles where only few events are available for the two-station measurements, as incompatible data are rejected during tomographic inversion (see Section 4.1). As we do not have such a *posteriori* control when we use individual profiles, we can use only the highest-quality phase velocities. The most objective approach corresponds to eliminating profiles where only a few events are available. We therefore use only profiles with more than five events for each of the Rayleigh and Love wave phase velocity measurements. The limit of five events is somewhat arbitrary, but it ensured, for example, that the phase velocities are not dominated by one event with several aftershocks. Only nine paths with both Rayleigh and Love dispersion curves were available (Fig. 11).

The results of the inversion for these profiles are presented in Fig. 12. We observe significant vertical anisotropy in the central Aegean Sea (profiles AGGI-LESB, AGGI-SKIR, LESB-SKIR, HIOS-PENT and PENT-SAMO shown as bold lines in Fig. 11). Along these profiles the Love wave-derived velocities are of the order of up to 10 per cent higher than those derived from Rayleigh waves in the mantle, down to our limit of resolution. A small or non-existent discrepancy is found in the northern (profile AGGI-LIMN) and southern Aegean (profiles SAMO-VELI and VELI-KOS1). Profile HIOS-VELI shows a high discrepancy, but the error bars are too high for this discrepancy to be reliably resolved.

Love–Rayleigh wave discrepancy is observed in many regions of the Earth. Globally, the lithosphere appears faster to Love waves than to Rayleigh waves. In global models, the velocity difference is about 4 per cent in PREM (Dziewonski & Anderson 1981) and 3 per cent in AK135 (Kennett *et al.* 1995). On a regional scale the Love–Rayleigh discrepancy may reach much larger values. For example,
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Absolute *S*-wave velocity distribution at depths of 50, 70, 100 and 200 km in the Aegean. The *S*-wave velocity maps are assembled from average models obtained by Monte Carlo stochastic inversion at each point in the region. The Hellenic Trench and the North Anatolian Trough are superimposed on the maps (see Fig. 1 for details). The resolution of lateral velocity variations varies from approximately 200 km to approximately 800 km depending on azimuth, with the best resolution generally being in the north–south direction.</figcaption></figure>

in Australia the Love–Rayleigh velocity discrepancy is up to 9 per cent (Debayle & Kennet 2000), close to our observations in Aegean Sea. This velocity difference cannot be explained by the orientation of olivine crystals in the mantle (Maupin 2003) or by the effect of small-scale lateral heterogeneities (Maupin 2002).

Even though an explanation of very high Love–Rayleigh discrepancies remains to be found, it is intriguing that the area of high discrepancy in the Aegean corresponds to an area of strong SKS splitting (Hatzfeld et al. 2001), also shown in Fig. 11, which is generally consistent with the present day deformation as measured by GPS (McClusky et al. 2000). The SKS delays are in places up to 2 s and can be interpreted as a 7 per cent anisotropy in the lithosphere, assuming the anisotropic layer is 100 km thick. However, this anisotropy would yield a smaller Love–Rayleigh discrepancy than we observe. Furthermore, our results show Love–Rayleigh discrepancy to a depth of at least 150 km, i.e. possibly into the

asthenosphere. One can speculate whether the *SKS* splitting and the Love–Rayleigh discrepancies are created by the same dynamic causes in the Aegean Sea, even though the physical mechanism of the anisotropy must be at least partly different.

## 5 DISCUSSION AND CONCLUSION

Teleseismic surface wave tomography provided a 3-D image of the lithosphere beneath the Aegean Sea down to 200 km. The resolution analysis shows that only large-scale lateral anomalies greater than 200 km can be revealed from this data set. The resolution is in most places better in the north–south direction than in the east–west direction.

In this study we obtained a well-constrained mean *S*-wave velocity model for the Aegean Sea. Our mean model shows lower
<figure><img src="image_11.png" /><figcaption>**Figure 10.** S-wave velocity perturbation cross-section (without any vertical exaggeration) oriented north–south, in the direction of the best resolution. A high-velocity zone (up to +3 per cent) is associated with the subducting slab. The low-velocity zone discussed in the text is indicated on the figure by LV. The position of the North Anatolian Trough is indicated as NAT on the profile. The geometry and the apparent discontinuity of the slab cannot be resolved by our data set. The velocity variations are probably underestimated because the tomographic inversion was damped toward a very smooth model.</figcaption></figure>

velocities between the depths of 20 and 80 km depth than in the model suggested by Papazachos *et al.* (1995).

Depth slices and cross-sections of shear wave velocities of the region under study show a high-velocity anomaly of 3 per cent associated with the subducting slab. The amplitude of this anomaly is generally in agreement with previous P-wave tomographic studies in the Aegean (Spakman *et al.* 1988; Papazachos & Nolet 1997; Tiberi *et al.* 2000), even though the velocity variations could be underestimated in our study because the solution of the 2-D tomographic inversion was damped towards a smooth model. The precise shape and position of the slab cannot be resolved because of an insufficient lateral resolution.

In addition, a large low-velocity anomaly located 500 km from the Hellenic Trench is observed at a depth of 50–100 km. This anomaly is 250 km wide with a velocity contrast of up to –3 per cent compared with the mean of the area and up to –4 per cent compared with IASP91. Its shape and size cannot be determined, as we applied strong smoothing and as the majority of profiles are east–west or northeast–southwest trending. The velocity decrease is resolved by the data, but the exact depth and amplitude of it are somewhat uncertain due to the non-uniqueness of the solution.

This low-velocity anomaly cannot be explained by the mantle wedge observed in the P-wave tomography study of Papazachos & Nolet (1997) in the Aegean subduction zone. The mantle wedge is situated beneath the volcanic arc and the Gulf of Corinth (200–300 km from the Hellenic Trench), at a depth of approximately 60 km. The lateral extent of this area, associated with the dehydration process of the slab, is approximately 70–100 km. An anomaly of such a small size cannot be resolved by our surface wave analysis, or would show up as a broader, but much smaller-amplitude, anomaly.

According to Goes *et al.* (2000), a –3 per cent shear wave velocity at 100 km depth can be explained by a temperature variation of approximately +200 °C, which can be considered as a good first-order estimation of a pure temperature anomaly beneath this area. However, small amounts of fluids would be able to produce a similar anomaly (e.g. Nolet & Zielhuis 1994).

The low-velocity zone observed in our study is located immediately below the zone of intense extensional deformation at the surface (McClusky *et al.* 2000). Moreover, this anomaly must be large, as it would otherwise not be detected by our analysis. Low-velocity anomalies in connection with surface extension were observed by Marone *et al.* (2004) in the Algero-Provencal and Tyrrhenian basins of the Mediterranean region and were interpreted as a shallow asthenospheric layer. In our study there is no evidence that the low-velocity zone is created by an asthenospheric layer in the northern Aegean. The size of this anomaly combined with its location immediately below the area of intense surface deformation support a model of a distributed deformation of the upper mantle in this particular area.

South of this area, there is evidence of very strong Love–Rayleigh discrepancy, and the resulting isotropic models are different down to a depth of at least 150 km. The difference between Love- and Rayleigh-derived models is up to 10 per cent, providing further evidence that an aligned olivine may not be the only cause of Love–Rayleigh discrepancy. It is, however, intriguing that the Love–Rayleigh discrepancy is strongest in the area of strong SKS wave splitting as measured by Hatzfeld *et al.* (2001), even though the values they obtain (azimuthal anisotropy of 7 per cent if 100 km thickness is used) cannot explain the anisotropy that we observe.
<figure><img src="image_12.png" /></figure>

**Figure 11.** Distribution of azimuthal anisotropy deduced from SKS splitting (Hatzfeld *et al.* 2001) and selected profiles for S-wave velocity inversion of both Rayleigh and Love dispersion curves corresponding to the same propagation path used in anisotropy analysis. Profiles with strong Rayleigh–Love discrepancy (AGGI–LESB, AGGI–SKIR, LESB–SKIR, HIOS–PENT and PENT–SAMO) are plotted as bold lines. Profiles with small or non-existent discrepancy (AGGI–LIMN, SAMO–VELI and VELI–KOS1) are indicated as dashed lines. The profile HIOS–VELI, here plotted as a dashed line, may have strong Love–Rayleigh discrepancy, but this is not fully resolved due to large uncertainties in both the Rayleigh and Love depth inversions.

## ACKNOWLEDGMENTS

This work was done while the first author was a PhD fellow-
ship student as a part of French Government collaboration pro-
gramme between the Joseph Fourier University of Grenoble, France
and the University of St Petersburg, Russia. This study was sup-
ported by contract ENV4-CT96-0277 in the program ‘Environ-
ment and Climate’ of the European Commission, and the pro-
gram IDYL of INSU-CNRS. We thank K. Makropoulos, B. Papaza-
chos, H. Lyon-Caen and K. Priesley for their efforts to make the
field experiment successful. We also thank the GeoForschungsZen-
trum Potsdam, Germany for the GEOFON data from GVD, KRIS
and SANT stations. We are indebted to D. Baumont and N.
Shapiro for providing the code of array phase velocity measure-
ments and Monte Carlo stochastic inversion respectively. We are
grateful to Eric Debayle and an anonymous reviewer for com-
ments that helped to significantly improve the manuscript. E.
Bourova thanks Professor Tetsuo Takanami, Dr Kazunori Yoshizawa
and Professor Yoshinori Furukawa (Hokkaido University, Japan)
for their support during the preparation of the revised version
of this manuscript. She also wishes to thank the Japan Soci-
ety for the Promotion of Science (JSPS) for the financial sup-
port. Most of the figures were plotted with the Generic Map-
ping Tools (GMT) data processing and display package (Wessel &
Smith 1991, 1995).

## REFERENCES

Armijo, R., Meyer, B., King, G.C.P., Rigo A. & Papanastassiou, D., 1996. Quaternary evolution of the Corinth Rift and its implications for the late Cenozoic evolution of the Aegean, *Geophys. J. Int.*, **126**, 11–53.

Backus G. & Gilbert F., 1968. The resolving power of gross Earth data, *Geophys. J. R. astr. Soc.*, **16**, 169–205.

Barmin, M.P., Ritzwoller, M.H. & Levshin, A.L., 2001. A fast and reliable method for surface wave tomography, *Pure. appl. Geophys.*, **158**, 1351–1375.

Baumont, D., Paul, A., Beck, S., Zandt, G. & Pedersen, H., 2002. Litho-spheric structure of the central Andes based on surface wave dispersion, *Geophys. J. Int.*, **107**(B12), 103–115.

Bokelmann, G.H.R., 2002. Which forces drive North America?, *Geology*, **30**(11), 1027–1030.

Bourne, S.J., Arnadottir, T., Beavan, J., Darby, D.J., England, P.C., Parson, B., Walcott, R.I. & Wood, P.R., 1998. Crustal deformation of the Marlborough fault zone in the South Island of New Zealand: geodetic constraints over the interval 1982–1994, *J. geophys. Res.*, **103**(B12), 30147–30165.

Bruneton, M., Farra, V., Pedersen, H.A. & the SVEKALAPKO Seismic Tomography Working Group, 2002. Non-linear surface wave phase velocity inversion based on ray theory, *Geophys. J. Int.*, **151**, 583–596.

Cotte, N., Pedersen, H.A., Campillo, M., Farra, V. & Cansi, Y., 2000. Off-great circle propagation of intermediate period surface waves as observed on a dense array in the French Alps, *Geophys. J. Int.*, **142**, 825–840.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** S-wave velocity models obtained by independent inversion of Rayleigh and Love dispersion curves corresponding to the same propagation path. A significant vertical anisotropy is observed in the central Aegean Sea (profiles AGGI-LESB, AGGI-SKIR, LESB-SKIR, HIOS-PENT and PENT-SAMO). The Love wave-derived velocities are 10 per cent higher than those deduced from Rayleigh waves in the mantle. The difference in velocity is approximately 0.5 km s⁻¹ at 100 km depth. Modest anisotropy is found in the northern (profile AGGI-LIMN) and southern Aegean (profiles SAMO-VELI and VELI-KOS1). The discrepancy corresponding to the profile HIOS-VELI is not well resolved.</figcaption></figure>

Debayle, E. & Kennet, B.L.N., 2000. Anisotropy in the Australian upper mantle from Love and Rayleigh waveform inversion, *Earth planet. Sci. Lett.*, **184**, 339–351.

Ditmar, P.G. & Yanovskaya, T.B., 1987. Generalization of Backus-Gilbert method for estimation of lateral variations of surface wave velocities, *Izv. Phys. Solid Earth*, **23**(6), 470–477.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **184**, 297–356.

Dziewonski, A., Bloch, S. & Landisman, N., 1969. A technique for the analysis of transient seismic signals, *Bull. seism. Soc. Am.*, **59**, 427–444.

England, P. & Jackson, J., 1989. Active deformation of the continents, *Annu. Rev. Earth Planet. Sci.*, **17**, 197–226.

Friederich, W., 1998. Wave-theoretical inversion of teleseismic surface waves in a regional network: phase velocity maps and a 3-D upper mantle shear wave velocity model for southern Germany, *Geophys. J. Int.*, **132**, 203–225.

Friederich, W., Wieland, E. & Stange, S., 1994. Non-plane geometries of seismic surface wavefields and their implications for regional surface-wave tomography, *Geophys. J. Int.*, **119**, 931–948.

Gautier, P. & Brun J.P., 1994. Crustal-scale geometry and kinematics of late-orogenic extension in the central Aegean (Cyclades and Evvia Island), *Tectonophysics*, **238**, 399–424.
Goes, S., Govers, R. & Vacher, P., 2000. Shallow mantle temperature under Europe from P and S wave tomography, *J. geophys. Res.*, **105**(B5), 11153-11169.

Hatzfeld, D. et al., 2001. Shear-wave anisotropy in the upper mantle beneath the Aegean related to internal deformation, *J. geophys. Res.*, **106**(12), 30737-30753.

Herrmann, R.B., 1987. *Computer Programs in Seismology, Volume IV: Surface Wave Inversion*, St Louis University, St Louis, MO.

Jackson, J.A. & McKenzie, D.P., 1988. The relationship between plate motions and seismic moment tensors, and the rates of active deformation in the Mediterranean and Middle East, *Geophys. J.*, **93**, 45-73.

Jobert, N. & Jobert, G., 1983. An application of ray theory to the propagation of waves along a laterally heterogeneous spherical surface, *Geophys. Res. Lett.*, **10**, 1148-1151.

Kalogeras, I.S. & Burton, P.W., 1996. Shear-wave velocity model from Rayleigh wave dispersion in the broader Aegean area, *Geophys. J. Int.*, **125**, 679-695.

Karagianni, E.E. et al., 2002. Rayleigh wave group velocity tomography in the Aegean area, *Tectonophysics*, **358**, 187-209.

Kennett, B.L.N. & Engdahl, E.R., 1991. Traveltimes for global earthquake location and phase identification, *Geophys. J. Int.*, **105**, 429-465.

Kennett, B.L.N., Engdahl, E.R. & Buland, R., 1995. Constraints on seismic velocities in the Earth from travel-times, *Geophys. J. Int.*, **122**, 108-124.

Levshin, A.L., 1989. Surface waves in media with weak lateral inhomogeneity, in, *Seismic Surface Waves in a Laterally Inhomogeneous Earth*, pp. 35-69, ed. Keilis-Borok, V.I., Kluwer Academic Publishers, Dordrecht.

Makris, J., 1978. The crust and upper mantle of the Aegean region from deep seismic soundings, *Tectonophysics*, **46**, 269-284.

Makris, J. & Stobbe, C., 1984. Physical properties and state of the crust and upper mantle of the eastern Mediterranean Sea deduced from geophysical data, *Mar. Geol.*, **55**, 347-363.

Marone, F., van der Lee, S. & Giardini, D., 2004. Three-dimensional upper-mantle S-velocity model for the Eurasia-Africa plate boundary region, *Geophys. J. Int.*, **158**, 109-130.

Maupin, V., 2002. The amplitude of the Love-Rayleigh discrepancy created by small-scale heterogeneities, *Geophys. J. Int.*, **150**, 58-64.

Maupin, V., 2003. Surface-waves in an heterogeneous anisotropic continental lithosphere, *Geophys. Research Abstracts*, **5**, 03645.

McClusky, S. et al., 2000. Global Positioning System constraints on plate kinematics and dynamics in the Eastern Mediterranean and Caucasus, *J. geophys. Res.*, **105**(B3), 5695-5719.

Meier, T., Dietrich, K., Stockert, B. & Harjes, H.-P., 2004. One-dimensional models of shear wave velocity for the eastern Mediterranean obtained from the inversion of Rayleigh wave phase velocities and tectonic implications, *Geophys. J. Int.*, **156**, 45-58.

Mercier, J.L., Sorel, D. & Simeakis, K., 1987. Changes in the state of stress in the overriding plate of a subduction zone: the Aegean arc from the Pliocene to the present, *Ann. Tecton.*, **1**, 20-39.

Nolet, G. & Zielhuis, A., 1994. Low S velocities under the Torniquist-Teisseyre zone: evidence for water injection into the transition zone by subduction, *J. geophys. Res.*, **99**(B8), 15813-15820.

Papazachos, B.C., 1969. Phase velocities of Rayleigh waves in southeastern Europe and eastern Mediterranean Sea, *Pure. appl. Geophys.*, **75**, 47-55.

Papazachos, C.B. & Kiratzi, A.A., 1996. A detailed study of the active crustal deformation in the Aegean and surrounding area, *Tectonophysics*, **253**, 129-153.

Papazachos, C. & Nolet, G., 1997. P ans S deep velocity structure of the Hellenic area obtained by robust nonlinear inversion of travel times, *J. geophys. Res.*, **102**(B4), 8349-8367.

Papazachos, C.B., Hatzidimitriou, P.M., Panagiotopoulos, D.G. & Tsokas, G.N., 1995. Tomography of the crust and upper mantle in southeast Europe, *J. geophys. Res.*, **100**(B7), 12405-12422.

Pedersen, H.A., Coutant, O., Deschamps, A., Soulage, M. & Cotte, N., 2003. Measuring surface wave phase velocities beneath small broad-band arrays: tests of an improved algorithm and application to the French Alps, *Geophys. J. Int.*, **154**, 903-912.

Shapiro, N. & Ritzwoller, M., 2004. Thermodynamic constraints on seismic inversions, *Geophys. J. Int.*, **157**, 1175-1188.

Shapiro, N., Campillo, M., Paul, A., Singh, S.K., Jongmans, D. & Sanchez-Sesma, F.J., 1997. On the origin of long period seismic-wave amplification in the Valley of Mexico, *Geophys. J. Int.*, **128**, 151-166.

Shapiro, N.M., Campillo, M., Singh, S.K. & Pacheco, J., 1998. Evidence of trapped seismic surface waves in the Mexican subduction zone, *Geophys. Res. Lett.*, **25**(1), 101-104.

Sodoudi, F., Bock, G., Hatzfeld, D., Kind, R., Hanka, W., Wylegalla, K., Stavrakakis, G. & Priestley, K., 2004. A receiver function study in Greece and the Hellenic subduction zone, *Geophys. Research Abstracts*, **6**, 01330.

Spakman W., Wortel, M.J.R. & Vlaar, N.S., 1988. The Hellenic subduction zone: a tomographic image and its geodynamical implications, *Geophys. Res. Lett.*, **15**(1), 60-63.

Taylor, S.R. & Toksoz, M.N., 1982. Measurement of interstation phase and group velocities and Q using Wiener filtering, *Bull. seism. Soc. Am.*, **72**, 73-91.

Thatcher, W., 1995. Microplate versus continuum descriptions of active tectonic deformation, *J. geophys. Res.*, **100**(B3), 3885-3894.

Tiberi, C. et al., 2000. Crustal and upper mantle structure beneath the Corinth rift (Greece) from a teleseismic tomography study, *J. geophys. Res.*, **105**(B12), 28159-28171.

Tsokaz, G.N. & Hansen R.O., 1997. Study of the crustal thickness and the subducting lithosphere in Greece from gravity data, *J. geophys. Res.*, **102**(B9), 20585-20597.

Wessel, P. & Smith, W.H.F., 1991. Free software helps map and display data, *EOS, Trans. Am. geophys. Un.*, **72**, 441.

Wessel, P. & Smith, W.H.F., 1995. New version of the Generic Mapping Tools released, *EOS, Trans. Am. geophys. Un.*, **76**, 329.

Wielandt, E., 1993. Propagation and structural interpretation of non-plane waves, *Geophys. J. Int.*, **113**, 45-53.

Yanovskaya, T.B., 1982. Distribution of surface group velocities in the North Atlantic, *Izv., Phys. Solid Earth*, **2**, 3-11.

Yanovskaya, T.B., 1997. Resolution estimation in the problems of seismic ray tomography, *Izv., Phys. Solid Earth*, **33**(9), 762-765.

Yanovskaya, T.B. & Ditmar, P.G., 1990. Smoothness criteria in surface wave tomography, *Geophys. J. Int.*, **102**, 63-72.

Yanovskaya, T.B., Kizima E.S. & Antonova, L.M., 1998. Structure of the crust in the Black Sea and adjoining regions from surface wave data, *J. Seism.*, **2**, 303-316.