Source parameters of the M 6 . 5 Skyros Island ( North Aegean Sea ) earthquake of July 26 , 2001

Teleseismic body wave modelling, time domain moment tensor inversion of regional waveforms and spectral analysis of the far-field P-wave pulses are used to derive the source parameters of the July 26, 2001 Skyros earthquake (M 6.5). Its epicentre is located south of the Sporades Islands in the North Aegean Sea (Greece). Previous focal mechanism solutions indicate motion on strike-slip faults. The time domain moment tensor inversion is applied for the first time to the regional waveforms of the recently established broadband network in Greece. Its application gave results which are highly consistent with teleseismic waveform modelling. The results of this study, in combination with the distribution of aftershocks, indicate left-lateral strike slip motion on a NW-SE striking fault with parameters: fault plane (strike = 151°, dip = 83°, rake = 7°) and auxiliary plane (strike = 60°, dip = 84°, rake = 173°), depth 12 km and M0 = 5.98e18 N m. Moreover, the time domain moment tensor inversion technique yielded a pure double couple source with negligible CLVD. The spectral analysis of the far-field P-wave pulses resulted in a fault length L ~ 32 km, stress drop ~ 9 bars and average displacement u ~ 30 cm. These values are in very good agreement with those estimated from empirical scaling relations applicable to the Aegean area.


Introduction
On July 26, 2001 an earthquake of M 6.5 occurred in the Aegean Sea, a few kilometres NW of the island of Skyros (fig.1).Despite its large magnitude and the fact that it was felt in a wide region, the earthquake did not cause extensive damage as its epicentre was located at sea (38.995°N, 24.382°E; Geophysical Laboratory, Aristotle University of Thessaloniki).Limited damage was only observed in the island of Skyros, concerning mainly old residential buildings, a 1000-year-old monastery located inside the island's castle and limited landslides and rock falls.
The aim of this paper is to investigate in detail the source parameters of the Skyros mainshock.We use teleseismic data from the Global Seismograph Network (GSN) as well as regional broadband waveforms recorded by the recently installed broadband network operated by the National Observatory of Athens (NOA).The focal mechanism and related parameters are determined through a comparative application of two different methodologies: teleseismic waveform modelling (Nábelek, 1984)

and time
Mailing address: Prof. Anastasia A. Kiratzi, Department of Geophysics, Aristotle University of Thessaloniki, GR 54124, Thessaloniki, Greece; e-mail: kiratzi@geo.auth.grdomain moment tensor inversion (Dreger, 2000).The latter methodology is applied for the first time in the study of Aegean earthquakes.Furthermore, the dimensions of the causative fault, the stress drop and the average displacement across the ruptured area, are determined from the far-field displacement spectra of the event.

Methodology and application
The technique of body-waveform modelling, described in detail by Nábelek (1984) and McCaffrey et al. (1991), was used to calculate the focal mechanism of the Skyros main-

N o r t h A e g e a n T r o u g h t
Sporades Isl.shock.The employed data consist of P and SH waveforms from stations of the Global Seismograph Network (GSN).All waveforms have a sampling frequency of 1 Hz and were recorded at epicentral distances ranging from 30°t o 90°.The nodal planes of the focal mechanism were constrained by the use of P-wave polarities from stations located at distances less than 30°.The applied methodology has been extensively described in past papers (e.g., Kiratzi and Louvari, 2001;Louvari and Kiratzi, 2001) and therefore we will only present a brief description of it.
A FIR (Finite Impulse Response) band pass filter was first applied to the data, in order to remove the high frequency noise.The corner frequencies of the filter were determined visually, by comparing the FFT amplitude spectrum of the signal with the corresponding spectrum of the noise.In most cases, the frequencies of inter-   est in the Long-Period waveforms covered the range from f 1 = 0.01 Hz to f 2 = 0.1 Hz.The aforementioned frequency limits were used as corner frequencies of the filter.All waveforms were integrated to displacement prior to the inversion.
We used the MT5 version of McCaffrey and Ambers (1988) algorithm to invert P and SH waveform data in order to obtain the strike, dip, rake, centroid depth, seismic moment and source time function of the examined event.The methodology assumes that the source can be represented as a point (the centroid) in space, although not in time.The time history of the displacement on the fault is represented by a source time function made up of a series of overlapping isosceles triangles, whose number and duration are defined by the user.The inversion routine yields amplitudes for each triangular shape.Amplitudes were corrected for the geometrical spreading, the epicentral distance (Langston and Helmberger, 1975), and the attenuation using Futterman's (1962) operator with t* = 1 s for P-and t* = 4 s for SH-waves.The inversion adjusts the relative amplitudes of the source time function elements, the centroid depth, the seismic moment and the source orientation (strike, dip, rake) to minimize the misfit between observed and synthetic seismograms.We refer to this solution as the «minimum misfit solution».The covariance matrix associated with the «minimum misfit solution» usually underestimates the true uncertainties associated with the source parameters.To find more realistic uncertainties, we followed the methodology of McCaffrey and Nábelek (1987) and Molnar and Lyon-Caen (1989) by fixing some of the source parameters at values close to, but different from those of the «minimum misfit solution» and allowing all the other parameters to vary during the inversion.The errors are determined by visual examination when the match of the observed to synthetic seismograms significantly deteriorates.
Synthetics were generated for a point source buried in a half-space.We used a P-wave velocity of 6.5 km/s, S-wave velocity 3.7 km/s and density 2.75 g/cm 3 .We also used a 300 m water layer overlying the half-space.

