Local earthquake tomography of the Erzincan Basin and the surrounding area in Turkey

In this study, selected travel time data from the aftershock series of the Erzincan earthquake (March, 1992, Ms=6.8) were inverted simultaneously for both hypocenter locations and 3D Vp and Vs structure. The general features of the 3D velocity structure of the upper crust of Erzincan Basin and the surrounding area, one of the most tectonically and seismically active regions in Turkey were investigated. The data used for this purpose were 2215 P-wave and 547 S-wave arrival times from 350 local earthquakes recorded by temporary 15 short-period seismograph stations. Thurber’s simultaneous inversion method (1983) was applied to the arrival time data to obtain a 3D velocity structure, and hypocentral locations. Both 3D heterogeneous P and S wave velocity variations down to 12 km depth were obtained. The acquired tomographic images show that the 3D velocity structure beneath the region is heterogeneous in that low velocity appears throughout the basin and at the southeastern flank, and high velocities occur at south and east of the basin. The low velocities can be related to small and large scale fractures, thus causing rocks to weaken over a long period of the active tectonic faulting process. The ophiolitic rock units mostly occurring around the basin area are the possible reason for the high velocities. The validity of 3D inversion results was tested by performing detailed resolution analysis. The test results confirm the velocity anomalies obtained from inversion. Despite the small number of inverted S-wave arrivals, the obtained 3D S velocity model has similar anomalies with lower resolution than the 3D P-wave velocity model. Better hypocenter locations were calculated using the 3D heterogeneous model obtained from tomographic inversion. Mailing address: Dr. Hüseyin Gökalp, Department of Geophysics, Karadeniz Technical University, 61080 Trabzon, Turkey; e-mail: hgokalp1@hotmail.com


Introduction
The 1992 March 13th Erzincan earthquake (M s =6.8) occurred in the eastern part of the North Anatolian Fault Zone (NAFZ) at the northern margin of the Erzincan Basin ( fig. 1). The basin where the epicenter of the 13th March 1992 event was located is a complex pull-apart basin created by the NW-SE trending right-lateral NAFZ and the NE-SW trending left lateral Ovacik Fault (Barka and Gulen, 1989). It is located between two segments at the NAF that suffered large destructive earthquakes in the past which affected Erzincan city and its surrounding towns and villages. The fault mechanism was dominantly a right lateral strike-slip, striking 122 SE and dipping 65 SW based on the modeling of the broadband P and SH waveforms (Fuenzalida et al., 1997). Some studies reported the focal depth of the main shock to be located in the upper crust as shallow as 10 km (Barka and Eyidogan, 1993;Fuenzalida et al., 1997) while the National Earthquake Information Service (NEIS) and Deschamps and Dufimier (1992) reported it as deep as 25 km.
Following the main shock, the aftershock activity occurring in the region decreased gradually over the following 2 months and the spatial distribution of aftershocks indicated that many secondary faults in the region were activated by the main event (Gürbüz et al., 1993). Two days after the main, a large aftershock (M s = 5.8) occurred about 40 km southeast of Erzincan near the town of Pülümür. Waveform inversion for this large aftershock showed both thrust and left lateral strikeslip components (Fuenzalida et al., 1997).
Eight days after the main shock, a temporary local seismic network (see fig. 1) was installed in the region by the General Directorate of Disaster Affairs (Ankara, Turkey) and Geo-ForschungsZentrum (Postdam, Germany). The travel time data available for the tomographic inversion were collected by the network during the aftershock activity following the 1992 Erz-incan Earthquake. Composite fault-plane solution for the aftershock sequence indicated that a strike-slip mechanism predominated for the aftershocks clustered in the vicinity of Erzincan but in the southeastern part of the Erzincan Basin the fault plane solutions were of normalfault type (Grosser et al., 1998). Some other portable local networks were independently established in the region by the Ecole et Observatoire des Sciences de la Terre de Physique de Globe de Strasbourg, France (Fuenzalida et al., 1997) and Marmara Research Center of Tubitak, Turkey (Aktar et al., 1996). Aktar et al. (2004) determined 3D P-wave velocity model of the Erzincan Basin by inverting P-wave arrival times of the 1992 Erzincan earthquake aftershock sequence. They found that strong lateral variations exist in the upper crust. The thickness of low-velocity soft sediments (Vp<2 km/s) in the Erzincan Basin is ≤ 3 km and Vp values beneath the sediments (>6-9 Erzincan earthquake located by a temporary local seismic network (solid triangles). The focal mechanism for the Erzincan earthquake and the largest aftershock (occurred near the Pülümür town) are taken from Fuenzalida et al. (1997). Hatched area denotes the study region and the grid net used in the inversion. The study area is also shown by rectangles in the small map of Turkey, lower left. km depths) range from 5.8 to 6.3 km/s. Recently, Kaypak and Eyidogan (2005) established the best 1D P-and S-wave velocity models for the basin using arrival times from the aftershock activity.
Tomographic studies of seismogenic structures have offered a robust tool to constrain both the active tectonic structures and the rupture process for several earthquakes (Amato et al., 1990;Eberhart-Phillips, 1990Foxall et al., 1993). This paper presents a detailed 3D P-and S-wave velocity structure of the Erzincan Basin and its surrounding area and an improved image of the aftershock hypocenters distribution relocated in the 3D velocity model. From the evaluation of seismic activity in the region with the 3D velocity structure obtained from tomographic studies, more information about structure of faults or seismogenic zones and its vicinity can be inferred. The results are compared and discussed with those by Aktar et al. (2004) providing new information on the complex structure and the tectonic process in this region.

