Three-Dimensional Seismic Velocity Structure of the Aegean Region of Turkey from Local Earthquake Tomography

SE111 ABSTRACT This study brings new insights to elucidate the three dimensional (3-D) seismic velocity structure of the Aegean region of Turkey by using Local Earthquake Tomography (LET). The study area has remarkable potential for geothermal resources. To provide the subsurface geological structure of seismogenic layers and geothermal areas, we develop new high-resolution depth-cross sections through Buyuk Menderes, Kucuk Menderes and Gediz grabens. Travel times of almost 265.000 readings (14.963 P-phase and 11.969 S-phase picks) from 2.085 well-located events recorded between 2007 and 2016 by a permanent seismic network of 75 broad-band seismometers were used to precisely interpret the 3-D Vp (lithological) and Vp/Vs (petrological) models of the study area. We detected four main layers down to 35-40 km depths with different Vp velocities ranging from 3.5 to 8.5 km/s. Mid-crustal interface (Conrad discontinuity) is discovered at 15 km depth in nearly all depth-cross-sections. Our results suggest an average of 25 km Moho depth in the Aegean region of Turkey. The depths range from around 18 km beneath SE of Aydin to 29 km beneath Aliaga, and approximate values of 19, 25 and 31 km beneath the Doganbey, Kutahya and Selendi-Demirci regions, respectively. The geothermal areas of the studied region are imaged by Vp and Vp/Vs anomalies. We propose the Aliaga, Denizli, Doganbey and Kosk areas as low Vp, low Vp/Vs anomalies which are an indicator of steam, CO2 or a mixture of both. Low Vp, high Vp/Vs models, suggesting geothermal fluids, are clearly visible near the Buharkent, Gumuskol, Guzelhisar, Kosk, Kuyucak, Saraykoy and Suzbeyli regions. We also report that the Bademler, Candarli, Kalekoy, Karahalli, Merdivenli, Ortakoy, Saruhanli, Yelki and Yuntdagikoseler regions might be good candidates for new potential geothermal resources.


Introduction
Turkey is one of the most seismically active areas in the world.Its tectonic and geological evolution has been managed by the repeated opening and closing of the Pa-leozoic and Mesozoic oceans.It is oriented within the Mediterranean Belt, whose complex deformation outcomes from the continental collision between the Eurasian and African plates.The border of these plates forms seismic belts marked by young active volcanos and faults, the next letting the transfer of water as well as heat.The distribution of hot springs in Turkey approximately parallels the young volcanism, fault systems and hydrothermally altered regions [Sanliyuksel and Baba 2011].
The seismic tomography technique is effective seismological tool to image the substructure of earth crust [Soldati et al. 2015, Cinar and Alkan 2016, Ibanez et al. 2016, Cinar and Alkan 2017, Parks et al. 2017].Moreover; geological/tectonic structures, geothermal resources, petrological properties of rocks and volcanoes can be identify using Local Earthquake Tomography (LET) method [Yang and Shen 2005, De Matteis et al. 2008, Kaypak 2008, Berger et al. 2011, Hofstetter et al. 2012, Yolsal-Cevikbilen et al. 2012, Del Pezzo et al. 2013, Muksin et al. 2013, Abril et al. 2016, Ozer and Polat 2017a, b, c].Studies [Kaypak and Gokkaya 2012, Komut et al. 2012, Erkan 2015, Bilim et al. 2016] on the crustal structure and geothermal potential of western Anatolia emphasise that different perspectives of the studies are still needed in the study area.Therefore; we used local earthquake tomography method to understand the crustal structure and geothermal potential of the Aegean region of Turkey.
Many studies performed to clarify the geological, geothermal and tectonic features of the region [Bozkurt et al. 1993, Bozkurt 2001, Gessner et al. 2001a, Gemici and Tarcan 2002, Okay 2008, Ersoy et al. 2010, Gessner et al. 2013, Ersoy et al. 2014, Ozer 2017, Ozer et al. 2017].Furthermore, some detailed geophysical studies were also applied in the region to shed light on geothermal and graben structures [Ilkisik 1995, Ates et al. 1999, Dolmaz et al. 2005, Drahor and Berge 2006, Duzgit et al. 2006, Sari and Salk 2006, Bilim 2007, Polat et al. 2008, Aktug et al. 2009, Isik and Senel 2009, Kaya 2010, Tufekci et al. 2010, Ekinci et al. 2013, Altinoglu et al. 2014, Gok and Polat 2014, Karakus 2015, Bilim et al. 2016].Salah et al. [2007] obtained 3-D tomographic images of the crust under the south-western part of Turkey by inverting a number of arrival time data of P and S waves.Mutlu and Karabulut [2011] determined crustal thicknesses as 28 ±2 and 33 ±2 km for the western Anatolia and the Aegean Sea, respectively.Komut et al. [2012] asserted that the thin Anatolian-Aegean lithosphere is being buoyed upwards by underlying mantle flow.The mantle flow might be related to active lithosphere delimitation; an operation that could also account for the continuing crustal extension.Salaun et al. [2012] focused on a low velocity zone (80-200 km depth) that reflects a slow and warm asthenosphere underlying a thin lithosphere.Kaypak and Gokkaya [2012] showed 3-D Vp and Vp/Vs seismic velocity structure for the upper 20 km of the crust beneath Denizli basin.Karabulut et al. [2013] defined that the crust-mantle Figure 1.Geological map and main tectonic units of the Aegean region of Turkey [compiled from Emre et al. 2013, Seghedi et al. 2015, Erkan 2015, Kaypak and Gokkaya 2012].Geothermal sources are given at Table 1 in detail.High-resolution topography and bathymetry data sets have been compiled from Farr [2007] and NGDC [2006].Black and white solid lines indicate the main fault zones (EAFZ: East Anatolian Fault Zone, KJ: Karliova Junction, NAFZ: North Anatolian Fault Zone, IFZ: Izmir Fault Zone).
boundary is flat beneath Menderes massif at ~25 km depths.Vanacore et al. [2013] claimed that high Vp/Vs ratio (>1.85) in western Anatolia might be indicative of partial melt in the lower crust affiliated with regional extension.Delph et al. [2015] studied S-wave characteristics and calculated average crustal shear wave velocity as 3.3-3.5 km/s in the region.
Turkey also exhibits high potential for geothermal resources in Europe [Baba A. 2015, Erkan K. 2015, Karakus H. 2015].Western Turkey, especially, holds enormous potential for improving the geothermal capacity of the country.The present geothermal regime, which reveals high temperatures [120-240 °C), is mainly concentrated in the Menderes massif.This huge potential is related to faults and heat flows from horstgraben systems [Kose 2007, Baba and Sozbilir 2012, Parlaktuna et al. 2013, Baba et al. 2014, Erkan 2015].
We aim to identify the 3-D seismic velocity structure of the Aegean region of Turkey down to 40 km by using a large earthquake (~25.000)database using for the first time to determine the seismic velocity structure.This paper also brings interpretations for possible new geothermal areas not drilled yet in the region by conducting vertical several profiles using a dense seismic array crossing BMG, KMG and GG.The LET technique was used for the first time to identify crustal structure variations mostly belong to lithological and petrological features of the geological units.