Inversion results
The focal mechanism of the 26 July 2001 earthquake was computed by inverting 24 Pand 22 SH -Long-Period waves with good azimuthal coverage (fig.2).The results are shown in table I. First-motion polarities of P-waves recorded at 8 stations located at epicentral distances 0°-30°are compatible with the inversion solution.The inversion yielded a source time function of boxcar shape and duration of ~ 8 s.Our solution is in good agreement with the solution published in the Harvard CMT catalogue, although the synthetics produced by the CMT solution (strike, dip and rake fixed, all other parameters free) indicate poor amplitude fitting for some stations like HIA, AAK and KDAK as shown in fig. 3.
Figure 4 summarizes some of the tests that were carried out, in order to investigate the errors of our minimum misfit solution.The fixed values for strike, dip, rake and depth tests are 135°and 170°, 70°and 90°, 20°and 0°and 16 km and 7 km, respectively.

Methodology
In addition to the teleseismic waveform modelling, we applied the time domain moment tensor inversion technique, which uses digital waveform data recorded at regional distances.The employed inversion code has the capability of revealing the seismic moment tensor using a single three-component station, although the use of a larger number of stations provides greater stability (Dreger and Helmberger, 1991).Generally, the methodology can be applied successfully if the following conditions are met: 1) the regional crustal model is sufficiently well known to explain the lowfrequency part of the recorded waveforms; 2) the event location can be well represented by the high-frequency hypocentral location, and 3) the source time history is synchronous for all of the moment tensor elements and may be approximated by a delta function (Dreger, 2000).Furthermore, the inversion procedure accounts only for the deviatoric tensor, neglecting any volumetric changes in the source area.
The applied methodology is based on a simplified general representation of the seismic source by considering a point source in both space and time where U n is the nth observed component of displacement, G ni, j is the nth component of the theoretical Green's function for specific forcecouple orientations, x corresponds to the sourcestation distance, z is the source depth and M ij is the scalar seismic moment tensor.The general force-couples for the deviatoric moment tensor are represented by the three fundamental faults, namely a vertical strike-slip, a vertical dip-slip and a 45° dip-slip (Jost and Herrmann, 1989).
Equation (3.1) is solved using linear least squares, assuming a constant source depth for each inversion.The estimated scalar seismic moment tensor, M ij , is then decomposed into the scalar seismic moment, a double-couple moment tensor and a compensated linear vector dipole moment tensor.The procedure of the decomposition is described in detail in Jost and Herrmann (1989).The optimum hypocentral depth is found iteratively through an examination of both an objective function and the variance reduction.The objective function, f, depends upon the RMS of the difference between the observed waveforms (d ) and the synthetic waveforms (s), modulated by the percent double couple ( pdc): (3.2) while the Variance Reduction (VR), is estimated through the relation Successful applications of the methodology result in small values of the objective function and large values of the variance reduction, indicating that both the waveform fit and the percent double couple are large.