Tectonic setting
The northward motion of the Arabian Plate towards Eurasia forces the Anatolian block, defined as lithospheric continental plate, to extrude to the west and causes ongoing convergence in Eastern Anatolia and the Caucasus Arc (McKenzie, 1972;Sengör et al., 1985;Dewey et al., 1986). This lateral tectonic escape (Burke and Sengör, 1986) occurs in between two strike-slip faults; the dextral North Anatolia (NAF) and the sinistral East Anatolian (EAF) faults, respectively, which meet at the Karlıova triangle in Eastern Anatolia (Chorowicz et al., 1999). The Erzincan Basin is located between three conjugate fault segments of NAF; left-lateral strike slip Ovacık Fault (OVF) and Northeast Anatolian Fault. The Erzincan Basin is a complex pull-apart basin and it opens as a result of interactions between a segment of the NAF zone and NE-SW trending OVF (Barka and Gülen, 1989). Fuenzalida et al. (1997) pointed out that the Erzincan Basin is wider and more complex than would be expected from a simple pull-apart basin, based on evidence obtained from satellite image and field observations. Barka and Gülen (1989) suggested a two-stage model to describe the basin evaluation. After a «classical» pull-apart opening of the basin (Burchfiel and Steward, 1966) the left-lateral movement along the Ovacik Fault cut the southeastern border fault and rotated the basin clockwise and hence opened it into a «wedge out» model (Grosser et al., 1998). The Erzincan Basin is mostly filled with fluvial Plio-Quaternery deposits with its lateral and transverse dimension being approximately 50 km and 15 km, respectively. Although the thickness of the basin sediments is unknown, Hempton and Dunne (1984) give an empirical relationship between the length and the sediment thickness of a basin so the thickness of the basin sediments has been estimated to be up to 3 km (Grosser et al., 1998). Several small volcanoes, which formed during the pull-apart opening, are aligned along the margins of the basin, generally along the northeastern margin (Barka and Gülen, 1989). These volcanoes and hot springs indicate high heat flow and thinning of the crust (Aydın and Nur, 1982).

Data
The aftershock data were collected by temporary 15-station network installed by Frankfurt University, the GeoForschungsZentrum, Postdam and the General Directorate of Disaster Affairs, Ankara. The seismic stations were equipped with short period vertical component seismometers. Array dimension is approximately 80×40 km. They recorded analogue data continuously on magnetic tapes that were digitized at Frankfurt University with a sampling rate of 128 Hz (Grosser et al., 1998). Some stations were moved to new locations and due to logistic problems, were operated at different time intervals at the same locations. For this reason, the data referred to in this paper are the aftershock data recorded by 15 stations, between 21/03/1992-16/06/1992. The code HYPOINVERSE (Klein, 1989) was used to determine locations for a total of 505 aftershocks. The majority of events clustered at the southeastern end of the Erzincan Basin with only a few events located close to the fault plane of the main event at the northeastern margin of the basin. Most of the aftershocks had a source depth between 5 km and 11 km with the maximum at 6.5 km depth.
In order to obtain the image of the crustal structure beneath the Erzincan Basin, 350 earthquakes located in the network area were selected with at least 7 P and 7 S arrival times and with horizontal and vertical standard errors ≤2 km for most of the event. The data were also weighted according to their accuracies inferred from the signal to noise ratio of the seismic signal and for arrival time residuals and hypocentral distance. Consequently, 2215 P and 547 S arrivals were selected for the inversion. A 1D crustal model by Grosser et al. (1998) was used for the preliminary location of the events. The velocities and depths of discontinuities in the model were V P1 =5.3 km/s, H 1 = 0.0 km, V P2 = 6.0 km/sn, H 2 = 4.0 km, V P3 = 8.0 km/sn, H 3 =50.0 km and V P /V S =1.78.