Tectonic, Geological and Geothermal Settings
The Aegean region of Turkey, located on the Alpine-Himalayan Belt, exhibits an extensional regime in a north-south direction, and this motion creates different sizes of east-west trending grabens such as Buyuk Menderes Graben (BMG), Kucuk Menderes Graben (KMG) and Gediz Graben (GG).It has been recognised that the tectonics of the study area are controlled by these grabens.Furthermore, a number of secondary strike-slip faults extending generally in NE-SW and NW-SE directions exhibit active seismic activity and manage the evolution of a few Miocene to recent sedimentary basins [Uzel and Sozbilir 2008, Sozbilir et al. 2011, Uzel et al. 2013].
The normal faults that have low angles (detachment faults) are related to a metamorphic core complex, the Menderes Massif, and the E-W and NE-SW basins symbolise the tectonics of the Aegean region of Turkey.
Moreover, some studies identify the existence of some NE-SW trending strike-slip faults controlling the NEtrending Miocene deposition.One of the important areas is the Menderes Extensional Metamorphic Complex (MEMC).The MEMC is bounded by NE-SW strike-slip faults, suggesting that it is not comprised of domal uplift through crustal-scale low-angle extensional detachment faults [Cemen et al. 2014, Ersoy et al. 2014].
The major E-W grabens and basins include the Kucuk and Buyuk Menderes grabens and the Gediz graben.The Quaternary Gediz Graben is a high-angle normal fault evolved at the northern end of the Alasehir detachment basin which contains deformed granodiorite coming from deeper parts.This progressive uphill replace from the undeformed granadiorite to a brittlely deformed detachment surface recommend that extension regime in a ductile deformation at depth and the ductilely deformed granitoids were carried to shallower depths.The content of the Alasehir shear zone mostly consists of mylonitic rock such as granodiorite, mafic and felsic dykes.The central Menderes massif is separated in the north and the south by the Gediz Graben and Buyuk Menderes Graben detachments, respectively, and is broken into pieces in the middle by the E-W trending Kucuk Menderes Graben.The north side of the Buyuk Menderes Graben is bound by a steeply dipping major fault that separates the Neogene sediments of the graben from the metamorphic sequences of the Menderes massif.These E-W trending grabens are bounded by high-angle to moderately dipping normal faults, some of which are seismically active [Isik et al. 2003, Ciftci and Bozkurt 2009a, Ciftci and Bozkurt 2009b, Oner and Dilek 2011].
The Miocene sediments and the remnants of late Miocene erosion surfaces (Figure 1) are flat-lying and parallel to each other, spreading along NE-SW direction and sometimes related to basaltic volcanic rocks [Gessner et al. 2013, Seghedi et al. 2015].The Miocene volcanic rocks are observed in Foca, Aliaga, Cesme, Ayvalik and surrounding areas.In the Foca-Aliaga-Izmir area, the Early Miocene volcanic rocks of the Yuntdag volcanics overlie the alluvial-fan and fluvio-lacustrine sedimentary units of the lower sedimentary succession.The Early-Middle Miocene sedimentation was controlled by NE-SW and NW-SE-striking strike-slip faults [Karacık et al. 2007, Altunkaynak et al. 2010, Ozkaymak et al. 2013].The NE-SW-trending zone of deformation along the western margin of the Menderes Core Complex is known as the Izmir-Balikesir Transfer Zone.The Izmir-Balikesir Transfer Zone symbolises the NE extension into western Turkey of a NE-SW trending shear zone [Erkul et al. 2005    tems such as Buyuk Menderes Graben, Kucuk Menderes Graben, Gediz Graben, Orhanli Fault Zone, Gulbahce Fault Zone and Denizli fault zone.Several geothermal studies have been conducted to understand the geothermal capacity of the Aegean region using different methods.Important geothermal areas (Table 1) of the study region have been compiled from different studies [Basel et al. 2009, Serpen et al. 2009, Erkan 2015, Kiyak et al. 2015].

Data, Algorithm and Checkerboard test
We supply a database containing phase readings available from AFAD (Earthquake Department of the Disaster and Emergency Management Authority) for the earthquakes located in the study area between 2007 and 2016.More than 25.000 earthquakes have been used for tomographic calculation (Figure 2).We mainly detect that earthquake depth distribution is concentrated in the first 25 km of the earth's crust.However, in rare cases some earthquakes are located down to 50 km.We perform some selection criteria to improve data quality.The first is total P and S picks for each event, which must be greater than or equal to 15.To guarantee enough azimuthal ray coverage for the whole data set, some additional selection criteria; the gap in the azimuthal coverage of each earthquake had to be lower or equal 200º; the maximum adopted residual time is 1.0 s.The LOTOS algorithm [Koulakov 2009], which based on the observed P-and S-wave travel times to invert simultaneously for a 3-D image of the P-and Swave velocities, hypocentres, and station corrections.LOTOS suggests two choices, inverting either for the Pand S-wave velocities or for the P-wave velocities and Vp/ Vs velocity ratio.This study outcomes were calculated with the second choice [Dzierma et al. 2012].The low values of residuals after the sixth iteration approve the high quality of the phase picking in this study.Additionally, an earthquake is refused if the horizontal distance to the nearest seismic station is more than predefined distance (e.g. 100 km).P-and S-phases were selected in an admissible correctness for inversion procedure (Table 2).The final data set comprises 14.963 Pand 11.969 S-phase arrivals from 2.085 local events recorded by 75 stations (Figure 3).The digital seismic stations consist of three-component broadband recorders.
Velocity models (1-D and 3-D) can be determined by many tomography algorithms using local earthquakes [e.g., Thurber 1993, Zhao 2009, Li et al. 2016].In the frame of this study, we have used LOTOS code [Koulakov 2009].It uses a grid search method to determine the source coordinates and origin times.The main aim is to obtain a maximum goal function (GF) in the 3-D space [Koulakov and Sobolev 2006].The GF shows the probability of a source location in the existing catalogue.The computation starts with initial earthquake locations in a 1-D velocity model at a preliminary step.This option makes initial location stage comparatively balanced and fast.Earthquake locations vary iteratively with the optimisation stage in the one-dimensional (1-D) model (Table 3) based on matrix inversion for Vp and Vs velocities, earthquake coordinates and origin times.After conducting an initial location, a major tomographic process optimises the 1-D model in the first iteration and then relocates the seismic 3-D velocity model by iterating the following models.For this case, several grid orientations (e.g.0°, 22°, 45°, 67°) increase the quality of the selected grid configuration for the tomographic calculations.The code uses the LSQR method [Paige and Saunders 1982] to conduct matrix inversion for P-, S-and source parameters (dx, dy, dz, dt).Additionally, two matrix blocks control the amplitude and smoothness of solution, which directly affect the results.These parameters and other free parameters are determined using synthetic modelling tests.Also; Model regularisation can be controlled in two ways: by choosing the weighting of the smoothness and by the number of LSQR steps in the inversion [Tryggvason et. al. 2002, Jeddi et al. 2016].
A basic image of the probable resolution of the data set is able to acquire by drawing the ray paths applied in the inversion.The ray path coverage is sufficient to understand seismic velocity structure under the central part of our study area.The local seismic rays is not concentrated on the eastern, western, northern and southern borders of study area.These indicate that sufficient numbers of rays covered in the central, densely  gridded, part of the study area (Figure 4).In addition, our checkerboard resolution tests show poor resolution on the eastern, western, northern and southern borders of study area.Hence, the vertical depth cross-sections were designed different length according to our tomographic model.We applied the checkerboard test, which shows theoretical capacity of inversion algorithm to rebuild amplitude recovery [Koulakov 2009] and calculate the impact of noise on resolution using synthetic travel times.In checkerboard tests, we defined periodic negative and positive velocity anomalies of their different dimensions and performed these tests to resolve the noise effect and optimum values of inversion.For the same source-receiver pairs, we computed synthetic travel times as in the case of the observed data to understand amplitude recovery.We then followed a similar process for real data analysis, introducing the absolute source location [Ozer and Polat 2017b].We designed four checkerboard tests using different block sizes, keeping the earthquake locations and observation number of earthquakes, and joining ±10 velocity anomalies considering 1-D seismic velocity model.We added random noise (normal Gaussian distribution) proportional to P-and S-phase reading uncertainties to reflect real conditions.
The checkerboard test synthetic travel times are calculated for the identical earthquake location and seismic station pairs.The real source location corresponding to the solution is acquired after six iterations of real data solution.Noise is added to the synthetic data by disturbing real readings.The random noise that have a realistic distribution histogram constructed traditionally used normal Gaussian distribution.The average amplitude of noise was determined at 0.15 s and 0.3 s for P and S synthetic travel time data, respectively.Once synthetic travel time is calculated, we discard the velocity model and earthquake location information.Then the code reconstructs the final earthquake location using the same process steps and inversion parameters as in the form of real data analysis (Figure 5).
To calculate the impact of the noise and inversion parameters on resolution in the study area, we conducted checkerboard tests using different sizes of anomaly shape.Periodical negative and positive anomalies (±10) of different sizes were described for four different test models.The sizes of anomalies with 5 km empty spacing in Model 1, 2, 3 and 4 were 50×50×40, 30×30×40, 20×20×40 and 10×10×40 km; respectively.0.15 s and 0.30 s random noises were added for P-and Stravel times, respectively.Note that the resolved areas continue down to 30 km depth.Resolution decreases OZER ET AL. when the box sizes become smaller towards the deeper parts.According to the checkerboard test results, velocity anomalies smaller than 20 km are not resolved (Figure 5).
The another synthetic model with realistic patterns of horizontal section for P wave velocity anomalies is shown in Figure 6.The shape of samples kept mostly the same, however the amplitudes were 5% higher in some parts.The seismic anomalies are described in prisms that may have a complex shape when seen in map view.The shapes and amplitudes of the synthetic patterns are reconstructed as we expected and the model is quite similar to the actual initial anomalies.Consequently, 3D VELOCITY STRUCTURE OF AEGEAN REGION FROM LET synthetic tests show that the tomographic models are comparably well resolved and our tomographic model has enough resolution to make interpretation using absolute seismic velocities.

Tomography Results
We have taken six Vp and Vp/Vs cross-sections (Figure 3) through vertical profiles which are almost parallel and perpendicular to each other in revealing crustal OZER ET AL.Farr [2007].The white parts indicate unsolved area because of lack of S-wave data.Smoothing factor are used to stabilize the tomographic inversion result.Optimum smoothing values can be determined trying different synthetic modelling.The inversion of this sparse matrix is calculated using the leastsquares QR (LSQR) method [Paige and Saunders 1982].Koulakov [2009] claimed that the most effective method to estimate optimal values of free inversion parameters is by conducting synthetic modelling that reproduces the real situation.This also lets qualitative guesses of amplitudes of seismic anomalies in the real Earth.It was decided to use smoothing factor as 0.4 after several synthetic tests.
velocity structures of the Aegean region of Turkey.There of them (profiles 1, 2, 3) were selected to bring new interpretations for potential geothermal resources along Buyuk Menderes, Kucuk Menderes and Gediz Grabens.And the last three slices cross horst and graben systems oriented N-S direction to reveal seismic velocity structure of the study area.The accuracy of the tomography produced from the 3-D seismic velocity model is shown with travel time residuals (RMS).They decreased by 44.44% from 0.27 s to 0.15 s for P-wave and 46.67% from 0.45 s to 0.24 s for the primary 1-D model (Table 2).
The absolute Vp values and Vp/Vs ratios usually refer to lithological and petrological features, respectively.To discover the seismic velocity structure of Aegean region of Turkey using local earthquake tomography method, we have produced vertical depth cross-section of the Vp and Vp/Vs models previously unavailable for the study area (Figure 7 and 8).The Vp ranges from 3.5 to 8.5 km/s, and the Vp/Vs change between 1.5 and 2. Four main seismic layers have been precisely detected down to 30 km.These are upper (~0-8 km), middle (~8-15 km), lower (~15-30 km) crustal and brittle layers (>30 km).The tomographic results indicate that the main Moho discontinuity may lie around a depth of ~29 km.
It is evident that some geothermal areas located near the Aliaga (profile 1), Sindirgi (profile 2), Saraykoy-Karahayit and Demirci-Simav (profile 3) and Denizli (profile 4) regions reduce seismic velocities in the first seismic layer through the active fault surface in the places of deformation.Similarly, we have observed low Vp anomalies along profiles 2, 3, 4, 6, at detachment fault (e.g., Figure 7), and at the edges of the basin formed by the IFZ (profile 1, Figure 7).

High Vp anomalies (>6.5 km/s)
High P-wave velocities (>6.5 km/s) detected in the upper crustal layer might be associated with geological features observed in shallow depths [Masson et al. 2000].While the high P-velocities are clearly associated with the Menderes massif around profile 2A, they represent Miocene volcanics towards profile 2B.Some high Vp values are also computed close to the surface along profiles 2, 3, 5 and 6 in Figure 7. Another clear high Vp zone surrounded by low Vp anomalies beneath Soke-Germencik is also visible between 3-15 km depths in profile 4A.It could represent a cycladic zone consisting of marbles, metavolcanic rocks by indicating Eocene metamorphism.All these type of velocities can be explained as magmatic materials apart from the deeper depths to the near surface [Baris et al. 2005, Uzel et al. 2015].Kaypak and Gokkaya [2012] have also reported a similar observation for Denizli in the Aegean region of Turkey, and interpreted this feature with paleotectonic blocks.Some tomograms are also traceable towards the deep edge of the basin near the same area (Denizli) as seen from profiles 3 and 4 (Figure 7).The Moho discontinuity, where high Vp values greater than or equal to 6.8 km/s, is detected at ~29 km depths (dashed black line in Figure 7) as an undulated shape, and this feature is coherent with previous studies [e.g., Bilim 2007, Kaya 2010, Vanacore et al. 2013, Karabulut et al. 2013, Bilim et al. 2016].

Low Vp, low Vp/Vs anomalies (Vp <4.5 km/s, Vp/Vs <1.65)
Most studies explain that this type of anomaly may indicate gas dominated rock, and this gas might be made of steam, CO 2 or a mixture of both as a result of magmatic events (Hauksson 2000, Kaypak andGokkaya 2012].The CO 2 arising from the decomposition of carbonate rocks is an outcome of the thermo-metamorphic or hydrolysis reaction in the marbles and limestone [Hauksson 2000, Kaypak andGokkaya 2012].It therefore seems reasonable to guess that this type of anomaly could also be an indicator of geothermal areas.In the present study, low Vp (<4.5 km/s) and low Vp/Vs (1.65) models are well correlated with the geological structures and petrological characteristics of the region.We identified low P-velocities and low Vp/Vs values at almost all cross-sections in Figures 7 and 8.In profile 1, we detect very low Vp and Vp/Vs values near Doganbey geothermal region (beginning of the profile) down to 6 km, in Menemen geothermal area (middle of the profile) down to 7 km, and at the S of Aliaga geothermal region down to 5 km.In profile 2, the Kosk geothermal area located at the bound of BMG is notably harmonious with these anomalies.The SW of Bigadic exhibits normal to low velocities, which indicates gas content.The beginning of profile 3 illustrates the release of CO 2 apart from magmatic material and low Vp, low Vp/Vs are especially 3D VELOCITY STRUCTURE OF AEGEAN REGION FROM LET traceable near the Denizli region down to 8-10 km, as also reported by Kaypak and Gokkaya [2012].The low Vp/Vs sections located near 29oE, 38oN (see profiles 3A, 4B, 5B, 6B in Figure 8) clearly reveal the enormous geothermal energy capacity of the Denizli region.
The Doganbey and Seferihisar areas demonstrate low P-velocities and low Vp/Vs rations at the beginning of profile 1 (Figure 8).We also described low P-velocity and normal (~1.75)-to-low Vp/Vs values near the periphery of Alasehir and Kula (profile 6, Figures 7 and 8) overlaying mostly the GG.(Vp <4.5 km/s,Vp/Vs >1.85)These types of anomaly are an indicator of geothermal liquid.The source of high Vp/Vs (>1.85) values is related to reduced S-wave velocities and evaluated as saturated, extremely fractured and huge liquid pressure.Accordingly, these anomalies might be associated with geothermal liquid [Hauksson 2000].In profile 1 of Figures 7 and 8, large scale low Vp and high Vp/Vs values can be followed along Suzbeyli (SW of the Menemen), Guzelhisar (SE part of the Aliaga) and Kizilcukur (Dikili) regions reported as rich geothermal resources [Erkan 2015].In addition, we explore some new resources at Çandarlı and Merdivenli regions located on the extinct Karadag volcano, Bademler and Yelki areas.

Low Vp, high Vp/Vs anomalies
The pipe type anomalies going down to several kilometres might be an indicator of new geothermal resources.However, additional detailed investigations need to be performed in the periphery of those areas.The high Vp/Vs ratio verified the opinion of a saturated and extremely fractured zone through the margin of BMG (N of Kosk area) in profile 2, Figure 8.Other high Vp/Vs values are also detectable down to 12 km between Akhisar and Sindirgi (profile 2), indicating ophiolitic rocks, water-saturated cracks and geothermal liquid.In profile 3, Gumuskol geothermal area, located on Miocene volcanic, display high Vp/Vs values by revealing the presence of liquid or liquid pressure down to deeper parts.The Hamamdere region, located 20 km from Gumuskol, may show volcanic units down to several kilometres.The NE of Simav area, formed by Lycian nappes, has hot geothermal liquid (100 °C) and shows low Vp and high Vp/Vs values.The result demonstrates the huge geothermal capacity of Simav region as clearly seen from the end of profile 3 (Figures 7 and 8).Profiles 4, 5 and 6 are specifically designed to investigate the geothermal features of Menderes Massif, which consists of BMG, KMG and GG.The BMG that is formed by normal fault systems contains many geothermal resources located between the Aydin and Denizli areas.We detect that the top layer of sedimentary units, which is represented by the lowest P-velocities (<4.5 km/s), lies down to 3 km, except W of the Guzelhisar-2 and Kosk areas, where thicknesses are present much deeper (profile 4, Figure 7).Some studies report the thickness of the top sedimentary layer as being 3.9 km at its maximum [Isik andSenel 2009, Baba andSozbilir 2012].In profile 4 of Vp and Vp/Vs, some geothermal areas (e.g.Germencik, Kuyucak, Kosk, Buharkent, Saraykoy) are characterised by low Vp and, high Vp/Vs values through BMG.The KMG has graben geometry 80 km long and 10 km wide.At the basement of KMG are formed pre-Miocene age metamorphic rocks and also andesitic and basaltic rocks.Several hot spots occur near Bayindir and Odemis in KMG [Gemici andTarcan 2002, Baba andSozbilir 2012].We report low P-wave velocities and high Vp/Vs ratios in most regions (e.g.Alasehir, Karahayit, Salihli, Sarigol) of the KMG except Doganbey, Seferihisar and W of Sapcilar (profile 5, Figures 7 and 8).As another important basin, the GG is 140 km long and 40 km wide.It comprises quaternary alkaline volcanism in the northern part down to 3 km by the exhibition of high heatflow values [Dolmaz et al. 2005].Normal faults also control the GG, which consists of gneiss, marble, schist and the Menderes core complex.The highest temperatures near Alasehir (215 °C) and Salihli (287 °C) grabens are observed in alkaline thermal waters.The Karahayit area contains different types of rocks (e.g.marble, schist, quartzite) located between BMG and GG near the Denizli province.This geothermal water is probably dominated by a mixture of cold waters, mineral dissolution and saturation reactions.In addition to graben systems, we also report that horsts in the GG are also good candidates for geothermal resources, as they exhibit low Vp, high Vp/Vs near Kursunlu, Caferbey and Urganli (85 °C) areas (profile 6, Figure 8).Surface temperatures are reported in hot springs from 50o to 150oC.We report Candarli, Yuntdagikoseler, N of Kalekoy, Saruhanli and Sarigol as new possible geothermal resources as they exhibit low Vp and high Vp/Vs values.However, we also recommend that these areas should be investigated in detail through further geophysical studies.

High Vp, high Vp/Vs anomalies (Vp >6.5 km/s, Vp/Vs >1.85)
These types of anomaly are observed especially in deeper layers in certain cross-sections.Hauksson [2000] asserted that the logical explanation for high Vp and high Vp/Vs at deeper parts is related to the existence of mafic rock such as gabbroic intrusion rather than the presence of liquid.W of Cine (profile 2), in the N and S part of Simav (profile 3), N of Soke (profile 4), and SW of Sarigol (profile 5) may indicate magmatic rocks at deeper parts of the above cross-sections (Figures 7 and 8).

Moho depth
The Moho depth of the study area is reported as being between 25 and 35 km in various studies [e.g.Kaya 2010, Tezel et al. 2010, Mutlu and Karabulut 2011, Kaypak and Gokkaya 2012, Karabulut et al. 2013, Vanacore et al. 2013].Some studies observe an increase in Pwave velocities between 6.5 and 7.2 km/s at Moho discontinuity [Hauksson 2000, Cambaz and Karabulut 2010, Salaun et al. 2012, Vanacore et al. 2013, Yolsal-Cevikbilen et al. 2014, Delph et al. 2015].We computed Pwave velocity as 6.8 km/s at an average of 25 km Moho depth from almost all depth cross-sections (Figure 7).Our tomographic model reveals that the Moho depth ranges from 20 to 30 km.The lowest values of Moho depth are generally observed beneath geothermal areas (Figure 9).Shallow Moho depth distributions (15-20 km) are observed beneath the western part of Aydin, SW of Izmir and SW of Kutahya (S1, S2, S3, respectively, in Figure 9).Several geothermal reservoirs exhibiting high temperature are located around these areas.We believe that these blocks are related to magmatism and the magma material is transported from the deeper depths to the surface areas [Gessner et al. 2013, Ersoy et al. 2014, Seghedi et al. 2015, Uzel et al. 2015].Deep Moho anomalies (30-35 km) are located near Aliaga, and W and S of Usak (D1, D2, D3, respectively, in Figure 9).Aliaga region exhibits many hot spots and this finding can be explained by the presence of Miocene age volcanic units [Gessner et al. 2013, Ersoy et al. 2014, Seghedi et al. 2015].D2 and D3 anomalies are considered as thick crustal structures in the higher topographic elevation.It is acceptable that we have not yet observed any magma contact or other magmatic activity at D2 and D3 anomalies in depth cross-sections.Conversely, the brittle zone (Vp<6.8km/s) is comparably thicker than other regions, and the ductile zone (Vp>6.8km/s) starts below depths of 30-35 km (Figure 9).
It is evident that some geophysical studies such as those of microgravity and magnetic conduction through grabens and mountainous areas, are in line 3D VELOCITY STRUCTURE OF AEGEAN REGION FROM LET with the present velocity values acquired in the frame of this study.Ates et al. [1999] conducted gravity and magnetic tests and obtained gravity values decreasing gradually from West to East.These results clearly fit with our results.Crust-mantle interface values prove shallower depths at the S1-S2 area and deeper depths in the D2-D3 region.The Bouguer gravity anomaly map shows a positive anomaly at S1 and S2 regions, while negative anomalies are observed at the D1-D2 and D3 areas of this study.Bilim et al. [2016] report the positive residual aeromagnetic anomaly, which exactly matches with the S2 and S3 areas in our study.Their Curie point depths also have good harmony with our findings, especially for S1 and S2 Moho depths [e.g.Aydin and SW of Izmir].

Conclusions
This study is a first attempt to perform the LET algorithm to develop the 3-D seismic velocity structure of the Aegean region of Turkey down to depths of ~35-40 km.This research also brings important features for present geothermal resources and some new potential/target reservoirs which have not yet been drilled.The results are obtained using more than 25.000 earthquakes recorded by 75 seismic stations operated by the Disaster and Emergency Management Authority (AFAD, Ankara).We have generated six depth cross-sections (Vp and Vp/Vs) along E-W and N-S directions in the study area to precisely illustrate velocity characteristics beneath grabens.We have identified four main crustal layers based on P-wave seismic velocity anomalies, which range from 3.5 to 8.5 km/s.The checkerboard tests indicate good resolution down to 40 km.Features deeper than 40 km are not well-calculated due to the insufficient or shallow depths of earthquakes.The lowest Vp velocity (<3.5 km/s) in shallow areas (<3-5 km) is identified as quaternary alluvium.Furthermore, we used the Vp/Vs ratio to interpret petrological features such as geothermal gas, and fluid contents beneath present or potential geothermal areas.While the gas-filled or gas-saturated zones could usually be represented by low Vp and low Vp/Vs anomalies, the low Vp and high Vp/Vs values may indicate high fluid contents through the fault zone.The high Vp/Vs ratio (>1.85) may indicate fluid contact at shallower depths.We report some intrusive and dense magmatic rocks with high Vp values (Vp > 6.5 km/s) at shallower depths (e.g.Bademler, Candarli, Karahalli, Kalekoy, Merdivenli, Ortakoy, Saruhanli, Yelki and Yuntdagikoseler).These locations may indicate new resources for geothermal fluid reservoirs.In addition to these, Aliaga, Denizli and its surroundings, and the Do-ganbey and Kosk areas could be potential sources for CO 2 , another gas or a mixture of them.This study reports seismogenic depth to be between 0 and 15 km and average Moho depth is range from 25 to 35 km in the study area.We observed that it changes between 20 km and 30 km depths beneath horst and graben structures.

Figure 2 .
Figure 2. Locations of all 25.000 earthquakes (hand-picked, M≥2.0) recorded between 2007 and 2016 for the study area.Stations, geothermal sources and event distributions are represented by filled blue triangles, green stars and light gray circles, respectively.Red lines show fault traces [digitized after Emre et al. 2013].Earthquakes generally concentrate in the first 35 km depths of the earth crust as clearly visible from vertical/horizontal depth cross-sections and histogram plot which denotes number of earthquakes versus depth.Black lines indicate the coast of study area.

Figure 3 .
Figure3.Locations of well located 2.085 earthquakes selected in the present study.Study area is well covered by recorded earthquakes (circle) and station (triangle) distributions.We selected 6 profiles (thick solid lines).3 of them (profiles 1, 2, 3) were chosen to bring new interpretations for potential geothermal resources along Buyuk Menderes, Kucuk Menderes and Gediz Grabens.And the last 3 slices cross horst and graben systems oriented N-S direction to reveal seismic velocity structure of the study area.High-resolution topography and bathymetry data sets have been compiled fromFarr [2007] and NGDC[2006].The average location uncertainty (RMS<0.30) is approximately ±2.5 km.The earthquakes with large location uncertainty (>5 km) were eliminated in this study.

Figure 4 .
Figure 4. Ray path coverage of selected earthquakes in horizontal planes for the P-wave inversion.Every path between an earthquake and a station was drawn as one straight line.The seismic stations and rays are shown by blue triangles and red lines respectively.

Figure 5 .
Figure 5. Periodic negative (blue color) and positive (red color) P-wave anomalies display ±10% velocity perturbations: 50, 30, 20 and 10 km box sizes for Model 1, 2, 3 and 4, respectively, with 5 km empty space.Input for each model was given at the top of each column.Synthetic (checkerboard) tests are applied to calculate the spatial resolution and to predict the optimal values of inversion parameters down to 40 km depth layers.Thick solid lines represent profiles.

Figure 6 .Figure 7 .
Figure 6.Synthetic P-wave model for testing the efficiency of our tomographic model and amplitude recovery.The synthetic model forms seven letter prisms described in one horizontal section.

Figure 8 .
Figure 8. Depth cross-sections of Vp/Vs models which denote petrological characteristics of the study area.Well resolved ratio values are changing between low (<1.65) and high (>1.85)values down to 30 km depth layers.Interpretations of low or high Vp/Vs ratios together with P-wave velocities of sub-crustal structures are discussed in detail in the text.Red stars and dark grey circle represent discovered geothermal areas and well located earthquakes, respectively.High-resolution topography data sets have been compiled from Farr [2007].The white parts indicate unsolved area because of lack of S-wave data.Smoothing factor are used to stabilize the tomographic inversion result.Optimum smoothing values can be determined trying different synthetic modelling.The inversion of this sparse matrix is calculated using the leastsquares QR (LSQR) method[Paige and Saunders 1982].Koulakov [2009] claimed that the most effective method to estimate optimal values of free inversion parameters is by conducting synthetic modelling that reproduces the real situation.This also lets qualitative guesses of amplitudes of seismic anomalies in the real Earth.It was decided to use smoothing factor as 0.4 after several synthetic tests.

Figure 9 .
Figure 9. Distribution of Moho depths for equal or greater than 6.8 km/s reference P-wave velocity.D1, D2 and D3 contours show increase in Moho depth layers while S1, S2 and S3 present decreased depths (or intrusive bodies) towards shallower parts.D3 is not well represented by a clear contour due to the weak data coverage.Result suggests an average of 25 km Moho depth and it ranges between 19 and 31 km depths beneath the study area.

Table 2 .
RMS (s) values of the P-and S-wave residuals after six iterations for Vp and Vs inversion procedures.
3D VELOCITY STRUCTURE OF AEGEAN REGION FROM LET

Table 3 .
The starting 1-D Vp velocity model used in the tomographic inversion[Ozer and Polat 2017].