Application and inversion results
In order to apply the time domain moment tensor inversion technique for estimation of the Skyros mainshock focal mechanism, an appropriate 1D velocity model had to be chosen.Unfortunately, due to the geological and seismotectonic complexity of the Aegean area, one single 1D model is not adequate to explain the low-frequency wave propagation along the different wave paths (see fig. 1 and http:// www.gein.noa.gr for more information on the NOA stations).Recently, Zahradnik (2001) tested the model proposed by Novotny et al. (2001), through a forward simulation of the same waveforms used in the present study (regional broadband data of NOA) and found that the examined 1D model can satisfactorily explain the frequency range 0.05-0.08Hz of the waveforms recorded at stations APE, PRK and ATH (fig.1).These stations are among the closest to the mainshock epicentre (their epicentral distances range from 135 to 245 km) and present the larger signal to noise ratios.
We also employed the model of Novotny et al. (2001) and used it as input in the frequencywavenumber integration code (FKRPROG) developed by Saikia (1994) in order to compute synthetic Green's functions for the examined epicentre-station distances.The synthetic Green's functions were then used in the time domain moment tensor inversion code to derive the focal mechanism of the examined event.
We found that the employed velocity model works very well for stations ATH, PRK and APE (variance reduction is larger than 90% for all three stations), although it is inadequate for the rest of the epicentre-station paths.Nevertheless, the applied methodology gives reliable results by using one single station, thus the combination of three stations provides high levels of accuracy in the derived focal mechanism.

VR d s dt
The results from the time domain moment tensor technique are presented in fig. 5 and  table I.The observed tangential, radial and vertical components at stations APE, PRK and ATH (solid lines) are compared to the synthetics (dashed lines).Table I shows that the two applied methodologies gave converging results.

Methodology
Source parameters such as moment (M 0 ), fault length (L), average displacement (u _ ) across the fault and static stress drop ( ), were determined for the Skyros mainshock, using the far-field amplitude displacement spectra.The data consist of Long Period P-waves (sampling frequency 1 Hz) recorded at teleseismic distances (30°-90°) from the GSN stations.We used a time window starting at the P-arrival and ending before the S-wave arrival.The displacement waveforms were corrected for the instrument response, attenuation and radiation pattern.For the latter correction we used an average radiation pattern over the focal sphere, which in the case of the Pwaves is assigned the value R = 0.51 (Fletcher, 1980).The far-field amplitude displacement spectra are characterized by 3 parameters: i) the low frequency level 0 , which is proportional to seismic moment; ii) the corner frequency, f c , and iii) the power of the high frequency asymptote.Following Brune (1970Brune ( , 1971)), we shall define corner frequency at the intersection Fig. 5. Results of the time domain moment tensor inversion.In the left part the comparison between observed (solid line) and synthetic (dashed line) regional broadband waveforms is shown.In the right part the beach-ball, of strike, dip and rake for the two nodal planes, the scalar seismic moment, M 0 , and the moment magnitude, M are also given.The time scale is shown below the vertical component, while information on the station name, examined frequency range, maximum observed amplitude (in cm) and variance reduction is shown below waveforms. of low-and high-frequency asymptotes in the spectrum.
Almost all far-field displacement spectra were characterized by a constant low-frequency level, 0 , and a fall-off above a corner frequency f c , at a rate proportional to f y .Spectra that did not show such a shape were not analysed.Determination of the spectral parameters ( 0 , f c ) was performed be eye fitting low-and highfrequency asymptotes to the observed spectra.Figure 6 shows indicatively the displacement amplitude spectra for 6 stations together with their best fitting results.
The scalar seismic moment was calculated from the relation (Keilis-Borok, 1959) (4.1)where 0 (P) denotes the low-frequency asymptote of the spectrum, the density, R the radiation pattern coefficient for P-waves from a double couple point source, R the epicentral distance and the P-wave velocity at the source.
In order to calculate the source dimensions and the stress drop, a circular fault of radius r was assumed.Following Madariaga (1976), the radius of the source can be estimated from relation where f cp is the corner frequency of the P-wave spectra and is the velocity of the S-waves.We have also assumed that the diameter of the circular fault area is equal to the observed fault length (Hanks and Wyss, 1972).
Stress drop was calculated from the relation of Keilis-Borok (1959) (4.3) and the average displacement, u , was calculated from the relation of Aki (1966) (4.4) where µ is the shear modulus (taken equal to 3.3 * 10 10 N/m 2 ) and A is the fault surface in km 2 .Average values <x> were computed for each parameter (stress drop, fault length, average displacement) following Archuleta et al. (1982): where N is the number of stations used.The basic reason for using this relation is that in the case of simple arithmetic average, the mean values are biased towards the larger values.The corresponding standard deviation of the logarithm SD(log< x >) and the multiplicative error factor, E x , were also calculated from the relations of Archuleta et al. (1982): (4.6) (4.7)