Local tomography technique
In local earthquake tomography, arrival times from local earthquakes are inverted to estimate velocity and hypocentral parameters (origin time and location) simultaneously. The SIMULPS code originally written by Thurber (1983) and subsequently modified by Eberhart-Phillips (1986, 1990 was used for 3D inversion. The study volume is parameterized in terms of velocities at the nodes of a 3D grid and linear velocity gradients are assumed between nodes. The velocity at any arbitrary point in the 3D velocity model is calculated by linear interpolation between the surrounding nodes. Starting with a priori values, the program performs a damped, iterative, least squares inversion for hypocentral parameters and the 3D velocity structure. Damping is used to suppress large perturbations in the calculated parameters. Thurber's method has been successfully used in many areas all over the world and is described in detail by Thurber (1983Thurber ( , 1993 and Eberhart-Phillips (1990. To solve the forward problem, approximate 3D ray tracing and pseudo bending (Um and Thurber (1987) were applied. Hypocenter locations were updated with the new velocity model at each iteration step. The damped least squares solution to the linearized problem is given by where M is the vector of model perturbation, d is the vector of arrival time residuals, G is a matrix of partial derivatives, θ 2 is the damping parameter and I identity matrix. Once the model perturbations are calculated and applied to the existing model, the earthquakes are located using a weighted least squares method. Pavlis and Booker (1980) have shown that this is mathematically equivalent to solving the simultaneous problem. In essence, since the problem is to solve two different sets of parameters (velocity model and hypocentral adjustments) simultaneously, the null data set resulting form an orthogonal transformation is created and then used for the velocity model inversion. The process is repeated until there is no statistically significant improvement to the model (measured using an F test on the rms residuals). An appropriated value of θ is selected, optimizing the data variance reduction without an overly complex velocity model.
The resolution of the resulting images and the standard errors are also estimated by the inversion procedure. The model resolution is a measure of the fit between the obtained model and «true» model. The model resolution R for a damped least squares problem is (Menke, 1989) (4.2) The diagonal elements of R represent the resolution of each model parameter. Values close to 1 indicate well-resolved parameters, whereas lower values show us that the computed velocities are smeared over a larger model volume. The model standard error provides an estimate of how the data error is mapped into the model error. Menke (1989) showed how the model standard error, S, is related to the data error σ d as follows: where N is a normalization factor that accounts for the volume influenced by α n , which is the influence coefficient used in the linear interpolation and depends on coordinate position. It is defined as (4.7) A large DWS indicates that the velocity at the grid point is based on a large body of data. The DWS is a superior indicator of ray distribution than the «hit matrix» since it is affected by the distance of the rays from the velocity node. Empirical studies indicate that a DWS of 50 or greater corresponds to a compact averaging vector with a spread of up to about 2 (Toomey and Foulger, 1989).

1D model
Choosing a suitable initial model helps the convergence toward the global minimum in the 3D inversion, strongly constraining the 3D velocity results. One of the most powerful approaches is to improve the inversion with the available independent knowledge of local structures (Eberhart-Phillips, 1990). Using a simple initial model in 3D inversion is better than using a complex initial model because it has the advantage that solution features are those required by the data and not due to the peculiarities of the initial model (Eberhart-Phillips, 1990).
Yusufhadjaev and Bayraktutan (1994) (i.e. YB) provided a 1D multi-layer velocity model to represent the (shallow) crustal structure beneath the Erzincan area, which was chosen as a starting point. The YB model was discretized in the vertical direction by assigning a velocity node in the middle of each layer and was used as an initial model in the 1D inversion of the travel time observations, performed using the same SIMULPS code used in the 3D inversions. The YB model did not yield the best results on the RMS of the The model covariance and therefore the standard error are directly related to the model resolution. High resolution, which is desirable, usually results in high standard error, which is undesirable. The trade-off between model resolution and standard error must be taken into account when choosing a damping value. The analysis of the resolution matrix is given by the spread function (Backus and Gilbert, 1967;Menke, 1989) S p , representing the spread of R relative to the diagonal element for each parameter, and the averaging vectors, indicating the volumes over the velocities which are averaged. A well-resolved node presents a compact vector (a small volume centered in the node) and its spread function has a low value given by (4.5) where r p is the averaging vector at the pth parameter, R pq is an element of the resolution matrix, and w(p, q) is the weighting function that corresponds to the linear distance between nodes in the grid (Menke, 1989).
A careful analysis of reliability must be conducted before any interpretation of tomographic results can be made. In fact, any tomographic image is only as good as its resolution estimates. When assessing the quality of the tomographic inversion, one must determine how well each node was illuminated (number of rays hitting the node) and how well it was resolved.
The very rough estimate of illumination of the model contained in the hit count matrix (reporting the number of rays that contribute to the solution of that node) is supported by the Derivative Weight Sum (DWS), a more sensitive measurement of the spatial sampling of the model space. It quantifies the relative ray density in the volume of influence at a model node, weighting the importance of each ray segment by its distance from the model node (Haslinger et al., 1999). In order to assess the compactness of the averaging vectors of R, the Derivative Weighted Sum (DWS) (Toomey and Foulger, 1989) is used defined as follows: travel-time residuals and the data variance reductions. Therefore, several 1D models were obtained by manually perturbing the layer thicknesses and layer velocities of the YB model and the above inversion procedure was repeated for each of them. The latter process resulted in many 1D inverted models among which the best one was chosen on the basis of data fitting performance, i.e. the RMS of the travel-time residuals and the data variance reductions. The corresponding 1D optimal model parameters are listed in table I. This model was utilized in the 3D tomographic inversion as an initial model. The damping parameters for the P and S arrival times, which were due to the trade-off analysis between the data and model variances, were set to 150 and 175, respectively. In the 1D velocity analysis, the V p /V s ratio was simply set to 1.78, i.e. Poisson's ratio was assumed 0.27. The choice of suitable damping parameter was made by constructing a trade-off curve and selecting the value that most reduced the data variance without causing a large increase in the solution variance (Eberhart-Phillips, 1986). The trade-off curve was computed using a large range of damping values, performing one iteration inversions.

3D velocity model with resolution analysis
3D inversions were performed mainly with four-grid spacings (20 km, 15 km, 10 km, 7.5 km) to obtain the information from the data with different scales concerning the velocity structure beneath the study region and to reach the best compromise between the spatial definitions of velocity images and resolution of model parameters. A trade-off analysis was also performed for the 3D inversion of the data. The trade-off curves were constructed for each model for 3D inversion to choose the suitable damping parameter. Beside the four main initial models, the data were also inverted with different initial velocity model, which has different grid spacing in the 3D grid mesh, and with a different damping parameters so as to compare the results and to verify the stability of the final model. The value of chosen damping parameter changes during the inversion steps.
The results of the inversion with a grid spacing of 10 km are given. This has the best spatial definition of the velocity structure beneath the study area, due to the fact that the approximate station spacing of the network is 10 km. From the analysis of the resolution parameters, this model also has better resolution than the others. To select the optimal value of damping factor, «trade-off» curves were constructed by performing the inversion with damping factor values ranging from 0.1 to 500 ( fig. 2). The optimal value substantially increases the resolution without introducing large velocity errors. Finally, the values of 200 and 50 were selected for the inversion of the P-and S-wave arrival times, respectively, which allows a good resolution with errors lower than 0.1 km/s especially for P-wave velocities. The value of the damping parameter is modified after each iteration step, according to the data variance, and generally increases (see Eberhart-Phillips, 1990). From the results of the inversion, the number of inverted parameters is 352, the final rms of the 3D inversion is 0.11386 s and the variance improvement with respect to the starting model is 33% and 9% for the P and S data, respectively. Figure 3a,b shows the computed 3D P-wave velocity model with the layers located at 0, 2, 3, 6, 8 and 12 km depths. In order to represent shallower crustal structure which is probably more heterogeneous than the deeper layers, a smaller vertical node spacing was chosen in the 3D model. Figures 4 and 5 show the S-N longitude and W-E latitude verti-  . 6a,b) has similar node distribution in 3D with the P-wave model, the nodes in the layer at 12 km depth are fixed. Despite both models having grid nodes in the layers of 12, 15, 20 km depths, the nodes in these layers remained fixed during the inversion due to poor ray sampling. S-wave velocities were set suitably using a V p /V s ratio of 1.78. Figure 6a,b shows the final S velocity model in five slices at different depths (0, 2, 3, 6 and 8 km). From a general analysis of both P and S velocity models, the presence of complex velocity variations beneath the region can be observed. Velocity images show similar anomaly patterns from the shallow to deep layers. Lowvelocity anomalies are present beneath the Erz-incan Basin, from NW to SE. A low velocity anomaly is more dominant in the SE of the basin, in the vicinity of the Ovacık fault. The station coverage and the resolution in the model layers are good for this area. Two high velocity anomalies are present at the southern (or southestern) edge of the basin. High velocity anomalies are also observed at the eastern of the basin, beneath Çayirli town. This pattern can be seen from the first layer to the layer at 8 km depth. No large velocity contrast is observed in the layer at 12 km depth; only a slight low velocity anomaly is shifted to the center of the basin. The same anomaly pattern was revealed with different grid spacings. The main features of the velocity pattern in these models are a NW-SE trending low velocity anomaly beneath the Erzincan Basin and high velocity anomalies surrounding the central low velocity anomaly at the southern and eastern flanks of the basin. No velocity contrast was observed in  the layers deeper than 12 km, probably due to poor ray sampling. The diagonal element of the resolution matrix R is the best and quickest estimate of the quality of the solution. It shows the degree of independence of the model parameter in the solution by showing the diagonal element of the full solution matrix; therefore, it is a diagnostic tool which ranges from 0.0 to 1.0. There is a major problem of how to use the diagonal element, DWS, spread function as diagnostic tools once they have been computed. One should fix reliability limits that tend to have good sampling, high value of diagonal element and DWS   Fig. 6a,b. a) Horizontal percentage S-wave velocity change relative to the 1D initial model. The grid spacing is 10 km in a horizontal direction. b) ∆Vp (%) (i.e. percent change of P velocity relative to 1D initial model) are shown in two vertical cross sections along the lines A-B and C-D in the first tomographic section (the layer at 0 km) for the model. two parameters through the model and making a comparison between them, it was decided to consider all the nodes with DWS ≥ 2000 (diagonal element ≥ 0.2) as fairly well resolved. From the analysis of the three resolution parameters, it can be concluded that the velocity anomalies obtained from the 3D inversion are well resolved. The central part of the Erzincan Basin is densely visited by crossing ray paths and therefore is better resolved. The ray path density is relatively poor near the edges resulting poor resolution there (fig. 7). The layers at 0, 2 and 3 km depths are less illuminated than the deeper layers, due to smaller vertical node spacing. It should also be noted that the diagonal elements of the resolution matrix are also lower for the nodes in the shallower layers but the spread function values for these nodes are slightly larger.
Despite the small quantity of S data, 3D Swave velocity models were determined from the inversion of S data ( fig. 6a,b) and relevant Swave velocity anomalies were obtained in the deeper layers. Figures 9 and 10 show the variations of the DWS and spread function, respectively, in the obtained 3D S velocity models. From the analysis of the two parameters, it can be concluded that the 3D S velocity model generally has lower resolution than the 3D P velocity model as expected due to the smaller number of S data with a greater quantity of uncertainties. For this reason, the 3D S velocity model is less resolved than the 3D P-wave velocity model.  The resolution analysis indicates that the best resolved region corresponds to the low velocity anomaly in the southeastern of the basin. The other relevant low velocity anomaly is located at the center of the basin where the 1992 March 13th Erzincan earthquake occurred, in the vicinity of Erzincan town, extending between 6 km and 8 km depth. A high velocity anomaly is found at the east of the basin, in the vicinity of Çayırlı town at 6 km depth. It is shifted towards the southeast of the basin where the low velocity anomaly exists in the layer at 6 km depth.
The 3D S velocity model is generally in agreement with the 3D P velocity model. For both models, the same anomaly pattern occurred at the same location indicating that there is a low velocity zone striking NW-SE through the basin. The low velocity anomaly located at the southeast of the basin is the best resolved in both the 3D P and S velocity models. These results show that both models are consistent. While a high velocity anomaly occurred in the south of the basin in the 3D P-wave velocity model, it is not observed in the S velocity model at same location probably due to both small quantity of the S-wave data and lack of ray paths crossing beneath this area.
The other obtained 3D P-wave velocity models with grid spacings of 20 km, 15 km and 7.5 km ( fig. 11) give the same high and low velocity pattern with various anomaly fluctuations at the same location as the model with grid spacing of  10 km. This result confirms that the anomaly patterns recovered in the P-wave velocity model with a grid spacing of 10 km is consistent and the model is stable. Furthermore, in the model with grid spacing 7.5 km, the layers were located at 1, 3, 6, 8, 12 and 16 km depths for the investigation of the target volume with a different vertical spa-tial scale. On the other hand, the obtained model with grid spacing of 10 km has slightly better resolution due to having lower spread function values and has better spatial definition due to having higher DWS values. From the evolution of DWS distributions for the models, the spatial resolution decreases with smaller grid spacing in the models Fig. 11. Horizontal percentage P-wave velocity change relative to the 1D initial model. The grid spacing is 7.5 km in a horizontal direction. and the model with a grid spacing of 7.5 km has the lowest spatial resolution. The results with the grid spacing of 7.5 km show that same anomaly patterns were recovered at the same locations as in the results with the grid spacing 10 km with more anomaly fluctuations. The observed main anomalies located at the center and southeast of the basin in all models were recovered in good resolution corresponding to the well-resolved areas of the models.

