Crustal Structure of the Eastern Anatolia Region (Turkey) Based on Seismic Tomography

: Here, we investigated the crustal structure beneath eastern Anatolia, an area of high seis-micity and critical significance for earthquake hazards in Turkey. The study was based on the local tomography method using data from earthquakes that occurred in the study area provided by the Turkiye Cumhuriyeti Ministry of Interior Disaster and Emergency Management Directorate Earthquake Department Directorate of Turkey. The dataset used for tomography included the travel times of 54,713 P-waves and 38,863 S-waves from 6355 seismic events. The distributions of the resulting seismic velocities (Vp, Vs) down to a depth of 60 km demonstrate significant anomalies associated with the major geologic and tectonic features of the region. The Arabian plate was revealed as a high-velocity anomaly, and the low-velocity patterns north of the Bitlis suture are mostly associated with eastern Anatolia. The upper crust of eastern Anatolia was associated with a ~10 km thick high-velocity anomaly; the lower crust is revealed as a wedge-shaped low-velocity anomaly. This kind of seismic structure under eastern Anatolia corresponded to the hypothesized existence of a lithospheric window beneath this collision zone, through which hot material of the asthenosphere rises. Thus, the presented results help to clarify the deep structure under eastern Anatolia.


Introduction
Active tectonic forces at the boundaries of lithospheric plates lead to strong seismic activity in tectonically active regions that affects the infrastructure and lives of people therein. In many cases, the plate interaction mechanisms, which are important for assessing long-term seismic hazards, are not completely clear. Promptly addressing problems related to collisional processes requires robust information about the main structural elements in the crust and upper mantle. For this reason, geophysical investigations, particularly seismic tomography studies, are in high demand for establishing geodynamic scenarios of lithospheric evolution in many regions of the world.
We implemented seismic tomography to study the crustal structure beneath eastern Anatolia, including the eastern Anatolian accretionary complex and the surrounding areas. Figure 1 shows the main tectonic features of the study region. The East Anatolian Plateau is bounded by two Neo-Tethyan sutures: the İzmir-Ankara-Erzincan suture (IAES, shown in Figure 1) in the north and the Bitlis suture in the south [1]. Moreover, the IAES separates the Pontides Unit (PU) in the north from the East Anatolian Accretionary Complex (EAAC) to the south [2] (Figure 1). The eastern part of Turkey is characterized by the junction of three large plates: the Anatolian plate in the west, the Arabian plate in the south, and the Iranian plate in the east. This system is divided by tectonic faults, in-cluding the North Anatolian fault (NAF) and the East Anatolian fault (EAF), which intersect in the study area and form the Karliova triple junction (KTJ) (Figure 1). This system of plates is pushed northward by the movement of the Arabian plate at a rate of ~18 mm/yr with respect to the Eurasian plate [3,4]. This movement causes westward displacement of the Anatolian block along the right-lateral strike-slip NAF at a rate of ~24 mm/year. To the south, the Anatolian plate is bounded by the EAF, which exhibits left-lateral displacement at a rate of ~9 mm/yr [4]. To the west of the NAF is the central Anatolian fault zone (CAFZ), which is a transtensional left-lateral strike-slip fault zone extending from central Anatolia to the Mediterranean Sea [5]. The region of eastern Anatolia features complex tectonic structures, which formed during its multistage evolutionary history. The general shortening regime of eastern Anatolia is mostly controlled by the Arabian-Eurasian collision, which began in the mid-Miocene [1,[7][8][9]. The continental collision between the Arabian and Eurasian plates began after the closure of the Tethys Ocean, which separated Gondwana and Laurasia. The formation of marginal basins and the accretion of island arcs and other continental fragments, as well as the subduction of intervening portions of the oceanic lithosphere, took place during the closure of the Neo-Tethys. The East Anatolian Accretionary Complex (EAAC, shown in Figure 1) was formed during subduction from the Late Cretaceous to the earliest Oligocene and included Neo-Tethyan material. The EAAC, having formed on the northward-subducting oceanic lithosphere, can be considered as a remnant of a large accretionary prism situated between the Pontides and the Bitlis-Portuge massif. There is no subducting slab beneath the EAAC.; This could be due to the detachment of the plate under the collision zone. [10] argues that the slab, whose subduction was generating the Pontide arc in the north, was not attached to the Arabian plate. Instead, the slab was possibly attached to the Bitlis-Poturge block before the break-off. The Bitlis-Portuge massif (BPM shown in Figure 1) consists of medium-to-highly metamorphosed units [11].
Several petrological studies suggest that subduction was the likely cause of volcanic activity in the region [10,12]. Although this issue is still actively discussed in the scientific community, most argue that the large volume of melt generated under eastern Anatolia was caused by hot asthenospheric upwelling. Most of the area of eastern Anatolia is covered by the products of recent volcanism (with ages from 20 Ma to the present day). According to the results of isotopic dating, the younger magmatism (6-4 Ma) occurred around Lake Van, whereas the older magmatism is mostly concentrated in the northern part of the plateau (eastward part of the study area). This suggests that the heating and melting zones migrated southward during the past 20 Ma [13,14].
In the high-elevation East Anatolian Plateau, an uplift of ~2 km is observed relative to the adjacent parts of the Arabian plate to the south (e.g., [9,15]). According to the results of receiver function analysis and surface wave tomography studies, the crust of the Anatolian plate gradually thickens to the east. The deepest Moho boundary is observed beneath eastern Anatolia (~50 km) [16]; beneath central Anatolia, the crustal thickness is ~ 37-47 km, and beneath western Anatolia, the crust is ~30 km thick [16,17]. According to the results of [18], in the latitudinal direction, the crust thickens northward from 42 km near the southern part of the Bitlis suture zone to 50 km along the northern edge of the Anatolian plateau in the vicinity of the North Anatolian fault. The thickness of the crust on the Arabian plate is ~40 km [19].
Some geodynamic models (see, for example, [11,[20][21][22]) suggest that the thickening of the eastern Anatolia crust and its uplift are associated with the delamination of the mantle lithosphere within the Arabian plate. Several regional seismic tomography models consistently show traces of slab detachment/break-off and asthenospheric upwelling beneath the Arabia-Eurasia convergence zone (see, for example, [23][24][25][26][27][28][29][30]). In all of the above studies, high-velocity anomalies beneath the margins of the Arabian and Eurasian plates indicate the subsidence of delaminated material, whereas low-velocity anomalies are observed directly under the collision zone, which can be explained by the absence of strong thinning of the mantle lithosphere and the filling of this gap with hot asthenospheric material.
As shown above, a large number of studies have been performed on the regional structures of Turkey, and the geodynamic processes in these areas were successfully investigated. However, despite the abundance of publications, the crustal structure of eastern Turkey is far from being perfectly understood, as many of the derived models demonstrate considerable inconsistencies with each other. In this study, we present new seismic tomography models of the crust and uppermost mantle (down to ~60 km), which were built using a new large set of travel time data. These results can help determine the links between the surface geology and deep structures controlled by complex collisional processes.

Data and Method
The models of seismic velocities were obtained using the arrival times of P-and Sbody waves from local events in the eastern Anatolia region recorded by stations in Turkey and the surrounding countries. The dataset was provided by the Turkiye Cumhuriyeti Ministry of Interior Disaster and Emergency Management Directorate Earthquake Department Directorate of Turkey. The magnitudes of most of the earthquakes in the study were generally between 1 and 7, and their focal depths were mostly between 5 and 35 km.
The dataset included the travel times of 114,800 P-waves and 88,962 S-waves from 13,029 earthquakes were recorded by 187 seismic stations.
To perform the tomographic inversion, we discarded data with fewer than a total of eight P and S picks per event. At the stage of determining the preliminary source locations in the initial 1D velocity model, we removed the data with residuals exceeding 1 s. After applying these two criteria, for the tomographic inversion, we obtained a dataset with 54,713 P-and 38,863 S-wave picks from 6355 events (almost 15 picks per event on average). Figure 2 demonstrates the coverage of the rays from this dataset. The inversion of the P-and S-waves arrival times were performed to obtain 3D velocity models and the locations of the events by using the nonlinear local earthquake tomography code LOTOS (Local Tomography Software) [42]. This code has been used to study various collision zones with similar data distribution geometries, such as the Tien Shan [43,44], the Caucasus [24], India [45], Calabria [46], and others.
The calculations began with locating the sources in the reference 1D velocity model using the grid search method. At this step, to make the calculations faster, the travel times of seismic waves were computed using a table of travel times calculated at a preliminary stage for all possible depths of events and epicentral distances.
Then, we performed several iterations of the tomographic inversion procedure, which included the relocation of sources, the calculation of the sensitivity matrix, and the inversion. The source relocation procedure in the first iteration was performed with the reference 1D velocity model and then with the updated 3D velocity models calculated in the preceding iterations. The search for the optimal source coordinates was based on the gradient optimization method (see [42]), which more rapidly calculates the most likely source locations. The algorithm of a 3D seismic ray tracing, based on the bending method proposed by [47], was used to calculate the travel times within the iteration.
The velocity models were parameterized using a set of nodes deployed according to the ray density. In map view, the nodes were installed regularly with a spacing of 25 km. In the vertical direction, the grid spacing was inversely proportional to the ray density, but the spacing could not be smaller than a predefined value of 10 km. The results were calculated using four grids (with orientations of 0°, 22°, 45°, and 67°) and then averaged them in one 3D model. The inversion was performed simultaneously for the P-and Swave velocities, as well as for the source parameters (dx, dy, dz, and dt) and for the station corrections. For the inversion, we used the LSQR algorithm by [48] and [49]. The inversion stability was controlled by damping and flattening. This procedure minimized the differences in the solutions from neighboring nodes. To derive the final models, we used five iterations in total.
The optimal initial 1D model was estimated after several runs of the full inversion. In the first run, we used a starting model defined as an average model derived from previous studies in eastern Turkey (e.g., [50]). Then, to determine the optimal reference model, we computed the average values of the P-and S-wave velocities at certain depths were computed after the full inversion procedure. These values were used as a new reference model for the next run. Thus, after several runs, we obtained the optimal reference model presented in Table 1. The P and S velocity values in the 1D model were defined at few depth levels. The velocity changes were linear between those depths. We avoided sharp interfaces, such as the Moho, to define artificial horizontal structures to obtain the final velocity model.

Seismic Tomography Models and Tests
The main local tomography results are the 3D models of P-and S-wave velocity anomalies, which are presented in horizontal slices at depths of 10, 30, and 50 km in Figure  3 and in four vertical sections passing through the whole study area and crossing the Bitlis suture zone in Figure 4.   Figure 3 shows that the derived P-and S-wave velocity structures correlated fairly well with each other, which may serve as an informal verification criterion. Indeed, at a regional scale, most of the geological structures (for example, sedimentary basins and igneous provinces) were similarly expressed by the P-and S-wave velocities. The vertical sections of Vp/Vs ratio models, shown in Figure 4, demonstrated similar structures as in the vertical sections of the main resulting dVp models. The fact that the independently derived dVp and dVs highlighted the same structures, the existing data provide sufficient coverage to reliably recover the structures discussed.
To assess the role of random noise in the data, we performed the so-called odd/even test, which consists of independent inversions of two data subsets distinguished by a random criterion, for example, using odd-and even-numbered events. Then, we compared the inversion results obtained using the same parameters and the same algorithms. If the images were strongly different, it means that the random factor in the data was essential. The results of the odd/even test in our case were demonstrated for both the P-wave and the S-wave velocity anomalies at a depth of 30 km in Figure 5. We could observe a good correlation between the major recovered structures. The main features of the velocity heterogeneities, which are important for our interpretation, could be traced in all the models. Synthetic modeling is an important step in tomography, as it provides estimates of the spatial resolution and reliability of the recovered model and allows the optimal values of free parameters (grid spacing, damping, and smoothing) to be determined. When performing synthetic modeling, we simulated the conditions of actual tomographic studies. Synthetic travel times were calculated for the existing source-receiver pairs through a predefined 3D synthetic model using the bending ray tracing algorithm and were then perturbed with random noise. To start recovering the model, we started processing with the absolute source determination results with the initial 1D velocity model, similar to the experimental data inversion. All the steps and controlling parameters during the recovery procedure were identical to the case of computing the main tomographic model based on the experimental data.
To assess the horizontal resolution of the model, we performed a series of checkerboard tests with anomalies of different sizes, as presented in Figure 6. We preset a horizontal checkerboard with rectangular positive and negative anomalies with lateral dimensions of 40 × 40 km, 60 × 60 km, and 80 × 80 km, and with amplitudes of ± 6% (Figure 6). In this particular case, the anomalies remained unchanged with depth. As shown in Figure  6, the checkerboard structure was robustly resolved for both the P-wave and the S-wave velocity models. For shallower depths (up to ~10-15 km), even the 40 km size pattern was robustly resolved in most parts of the study area. At a depth of ~30 km, the anomalies at 60 × 60 km were more resolvable than those at a depth of 40 km. Additionally, in the 60 km depth section, the resolvable pattern size was approximately 60 km. However, the models with smaller anomalies were strongly smeared at this depth. Based on this test, we could conclude that for most parts of the study area, we had sufficient resolution to reveal the main structures, which will be discussed later. The parts of the model where the checkerboard was not correctly restored should be considered with prudence in the interpretation.
To estimate the vertical resolution, we performed a synthetic test with anomalies defined along vertical profiles. Because of the trade-off between the source depth and velocity distribution, the vertical resolution is typically much poorer than the horizontal resolution in cases of passive source tomography studies. In Figure 7, we present vertical checkerboard tests associated with the vertical profiles used to present the main results. Across these sections, the anomalies were defined within a 100 km thick band. In the horizontal direction, the anomalies had a size of 90 km. In the vertical direction, the signs of the anomalies changed at depths of 20 km and 70 km. For the shallow part of the model, the recovery result of this test appeared to be robust; the vertical resolution was sufficient to detect the change in the anomaly sign at a depth of 20 km depth, but the model did not have a sufficient resolution to recognize anomaly sign changes at 70 km. It may also be observed that the anomalies at the ends of the profiles were smeared due to the poor data coverage along the margins of the study area.

Discussion and Interpretation
The horizontal sections of the obtained seismic tomography model demonstrated a strong predominance of high-velocity anomalies in the southern part of the studied region beneath the Arabian plate. To the north of the Bitlis suture zone, we observed low-velocity anomalies that mark a boundary between the Arabian plate and the East Anatolian plate. Horizontal sections at a depth of 10 km (Figure 3) showed a clear correlation between the low-velocity anomalies and known tectonic features, such as the NAF, the EAF, and the MAFZ, which are highlighted in the Figures with black dotted lines. Deeper horizontal sections at depths of 30 and 60 km show that the entire study area to the north of the Bitlis suture zone is associated with low-velocity anomalies.
The four vertical profiles shown in Figure 4 are oriented north-south and are nearly perpendicular to the Bitlis suture zone. The first section crosses the Anatolian block, the second section crosses the NAF and EAF, the third section passes through the KTJ, and the fourth section passes through the EAAC and crosses the volcanic province around Lake Van.
In all vertical sections, the Bitlis suture zone was associated with a low-velocity anomaly at the surface, which is most likely due to the accretion of sedimentary and upper crustal rocks at the boundary between the Arabian and East Anatolian plates. Moreover, the vertical sections demonstrate some clear structural differences between the high-velocity Arabian plate to the south and the low-velocity East Anatolian plate to the north of the Bitlis border. Compared with the crust of eastern Anatolia, the crust of the Arabian plate was denser and more consolidated; therefore, high P-and S-wave velocities prevail.
In addition, according to the results shown in the vertical sections, a shallow highvelocity anomaly with a thickness of ~10-12 km was observed to the north of the Bitlis suture boundary, which was underlain by a thick low-velocity anomaly. Based on our model, it was difficult to evaluate the actual thickness of this deeper layer, as, according to the results of the synthetic test shown in Figure 7, the anomalies were smeared at great depth in the obtained models. This low-velocity anomaly has an inclined shape that dips northward from the Bitlis suture zone underneath the East Anatolian plate. This feature seems to be consistent with available receiver function results showing the gradual northward thickening of the crust [16,17]. According to [18], in the Bitlis suture zone, the Moho depth is 42 km, while in the northern part, for example, under the NAF, the Moho reaches 50 km in depth [19].
Beneath eastern Anatolia, the crustal structure appears to be atypical. In the upper part of the velocity sections, the velocities appear to be higher than those of regular continental crustal structures, whereas, in the lower crust, we observe prominent low-velocity anomalies. The higher-velocity structures in the upper crust might be associated with the widespread distribution of Neogene-Quaternary volcanism. Strongly consolidated intrusions and massive effusive flows have an integral effect of increasing the velocity in the upper crust.
The low velocities of seismic waves in the lower mafic crust can be explained by active shortening, which led to the entrainment of highly weakened material in the lower crust. This material was immersed to great depth, where it underwent phase changes and melting. An important additional factor affecting the lower crust is the absence of strong thinning of the mantle lithosphere beneath eastern Anatolia (see, for example, [20,21]). As a result, the lithospheric window is filled with overheated asthenospheric material that reaches the bottom of the crust. In this case, the low velocities in the lower crust in our tomography models may represent anomalous heating of the crust and possible melting due to elevated heat flow from the asthenosphere. Our findings are strongly consistent with the modeling results by [51], who presented a rheological model of eastern Anatolia and showed that the crust consists of a strong brittle upper crust (down to depths of 10-15 km) and a weak ductile lower crust (depths from 10-15 km to 40 km).
In the southern part of vertical profiles 1, 2, and 3 ( Figure 4), there was a high-velocity anomaly related to the Arabian plate, which has north immersion trending underneath the observed low-velocity anomaly. This positive anomaly was divided into two parts: one has an isometric shape and was located directly below the Bitlis suture zone, and another larger anomaly was located to the south. We propose that the former anomaly represents a remnant of the Neo-Tethys oceanic lithosphere remaining under the Bitlis suture zone, the main part of which broke off and plunged into the mantle during subduction. The larger part belongs to the crust of the Arabian plate. A schematic structure of the crust beneath eastern Anatolia according to the results obtained in this work is shown in Figure 8. In the southern part of the vertical sections, we can see that the crust of the Arabian plate is associated with a high-velocity anomaly reaching depths of more than 60 km, although the Moho boundary reaches only 40 km, according to [19]. This structure can be explained by the ongoing northward dipping of the Arabian plate. Since the crust of the Arabian plate does not differ structurally (thick sedimentary layer, many faults, etc.), the velocity anomaly throughout all of the investigated volume (up to a depth of 60 km) looks uniformly positive. In the northern part of vertical profile 1, a high-velocity anomaly stands out, which most likely belongs to a denser crust under the Pontides.

Conclusions
Using the LOTOS algorithm, we obtained a 3D model of heterogeneous crustal Pand S-wave velocities beneath eastern Anatolia down to a depth of ~60 km. We used a large amount of data provided by the T. C. Ministry of Interior Disaster and Emergency Management Directorate Earthquake Department Directorate of Turkey. The obtained tomographic models were verified using different tests showing high reliability and fair resolution.
The tomographic models demonstrate a clear separation of the study region into two tectonic units. To the south of the Bitlis suture zone, the crust of the consolidated Arabian plate is marked as a high-velocity anomaly. The low-velocity patterns to the north of the Bitlis suture are mostly associated with eastern Anatolia. According to the presented vertical sections, beneath the East Anatolian plate, we observed a thin (~10-12 km) high-velocity anomaly in the upper crust, which might be associated with massive Neogene-Quaternary magmatic activity in this area. In the lower crust, we observed an inclined lowvelocity anomaly having the smallest thickness beneath the Bitlis zone and the largest thickness beneath the northern part of the section (>50 km). This feature seems to be consistent with the general northward thickening trend of the crust, as revealed from existing receiver function studies.
It is most likely that the upper crustal high-velocity anomaly corresponds to the strong felsic layer of the crust, while the low-velocity anomaly is associated with the weak lower mafic crust. Thus, the presented results confirm the hypothesis of the existence of a lithospheric window beneath eastern Anatolia, through which hot material of the asthenosphere rises, and these findings help clarify the deep structure under eastern Anatolia.
Author Contributions: G.P. analyzed the data; I.M. and I.K. obtained seismic tomography models and contributed to preparing figures; All authors participated in discussions of the results and writing the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding:
The reported study was funded by Russian Foundation for Basic Research, project number 19-35-60002.

Data Availability Statement:
The full directory of LOTOS code with data corresponding to this study is available at https://doi.org/10.5281/zenodo.3944178.

Conflicts of Interest:
The authors declare no conflict of interest.