Source parameters of the mainshock from the spectral analysis
The distance, azimuth, the radiation pattern of the P-waves, R (P), the low frequency asymptote, 0 , and the corner frequency, f c , of the stations used are listed in table II.Table III lists the average values and the multiplicative error factor of the stress drop, the radius of the circular fault and the average displacement across the ruptured area.In our calculations, we used 6.5 km/s for the P-wave velocity, 3.7 km/s for the S-wave velocity, 2.6 g/cm 3 for density and 3.3 * 10 10 Nt/m 2 for the shear modulus, µ.The scalar moment, based on eq.(4.1), is M 0 = 8.06.10 18 N m.This value is larger compared to the values obtained by the teleseismic body waveform modelling and the time domain moment tensor inversion.Nevertheless, the previously determined values for the seismic moment of the examined event lie within the range of ± 1 standard deviation of the present estimation from the spectral analysis.

Conclusions
We investigated the source parameters of the Skyros M 6.5 earthquake that occurred in the Northern Aegean Sea using body wave teleseismic waveform modelling, time domain moment tensor inversion using regional waveforms, as well as spectral analysis of the farfield displacement P-wave spectra.We used the first two methodologies, in parallel, both to ensure reliable estimation of the focal mechanism, for further use in more detailed studies of the source process during the Skyros event (e.g., slip distribution studies), and to validate the results of the time domain moment tensor inversion, which is applied for the first time in the Aegean area.The converging results of both methodologies were quite encouraging for further application of the time domain moment tensor inversion scheme to study moderate size Aegean earthquakes.
Figure 7 shows the focal mechanism of the mainshock together with previously determined focal mechanisms of earthquakes of the area (Papazachos et al., 1998 and references therein).In the same figure, the aftershock distribution is also shown for comparison (Melis et al., 2001).It is seen that the NW-SE striking plane is the fault plane which implies sinistral strike-slip motion contrary to what is expected from the strands of the North Anatolian Fault zone that terminates in this area, which is well known for its dextral strike-slip character.This observation of sinistral strike-slip motion in the Northern Aegean Sea is documented for the first time from the study of an earthquake recorded by the modern global and regional networks.The significance of the NW-SE striking structures, inherited from the old compressional phase of the Aegean tectonics will be further discussed in another paper.However, this event gave strong evidence that these structures can be activated from the modern stress field presently acting in the Aegean area.

MFig. 1 .Fig. 2 .
Fig. 1.Regional map showing the sites of the broadband stations, operated by the Geodynamic Institute of the National Observatory of Athens, and the epicentre location (asterisk) of the Skyros July 26, 2001 earthquake.

Fig. 3 .
Fig. 3. Comparison of the «minimum misfit solution» obtained in this study (upper part) to the Harvard CMT solution (lower part).The solutions are not very different but our solution gives a considerably better fitting (HIA and AAK for instance).

Fig. 4 .
Fig. 4. Presentation of the tests performed in order to define the uncertainties in strike, dip, rake and depth (as discussed in the text and listed in tableI).

Fig. 6 .
Fig. 6.Far-field amplitude displacement spectra used to estimate the fault dimensions and other parameters of the Skyros mainshock.The low and high frequency asymptotes (fitted by eye) are depicted as straight lines.The corresponding values of the low-frequency part of the spectrum, 0 , and the corner frequency, f c , are also shown for each spectrum.

Fig. 7 .
Fig. 7. Regional map showing the focal mechanism of the Skyros mainshock, the distribution of the aftershocks, the orientation of the fault that moved during the earthquake, and other focal mechanisms previously determined.The activation of an NW-SE striking fault is evident.

Table I .
Source parameters of the July 26, 2001 Skyros earthquake derived from the two applied methodologies.In the case of the teleseismic waveform modelling, the numbers in the brackets indicate the lower and upper bounds of the uncertainties.

Table II .
Station parameters and spectral parameters obtained from the far-field displacement spectra of P-waves for the July 26, 2001 Skyros earthquake.

Table III .
Mean values and multiplicative error factors of the seismic moment, M 0 , the fault radius, r, the stress drop, and the average displacement, u, as defined by means of the far-field P-wave spectral analysis.