Earthquake location
The initial HYPOINVERSE (Klein, 1989) locations were improved in two stages. First, the whole arrival time data was relocated using the 1D model based on early simultaneous inversion results. Second, the events were relocated using the 3D velocity model obtained from the simultaneous inversion (fig. 2). The seismicity relocated in the 3D velocity model shows a better defined aftershock distribution concentrated in a more compact volume. Figure 12 shows the comparison between relocations of the aftershocks using the final 3D velocity model and 1D velocity model. The most intensive local seismicity is observed beneath the southeast of the Erzincan Basin where the low velocity anomaly was observed in both the 3D P and S velocity model. This low velocity volume within which most of the aftershocks occurred in the depth range of 5 to 11 km probably corresponds to the nucleation of the Pülümür main shock. Overall, the 1D to 3D hypocentral shifts are modest but reduce the scattering. The largest amount of hypocentral shifts are about 2 km in the northwestern part of the basin, whereas the hypocentral shifts are about 500 m in the southeast of the basin where intensive aftershock activity occurred.

Discussion and conclusions
The recorded aftershock sequence of the 1992 Erzincan earthquake by a temporary network provided the study of the seismogenic structure of Erzincan Basin and the surrounding area. By inverting aftershock data using Thurber's simultaneous inversion method (1983), a detailed 3D P and S velocity structure in the Erzincan Basin was revealed and also accurate images of the location and shape of this seismically active zone were obtained. Four initial models with different grid spacing (20 km, 15 km, 10 km, and 7.5 km) were used in the inversion to both recover all the avail- able information from the data in different spatial scales and to check the stability of the solution of the inversion. It can be concluded that the inversion results for all the initial models gave similar anomaly patterns. This comparison shows that the obtained inversion solution is stable. The 3D velocity model with a grid spacing of 10 km is consistent with the approximate station spacing of the network so it is the best model presenting velocity structures beneath the region. It has better resolution from the evaluation and comparison of resolution parameters for all models.
In this study, both the upper crustal structure within the Erzincan Basin was modeled to a depth of 12 km and the alignment of the aftershocks was determined. The obtained tomographic images show the existence of a complex velocity structure beneath the basin. Low velocity anomalies are dominant both within the basin and southeast of the basin where the second main shock (Pülümür earthquake) occurred. This imaged zone with low velocity anomaly is low velocity corridor which extends in the NW-SE direction parallel to the main fault (NAF) and to a depth of 12 km. The high velocity anomalies were located in the east and south of the Erzincan Basin extending to a depth of 8 km. Resolution analysis confirmed the reliability of the overall features of the imaged velocity distribution beneath the Erzincan Basin. On the other hand, the obtained 3D S velocity model generally shows a similar velocity anomaly pattern to the P velocity model. Due to lower amount of the S data with more ambiguous phase reading in the inversion, the calculated 3D S velocity model has a lower resolution than the 3D P velocity model as expected.
The 3D P and S velocity models reflect the 3D geologic and tectonic structure in the Erzincan Basin. One of the two main local low-velocity anomalies through the basin in the tomographic images seemed to be located beneath the epicenter of the Erzincan earthquake (March 13, 1992) and probably corresponded to main shock nucleation. Likewise, the other low velocity anomaly in the southeast of the Erzincan Basin was located at the beneath epicenter of the largest aftershock (March 15th, 1992) that occurred near Pülümür town and probably corresponded to its nucleation. Aftershock areas were imaged as low velocity anomalies. Generally, the acquired low velocity areas are considered to be associated with the weak materials (or a weakness zone) fractured by the active faults acting in the region for a long time. Furthermore, the relocated hypocenters based on the 3D model are mostly located at the southeast of the basin and confined at a depth range of 2 to 12 km. All these findings show the existence of an active fault zone in the region where the strength rocks deformed by the faults acting in the region are, causing a weakness zone.
Although the use of S waves in the inversion is generally limited mostly due to lack of the three component data, it can provide information on structure as in the earthquake source region that P waves alone do not. The lack of three-component readings results in less well-distributed S data and the inversion can produce only negligible velocity changes where there is low resolution (Thurber, 1993). The resolution for S-wave velocity model solutions is quite different than for P velocity models. One reason for the difference in resolution pattern is having a different set of stations besides using a smaller number S data in the inversion. While a low-or a high-velocity zone exists in the P-wave velocity model, it cannot be found in the V s solution. Therefore while the Pand S-wave velocity models are similar, we can not put much importance on different shapes and locations of specific velocity features in the P and S solutions (Eberhart-Phillips, 1993). In this study, P-and S-wave velocity models show an identical anomaly pattern. Two local low-velocity anomalies through the basin and a local highvelocity east of the basin where the best resolved regions are, determined in the both P-and S-wave tomographic images. Due to ray geometry between hypocenters and stations, most S waves had to travel mainly through the cracked weakness zone, therefore the travel-time difference between P-S waves arrivals increased with clearer S phase identifications, corresponding to the S data used in the inversion, have been performed. This is a possible explanation for obtaining similar anomaly pattern in both P and S velocity models.
A high velocity anomaly located south of the basin can be related to the carbonatic and ophiolitic units in the region. It most likely corresponds to the ophiolitic melange and the Upper Triassic to Lower Cretaceous limestone unit of Munzur Mountain (Özgül and Tursucu, 1984;Westaway and Arger, 2001) which is a deep rooted structure formed at the southern border of the basin. Likewise, the other high velocity anomaly that occurred at the east of the basin can be interpreted as metamorphic and ophiolitic units.
It is worth noting that the shape of the basin is imaged by the low velocity zone as a boat shape extending in the NW-SE direction at shallow depths between the surface and 8 km depth, assuming relative low velocity variations in the central area which are well resolved in the tomographic images corresponding to sediment filling the basin and fractured rock units. This low velocity zone begins to disappear at a depth of 8 km. This depth probably corresponds to the bottom of the basin and to a transition point from basin sedimentary and fractured fault zone units with low velocity to bedrock with high velocity. These results are consistent with the P-wave velocity model obtained by the results of the tomographic inversion given by Aktar et al. (2004). The overall patterns of the P-wave images are similar, but there are some differences in the details. In this study, the basin was determined to have a low velocity zone extending in the NW-SE directions with a deeper depth range 0-12 km than the tomographic results of Aktar et al. (2004). Two high velocity anomalies imaged by this study were located at the same place as the tomographic images of Aktar et al. (2004), however, there is no indication of velocity variation due to the existing small volcanoes aligned along the northeastern margin as in the results of Aktar et al. (2004) since the resolving power of the inversion too lower to detect such a relative small size volcanic structures taking place in the region. In this study, beside the 3D P-wave velocity model, a 3D S-wave velocity model is also presented, whereas the result of the S data inversion does not exist in the study of Aktar et al. (2004). The presented results confirmed those by Aktar et al. (2004) but in a smaller region due to lower resolving power.
A detailed analysis of the resolution throughout the model has guaranteed image fidelity, allowing a definition of the region of the model where the computed anomalies are meaningful. Since high velocity anomalies are more dominant in the deeper layers (at depths of 8 and 12 km) the shape of the basin was discovered to be deeper on the northwest and narrower east of the basin.
The improved hypocenter locations show that most of the aftershock sequences occurred southeast of the Erzincan Basin in the depth range 5-11 km. This probably corresponds to the nucleation of the Pülümür main shock.