Isostatic Moho undulations and estimated elastic thicknesses of the lithosphere in the central Anatolian plateau, Turkey

Accepted Manuscript Isostatic Moho undulations and estimated elastic thicknesses of the lithosphere in the Central Anatolian plateau, Turkey Bülent Oruç, Oya Pamukçu, Tuba Sayın PII: DOI: Reference:

S1367-9120(18)30454-1 https://doi.org/10.1016/j.jseaes.2018.11.001 JAES 3692

To appear in:

Journal of Asian Earth Sciences

Received Date: Revised Date: Accepted Date:

3 May 2018 31 October 2018 1 November 2018

Please cite this article as: Oruç, B., Pamukçu, O., Sayın, T., Isostatic Moho undulations and estimated elastic thicknesses of the lithosphere in the Central Anatolian plateau, Turkey, Journal of Asian Earth Sciences (2018), doi: https://doi.org/10.1016/j.jseaes.2018.11.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Isostatic Moho undulations and estimated elastic thicknesses of the lithosphere in the Central Anatolian plateau, Turkey Bülent Oruça*, Oya Pamukçub, Tuba Sayına

a

Kocaeli University, Engineering Faculty, Department of Geophysical Engineering, İzmit, Kocaeli, Turkey

b

Dokuz Eylül University, Engineering Faculty, Department of Geophysical Engineering, Buca, Izmir, Turkey

Abstract In this paper, the flexural isostatic deformation of the crustal structure and an average effective elastic thickness (EET) of the lithosphere in the central Anatolian plateau has been estimated. We have used the Vening-Meinesz model to calculate the isostatic Moho depths or crustal thicknesses due to topographic loads. The maximum Moho deepening is under the Alpine Orogenic belt throughout the southern part of Anatolian plate, whereas the Moho depths are relatively shallow in the northwestern and central parts of the region. An average EET of central Anatolia’s lithosphere is calculated by multispectral coherence between EGM08 Bouguer anomalies and ETOPO1 topographic data sets. Application of spectral coherence estimation provides an average EET value of 21.3 km and hence points out that the lithosphere is deflected. In addition, density contrasts of the horizontally and vertically extended in the uppermost crust were estimated using compact gravity inversion. We interpret the uppermost crustal blocks in terms of density contrasts may contribute to the total crustal load in addition to topographic loading.

Keyword: Isostatic Moho, regression analysis, effective elastic thickness, gravity inversion

1.

Introduction

Isostasy is defined as the state of equilibrium in response to mass loads which can be defined as surface and subsurface forces. The dynamic structures of the earth and load distribution are the processes which influence isostatic compensation. Crustal deformation is thus concerned with responses of the Earth's crust to tectonic processes. Several hypotheses were proposed to explain the isotactic compensation mechanism (Pratt 1854, Airy 1855; Vening Meinesz 1931, 1940). Airy isostasy model assumes no flexural support and proposes that topographic loading

is compensated by buoyancy forces resulting from lateral variations in crustal thicknesses. Vening Meinesz, or flexural isostasy model modifies the Airy isostasy theory by introducing a regional isostatic compensation model based on a thin plate flexure. A method to calculate the crustal thicknesses by considering a regional compensation scheme of topographic masses was developed by Vening Meinesz (1931). Parker (1972) has described a forward algorithm of potential field anomalies caused by a deeper interface with Fourier transforms. Oldenburg (1974) has rearranged Parker’s method and developed a gravimetric inversion technique in order to estimate for a deeper interface depth. Deflections at the Moho boundary are useful to reveal the crustal deformation in different tectonic regimes concerning the crustal stresses and seismic risk analysis (Pamukçu and Yurdakul, 2008; Pamukçu and Akçığ, 2011; Naliboff et al., 2012; Oruç et al., 2017) in addition to the geodynamic modeling (Beaumont et al., 2000; Grobys et al., 2008; Pfiffner et al., 2000; Valera et al., 2011).

The effective elastic thickness, EET, of the lithosphere is defined as a cumulative thickness of an equivalent elastic plate under the vertical loads (topographical masses) and is a proxy for lithospheric strength. Loads are supported by the lithospheric strength, characterized by EET. Forsyth (1985) developed a coherence method by evaluating the statistical relationships between the topography and Bouguer gravity spectra. The spectral coherence between both data has been widely used to estimate an average value of EET (Zuber, 1989; Betchel, 1990; Macario et al., 1995; Gomez-Ortiz et al., 2005; Oruç et al., 2017). The EET is a critical factor for understanding lithosphere dynamics in addition to the calculation of the isostatic Moho depths in a manner of flexural isostasy model. The land gravity and magnetic surveys have been carried out by many researchers (Ateş et al., 2005; Aydemir, 2009; Aydemir, 2011; Bilim et al., 2015) in mapping geological structures in the central Anatolian region. The three-dimensional gravity model of the Haymana, Tersakan and Tuzgolu basins has been estimated by Aydemir and Ateş (2006, 2008). Aydemir (2011) used the first-order vertical derivative of the gravity anomalies to analyze the hydrocarbon generation potential in the Haymana basin. In this work, we generate a flexural model of central Anatolian’s Moho that satisfies regional isostasy model. For this purpose, an average EET is estimated from Bouguer-topography coherence as a function of wavenumber in order to examine the rigidity of the elastic lithosphere. The relationships between Moho depths and surface topography with regression

analysis are investigated. A two-dimensional density contrast distribution of uppermost crust is also inferred in order to assess the crustal deformation.

2.

Brief outline of the tectonic history of the central Anatolian region

The neotectonic deformation of Turkey is controlled by the collision of the Arabian, African shaped by the Cretan and Cyprus Arc, Dead Sea Fault system. Göncüoğlu et al. (1997) has pointed out that the present distribution of Turkish terranes is controlled by Alpine orogenic belt. The neotectonic structures of central Anatolian region are major fault zones and intracontinental basins during since Late Miocene–Pliocene times (Dirik, 2001). The stresses between the plates are compensated by the westward movement of Anatolian microplate along North Anatolian Fault System (NAFS) that separates the Eurasian plate from the Anatolian microplate in northern Turkey, and East Anatolian Fault System (EAFS) (Barka and Cadinsky-Cade, 1988). In addition to NAFS and EAFS, second-order contractional and extensional fault zones, deforming and dividing them in a number of smaller blocks. Most of second-order structures are mostly dextral and sinistral strike-slip faults which are controlled by tectonic regime of central Anatolian region (Koçyiğit and Deveci, 2007). Inonu-Eskisehir Fault Zone (IEFZ) and Tuzgolu Fault Zone (TGFZ) compose the transitional zone between Eastern Anatolia compression zone and Western Anatolia extension zone, Ezinepazarı Fault (EZF) and the eastern flank of the Isparta Angle (OIA) and Inner Isparta Angle (IIA). The IEFZ also forms the north-norheastern boundary of the western Anatolian region. The eastern flank of the OIA and IIA are reverse V-shaped morphotectonic structures and extend to the western part of the Mediterranean region (Figure 1a). Koçyiğit (2005) has pointed out that continental extension in the OIA and IIA is multidirectional from focal mechanism solutions of earthquakes. According to GPS-Positions, counter-clockwise rotation of Anatolia region with respect to Eurasia is not limited to central Anatolian region (Reilinger et al., 2006, Aktuğ et. al, 2013).

3.

EGM08 Bouguer anomalies and ETOPO1 topography

The regional Bouguer gravity anomalies (EGM08) and surface elevation data (ETOPO1) have been compiled in global data sets. EGM08 Bouguer gravity data have been a combined satellite-ground gravity model with regular coverage on earth’s surface and include spherical

harmonic coefficients up to degree and order 2190 and order to 2159 (Pavlis et al., 2008, 2012). The gravity data have been corrected for the long wavelengths (300 km) using GRACE satellite data and expanded to its full ~10 km spatial resolution with an accuracy of ~1 to ~3 mGal over wide parts of the Earth's surface (Pavlis et al., 2008, 2012). We downloaded regional Bouguer gravity anomaly grids for the central Anatolian region derived from the EGM08 Gravitational Model by the National Geospatial-Intelligence Agency (NGA). On the other hand ETOPO 1 (Amante and Eakins, 2009) a global topographic data set with a resolution of 1 arc min were compiled for estimating regional isostatic Moho undulation. Figure 1b shows EGM08 Bouguer anomalies of the central Anatolian region and surroundings. Negative anomalies trend NE-SW and SE-NW, and cover the higher elevation regions in the most eastern, and partly western and southeast part of the study area. Thus one can conclude that negative anomalies are due to crustal roots of lower density, compared to the lithospheric mantle. Positive anomalies mostly correspond to a low topographic relief (Figure 1c) and can be interpreted in terms of crustal thinning beneath these areas. Thus considerable positive anomalies have been identified in areas of IEFZ, outer and inner Isparta angle, EZF and TGFZ cut through basins. Figure 1 4. The average estimate of effective elastic thickness The EET is directly related to the mechanical thickness, defined as a cumulative thickness of rocks that deforms elastically. An average value of EET must be known for estimating the isostatic Moho depths in addition to the mechanic behavior of the lithosphere. The wavelength-dependence of the Bouguer anomalies can be obtained from the correlation between Bouguer anomalies and topography data. The coherence function can be used to correlate both data sets as a function of wavelength.

The spectral coherence from Bouguer anomalies and topography data is defined as

C 2 (k )  (k )  E0 (k ) E1 (k ) 2

(1)

where E1 and E0 is the power of Bouguer anomalies and topography data respectively. C is the cross power between both data (Forsyth, 1985). The theoretical form of Equation (1) based on the ratio (f) of surface to subsurface loading is given by Forsyth (1985):

2  Theo (k ) 

(1  ( f  ) 2  ) 2 . (1  ( f  ) 2 )(1  f 2 2 )

(2)

Forsyth (1985) showed that the loading ratio is related to many geological processes and there is no significant correlation when averaged over a wide region. The other parameters





and

is given as   1  D(2k ) 4  m g and   1  D(2k ) 4 ( m  c ) g . ρc and ρm are densities

of the crust and lithospheric mantle, respectively..

Thompson (2000) developed a method based on the multitaper spectral analysis in order to compute Fourier transform of the Bouguer gravity and topography data for spectral coherence function. The most important property of the multitaper spectral analysis is an estimate of the overall variance that reduces in the calculation of power spectral density by incorporating different tapers in the data sets. The method also finds ideal bandwidth of the taper window of the central lobe (W) by minimizing the spectral leakage. In that case, the choice of W is important since it influences the resulting power spectra. An ideal bandwidth is given as NW and, where N is the number of samples within the data window (Percival and Walden, 1993). Simons et al. (2000) have suggested that an appropriate value for NW can be found by comparing the results between the observed and predicted spectra for different NW values. Figure 3 shows the examples of the estimation of EET obtained from EGM08 BouguerETOPO1 topography coherence in the central Anatolia region. The method works in iteration steps to obtain a series of predicted coherence values for various EET values such as 0-50 km. It was assumed that the f is 1. Thus the subsurface load which is produced by abnormal mantle beneath the region is equal the surface loads (topography). Figure 2a, Figure 2b and Figure 2c shows the predicted curves which are compared with the observed curve, and the EET value that minimizes the root mean square error (RMS) is estimated for NW=2, NW=3 and NW=4. Both coherences as a function of wavenumber based on different EET value are shown in Figure 2d, Figure 2e and Figure 2f. Figure 2f displays the best fitting model by RMS minimization for NW=4 and EET is 21.3 km.

It is clear, from the results in Figure 2f, that the results of the coherence analysis fit a flexural model with EET=21.3. Accordingly low elastic thickness value characterizes a weak and thinned lithosphere in the region. In addition, it is understood that the surface loads create compensation roots in a manner of regional sense.

Figure 2

5. Estimating the isostatic Moho depths In regional isostasy model, the loads of the compensating masses are distributed entirely within the crust. Thus, compensation is achieved with varying thicknesses or depths of Moho interface. Indeed, a flexible crust deforms due to the elasticity of rocks that supports shortwavelength topographical masses and forms the shape of the flexure. However, longwavelength lateral variations of Moho undulations are originated from surface displacements with long wavelengths and/or ductile flow in hot lower crust. The Moho interface is deflected downward a distance w by an applied vertical load such as excess topography (Figure 3). Two parameters such as flexural rigidity and effective elastic thickness are effective in deflecting the Moho interface. Flexural bending (w) in the x-y plane is described by Nadai and Hodge (1963):

D4 w( x, y)  g ( m  c )w( x, y)  L( x, y)

(3)

where the g is the acceleration of gravity (~9.81 m/s2), D is the flexural rigidity of the lithosphere and is defined as D  E  EET 3 [12(1   2 )] where E and  are Young's modulus and Poisson ratio (~0.25), respectively. Equation (3) also means that the sum of bending resistance and buoyancy forces equals to load L. Figure 3 The exact solution to Equation (3) for a point load was discussed by Hertz (1895) and is expressed in terms of Bessel-Kelvin functions (Gruninger, 1990; Abd-Elmotaal, 1993). On the other hand, Vening-Meinesz (1940) introduced a method based on two parameters such as maximum deflection and degree of regionality. The maximum deflection w 0 is the distance directed from the compensation level to the Moho under the load at origin and written as w0 

mL 8(  m   c ) 2

(4)

where mL denotes the mass of the load generated by the topographic masses, and ℓ is the socalled “degree of regionality” (Vening Meinesz, 1940)

4

D g (m  c )

(5)

Gruninger (1990) and Abd-Elmotaal (1993, 1995a) have showed that the difference between the exact solution and the approximate solution can be neglected in all practical applications. Vening Meinesz (1940) introduced a simple method based on the Moho deflections at a distance r from the load at the origin in order to obtain an approximate solution of Equation (3). The solution follows the polynomial equations, as follows 6

4

2

w1 r r r  c1    c2    c3    c4 w0    4

0

r 2 

(6)

2

w2 r r  c5    c 6    c 7 w0  

0

r  rmax 

(7)

where rmax is the radius of curvature at compensation level or maximum extension of the regional compensation for one single load/relief and its formulation (Vening Meinesz, 1940) is rmax=2.905  . Hence the vertical deflections of the bending by Equations (6) and (7) are

w  w1  w2

(8)

The coefficients of the polynomials in Equations (6) and (7) are given by Vening-Meinesz (1940), Moritz (1990) and Abd-El Motaal (1991, 1993). For determining the Vening-Meinesz compensation depths, the mass mL due to a vertical column of the topography of density ρC, height h and infinity small area ds, can be written as:

mL   c hds

(9)

Thus the bending due to topographic column (Moritz, 1990) will be

ZmL  w(r ) c hds

(10)

The bending due to the entire masses can be given by Moritz (1990) and Abd-Elmotaal (1995a) w( x, y)    c h( x , y ) w(r )ds s

The total bending is closely related to the Moho depth, as shown in Figure 3

(11)

hM  T0  w

(12)

where T0 is the crustal thickness for a topography with respect to sea level (zero heights).

The topographic loads are represented by the ETOPO1 surface elevation data in Figure 1b. A density contrast of 630 kg/m3 is assumed according to standard values of densities (ρ m=3300 kg/m3, ρc=2670 kg/m3). The coise of T0 also plays a key role in mapping Moho undulations rather than EET. Although there is no predictive value of T0 in central Anatolian region, some studies provide an approximation. For instance, Pamukçu and Yurdakul (2008) have found the value of 33 km from the Bouguer-topography admittance analysis in western Anatolia that partly covers central Anatolia. Vanacore et al. (2013) have showed that Moho depths increase from west to east from less than 30 km to greater than 50 km beneath Eastern Anatolian’s higher elevations. As a result, the compensation depth (normal crustal thickness) in central Anatolian region was taken as 35 km. On the other hand, it is difficult to determine the degree of regionality (in km) due to the heterogeneous structure of the lithosphere. This parameter has a length dimension and its values (in S.I. units) ranges from 10 to 60 km. Small values of this parameter causes Moho surface with minima more narrow and deeper, the opposite cause wide and shallow minima (Moritz, 1990; Corchete et al., 2010). In the calculations the value of l=10 was used to map the narrow and deeper gradient zones of the Moho undulation. The EET value which is necessary to calculate the flexural rigidity (D) and hence isostatic Moho depths, was taken as 21.3 km from Bouguer-Topography coherence. The isostatic Moho depths were mapped as an inverse reflection of the topography (Figure 4). The geometry of the Moho undulation is generally smooth with an average depth of 40 km and has long-wavelength undulation beneath the central part of the region corresponding to the lower topography. However, the Moho undulations occur prominently along the N-S direction in the western part of the region. The maximum isostatic Moho deepening is under the Alpine Orogenic belt throughout the southern part of Anatolia, while the minimum Moho depths are beneath mostly the central part and westward displacement of the crust in the northwest part of the region (Figure 4). The northwest-southeast zone of region is a transitional zone of Moho undulation. Locations with the only significant change of isostatic Moho undulation appear beneath the zone between central and southeastern part of the region. This is because of the existence of many basins in the central Anatolia such as the Tuzgolu basin, Haymana Basin, Tersakan Basin (Aydemir, 2008; Aydemir, 2011; Ateş, 2006). This

steep belt of deepening Moho appears in both the CAFZ and the OIA. In general it is understood that the higher elevations in the region not only pushes down the crust but its surroundings from the results of isostatic Moho depths. On the other hand, the depressions induce sedimentary basins in the region that works as a negative load that compensate the higher elevations (positive loads). As shown in Figure 4b, it is clear that the flexural deformation works as a low-pass filter for the horizontal distribution of all the loads. As a result, topographic loads in the Alpine Orogenic belt throughout the southern part of Anatolia provide the deeper and narrow depressions for an elastically deformed lithosphere but wide and shallow ones may be caused by many basins or rigid lithosphere. The isostatic Moho depths (37.7-43 km) are generally consistent with previous studies including seismic experiments. Accordingly Marone et al. (2003) determined that the Moho depths in the region vary between 36 and 40 km. Kuleli et al. (2004) found that the crustal thickness varies between 37 and 42 km in the region with controlled seismic sources. According to Tezel et al. (2013), the crustal thicknesses in central Anatolia vary between 31 km and 38 km using receiver function. Vanacore et al. (2013) declared that Moho depths range from about 37 to 47 km using receiver functions. Figure 4 The statistical relationship between isostatic Moho depths and topography data was examined using a regression analysis. Simple linear regression analysis was performed in a least squares sense,

as

shown

in

Figure

5.

The

linear

regression

model

has

the

form:

Topography= 0.503×IsostaticMoho-19.041 in the central Anatolian region. A strong positive linear relationship shows that regional isostatic mechanism is generally achieved. A positive correlation coefficient (0.8497) indicates that as regional isostatic Moho depths increases the topography data tend to increase. Thus, it is implied that topography data enable us to estimate the Moho depths within the study area in a manner of regional isostasy concept. Figure 5

5. Linear inversion for estimating density contrast distribution in the uppermost crust In this section, the density contrast distribution within the uppermost part of crust was infered. A regional-residual separation process was applied to extract the wavelengths coming from uppermost crust in the space domain. To this end, residual gravity anomaly map based

on removing a third-order polynomial trend surface from the EGM08 Bouguer anomalies were obtained (Figure 6a). The inverse solution of potential field data is non-unique and unstable. Thus the problem is defined as ill-posed. A way to reduce the instability and to guarantee the uniqueness of the inverse solution is to introduce a priori knowledge in relation to the unknown model parameters. Indeed, a linear gravity inversion has difficulties in constructing the anomalous sources in the number of iterations. Last and Kubik (1983) introduced a concept of compactness for the causative bodies that limit the lower and upper bound of the density contrasts. With matrix notation, the parameter vector for model parameters that represents density contrast of the blocks is rewritten as: 1

  m2  1 T  dg P  W G  GWm G  1   e2   1 m

T

(14)

where dg is the data vector whose elements are the observed gravity data, G is the system matrix whose jth element contains the kernel function of the jth block on the ith observation 1 2 point, Wm is the parameter weighting matrix (the density contrast of each block),  m is the 2 variance of the density contrasts, and  e is the variance of the deviation between observed

and computed gravity anomalies. It should be noted that Equation (14) was defined as an underdetermined system. In the case of any block that crosses the density contrast barrier (Po) will be set to Po and the algorithm automatically freezes this block in the next iteration by assigning a very small weight to this block (Last and Kubik, 1983). The iterative procedure stops when a minimum area of the density contrast distribution is reached. The stopping criteria in inversion algorithm are based on the fit between the observed data and calculated data by the density contrast model. In Figure 6c, the subsurface density contrasts were estimated from the inversion of the profile data (Figure 6b) which is extracted from residual gravity anomaly map. This interpretive model was built by subdividing into 882 cells (9×98) in the inversion procedure. As a solution from the compact gravity inversion by using Equation (14), the result has the boundaries of density contrasts and converged after 100 iterations with an RMS of 0.42. It should be noted that the horizontal variations in density contrast distribution are well correlated the lineaments

in the region. The uplifts and depressions along the AB profile have been sequentially mapped as an indicator of strong heterogeneities in the uppermost crust. Density contrasts up to ±0.5 g/cm3 in the subsurface were mostly imaged the boundaries between uplifts and depressions. The central part of the region displays relatively high (positive) density contrasts, while low density contrasts are imaged in the zone between IIA and OIA and, which might be resulted from the sedimentary basin with lower densities than those of surrounding rocks.

The

problematic areas from the aspect of rigidity within the uppermost crust are those small fields comprising low density contrasts.

Figure 6 6.

Conclusion

In this study, flexural isostasy model were applied to estimate the Moho depths for the first time in the central Anatolian region. It was found that the isostatic state is primarily controlled by young geological provinces characterized by low EET (21.3 km). We conclude that partial melting probably occurs throughout the lithospheric mantle in addition to a relatively shallow lithosphere-asthenosphere boundary (LAB). The isostatic Moho depths can be up 43 km in the southeast part of the central Anatolian region with a higher elevation near the CAFZ. The Moho depths are well correlated with topography data according to a strong positive correlation of regression analysis. The correlation reveals that thin crust is represented in the low-relief plateau whereas thicker crust is related to the high-relief plateau. Thus, one can conclude that the higher elevation is isostatically supported by deeper Moho depths. However, the isostatic conditions may be potentially affected by uppermost crustal heterogeneity in terms of lateral density contrast variations. Thermal structures in the lithosphere are also crucial to investigate isostasy dynamics. On this basis, the rheology of the lithosphere is of primary importance for lithospheric strength and its relation to seismic activities. It should be focused on LAB undulation, thermal modeling, and spatial variations of effective elastic thicknesses of the lithosphere, rheological modeling and regional stress tensor using focal mechanisms of earthquakes in order to better understand dynamic stresses in the region.

Acknowledgements

We are grateful to David Gomez-Ortiz for providing the spectral coherence programs used in this study. We thank the anonymous reviewers for their insightful comments and suggestions. This study was supported by the Scientific and Technological Research Council of Turkey (TUBITAK), Project Number 115Y217.

References Abd-Elmotaal, H., 1991. Gravity anomalies based on the Vening Meinesz isostatic model and their statistical behaviour. Mitteilungen der geoda¨tischen Institute der Technischen Universitat Graz, Folge 72, Graz Technical University

Abd-Elmotaal, H., 1993. Vening Meinesz Moho depths: traditional,

exact

and

approximated. Manuscripta Geodaetica, 18, 171-171.

Abd-Elmotaal H., 1995a. Theoretical background of the Vening Meinesz isostatic model. In: Su¨nkel H, Marson I (eds) Gravity and geoid. IAG Symposium 113. Springer, Berlin Heidelberg New York, 268–277

Abd-Elmotaal, H. A., 2004. Isostatic response of the earth's crust derived by inverse isostasy. Journal of Geodynamics, 37, 139-153.

Airy, G.P.1855. On the computation of the effect of the attraction of the mountain masses as disturbing the apparent astronomical latitude of stations in geodetic surveys. Philosophical Transactions of the Royal Society of London, 145, 101-104. Aktuğ, B., Parmaksız, E., Kurt, M., Lenk, O., Kılıçoğlu, A., Gürdal, M.A., Özdemir, S., 2013. Deformation of central Anatolia: GPS implications. Journal of Geodynamics, 67, 78-96.

Amante, C., Eakis, B.W., 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M. Ateş, A., Bilim, F., Büyüksaraç, A., 2005. Curie point depth investigation of central Anatolia, Turkey. Pure and Applied Geophysics, 162, 357-371. Ateş, A., Bilim, F., Büyüksaraç, A., Aydemir, A., Bektaş, O., Aslan, Y., 2012. Crustal structure of the Turkey from aeromagnetic, gravity and deep seismic reflection data. Surv. Geophys. 33, 869-885. Aydemir, A., Ateş, A., 2006. Structural interpretation of the Tuzgolu and Haymana Basins, Central Anatolia, Turkey, using seismic, gravity and aeromagnetic data. Earth, Planets and Space, 58, 951-961.

Aydemir, A., 2008. Hydrocarbon potential of the Tuzgolu (Salt Lake) Basin, Central Anatolia, Turkey: A comparison of geophysical investigation results with the geochemical data. Journal of Petroleum Science and Engineering. 61, 33-47.

Aydemir, A., 2009. Tectonic investigation of Central Anatolia, Turkey, using geophysical data. Journal of Applied Geophysics, 68, 321-334.

Aydemir, A, 2011. An integrated geophysical investigation of Haymana Basin and hydrocarbon prospective Kirkkavak Formation in Central Anatolia, Turkey. Petroleum Geoscience, 17, 91-100.

Banks, R. J., Parker, R. L., Huestis, S. P. 1977. Isostatic compensation on a continental scale: local versus regional mechanisms. Geophysical Journal International 51(2), 431-452.

Barka, A. A., Kadinsky‐ Cade, K., 1988. Strike‐ slip fault geometry in Turkey and its influence on earthquake activity. Tectonics, 7, 663-684.

Beaumont, C., Muñoz, J. A., Hamilton, J., Fullsack, P., 2000. Factors controlling the Alpine evolution of the central Pyrenees inferred from a comparison of observations and geodynamical models. Journal of Geophysical Research: Solid Earth. 105, 8121-8145.

Bechtel, T.D., Forsyth, D.W., Sharpton, V.L., Grieve, R.A.F., 1990. Variations in effective elastic thickness of North American lithosphere. Nature, 343, 636-638.

Bilim, F.; Aydemir, A.; Ates, A., 2015. Determination of block rotations and the Curie Point Depths of magnetic sources along the NW-SE-trending Suluklu-Cihanbeyli-Goloren and Sereflikochisar-Aksaray Fault Zones, Central Anatolia, Turkey. Geodinamica Acta, 27, 203213

Bozkurt, E., 2001. Neotectonics of Turkey–a synthesis. Geodinamica acta, 14, 3-30. Büyüksaraç, A., Jordanova, D., Ateş, A., Karloukovski, V., 2005. Interpretation of the gravity and magnetic anomalies of the Cappadocia Region, Central Turkey. Pure and Applied Geophysics 162, 2197–2213.

Corchete, V., Chourak, M., Khattach, D., 2010. A methodology for filtering and inversion of gravity data: an example of application to the determination of the Moho undulation in Morocco. Engineering 2, 149.

Dirik, K. 2001. Neotectonic evolution of the northwestward arched segment of the Central Anatolian Fault Zone, Central Anatolia, Turkey. Geodinamica Acta 14, 147–158.

Flesch, L. M., Haines, A. J., Holt, W. E., 2001. Dynamics of the India‐ Eurasia collision zone. Journal of Geophysical Research: Solid Earth 106, 16435-16460.

Forsyth, D. W. 1985. Subsurface loading and estimates of the flexural rigidity of continental lithosphere. Journal of Geophysical Research: Solid Earth 90, 12623-12632.

Gomez-Ortiz, D., Tejero, R., Ruiz, J., Babın-Vich, R., Gonzalez-Casado, J.M., 2005. Estimating the effective elastic thickness of the Iberian Peninsula’s lithosphere based on multitaper spectral analysis. Geophys. J. Int. 160, 729–735.

Grobys, J. W., Gohl, K., Eagles, G., 2008. Quantitative tectonic reconstructions of Zealandia based on crustal thickness estimates. Geochemistry, Geophysics, Geosystems 9.

Gruninger W. 1990. Zur topographisch-isostatischen Reduktion der Schwere. Dissertation an der Universita¨t Karlsruhe, Karlsruhe Technical University

Gürsoy H., Piper J.D.A., Tatar O., Temiz H. A., 1997. Paleomagnetic study of the Sivas Basin, Central Turkey: crustal deformation during lateral extrusion of the Anatolian block, Tectonophysics 271, 89–105.

Hertz, H., 1895. Gesammelte werke, vol. 1. Johann Ambrosius Barth, Leipzig, 155.

Koçyigit, A., 2003. Neotectonic characteristics and Seismicity of Central Anatolia. Türkiye Petrol Jeologları Derneği Bülteni (special issue) 1–26.

Kuleli, S., Tuksoz, M. Gurbuz, C., Gok, R., Schultz, C., 2004. Crustal structure study in Turkey with controlled seismic sources, EOS, Trans.Am.Geophys.Un., 85 (47), Fall Meet.Suppl., Abstract S13B-1058.

Last, B. J., Kubik, K., 1983. Compact gravity inversion. Geophysics 48, 713-721.

Lewis, B. T., Dorman, L. M. 1970. Experimental isostasy: 2. An isostatic model for the USA derived from gravity and topographic data. Journal of Geophysical Research 75(17), 33673386. Marone, F., van der Meijde, M., van der Lee, S., Giardini, D., 2003. Joint Inversion of Local, Regional and Teleseismic Data for Crustal Thickness in the Eurasia-Africa Plate Boundary Region. Geophysical Journal International 154, 499-514.

McKenzie, D., Bowin, C. 1976. The relationship between bathymetry and gravity in the Atlantic Ocean. Journal of Geophysical Research 81, 1903-1915.

Macario, A., Malinverno, A. and Haxby, W.F., 1995. On the robustness of elastic thickness estimates obtained using the coherence method. Journal of Geophysical Research, 100, 15163- 15172.

McNutt, M. K., Parker, R. L. 1978. Isostasy in Australia and the evolution of the compensation mechanism. Science, 199, 773-775.

McNutt, M. K. 1983. Influence of plate subduction on isostatic compensation in northern California. Tectonics, 2, 399-415.

Moritz, H. 1990. The inverse Vening Meinesz problem in isostasy. Geophysical Journal International 102, 733-738. Nadai, A., Hodge, P. G., 1963. Theory of Flow and Fracture of Solids, vol. II. Journal of Applied Mechanics 30, 640.

Naliboff, J. B., Lithgow-Bertelloni, C., Ruff, L. J., de Koker, N. 2012. The effects of lithospheric thickness and density structure on Earth′s stress field. Geophysical Journal International 188, 1-17.

Oruç, B., Sönmez, T., 2017. The rheological structure of the lithosphere in the eastern Marmara Region, Turkey. Journal of Asian Earth Sciences, 139, 183-191.

Oruç, B., Gomez-Ortiz, D., Petit, C. 2017. Lithospheric flexural strength and effective elastic thicknesses of the Eastern Anatolia (Turkey) and surrounding region. Journal of Asian Earth Sciences 150, 1-13.

Pamukçu, O., Yurdakul, A. 2008. Isostatic compensation in western Anatolia with estimate of the effective elastic thickness. Turkish Journal of Earth Sciences, 17(3), 545-557.

Pamukçu, O. A., Akçığ Z., 2011, Isostasy of Eastern Anatolia (Turkey) and discontinuities of its crust: Pure and Applied Geophysics 168, 901–917.

Pavlis, N. K., Holmes, S. A., Kenyon, S. C., Factor, J. K., 2008. An earth gravitational model to degree 2160: EGM2008. EGU General Assembly, 2008 (4), 4-2.

Percival, D. B., Walden, A. T., 1993. Spectral analysis for physical applications. Cambridge University press.

Pfiffner, O. A., Ellis, S., Beaumont, C., 2000. Collision tectonics in the Swiss Alps: Insight from geodynamic modeling. Tectonics 19, 1065-1094.

Reilinger, R. E., McClusky, S. C., Oral, M. B., King, R. W., Toksoz, M. N., Barka, A. A., Sanli, I., 1997. Global Positioning System measurements of present‐ day crustal movements in the Arabia‐ Africa‐ Eurasia plate collision zone. Journal of Geophysical Research: Solid Earth 102, 9983-9999.

Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., Nadariya, M., 2006. GPS constraints on continental deformation in the Africa‐ Arabia‐ Eurasia continental collision zone and implications for the dynamics of plate interactions. Journal of Geophysical Research: Solid Earth, 111.

Simons, F. J., Zuber, M. T., Korenaga, J., 2000. Isostatic response of the Australian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper spectral analysis. Journal of Geophysical Research: Solid Earth 105, 19163-19184. Şengör A. M. C., Yılmaz, Y., 1981. Tethyan evolution of Turkey: a plate tectonic approach. Tectonophysics 75, 181193203-190199241.

Tezel, T., Shibutani, T., Kaypak, B., 2013. Crustal thickness of Turkey determined by receiver function. J Asian Earth Sci 75, 36-45.

Thompson, D. J., 2000. Multitaper analysis of nonstationary and nonlinear time series data. Nonlinear and nonstationary signal processing 317-394.

Valera, J. L., Negredo, A. M., Jiménez-Munt, I., 2011. Deep and near-surface consequences of root removal by asymmetric continental delamination. Tectonophysics 502, 257-265.

Vanacore, E.A., Taymaz, T., Saygin, E., 2013. Moho structure of the Anatolian Plate from receiver function analysis. Geophys. J. Int., 193, 329-337.

Vening Meinesz F.A.,1931. Une nouvelle methode pour la reduction isostatique regionale de l’intensite de la pesanteur Bull. Geod 29, 33-51.

Vening Meinesz F. A., 1940. The earth's crust deformation in the East Indies. Koninklijke Nederlandsche Akademie van wetenschappen.

Walcott, R. I., 1970. Isostatic response to loading of the crust in Canada. Canadian Journal of Earth Sciences, 7(2), 716-727.

Watts, A. B., Cochran, J. R., 1974. Gravity anomalies and flexure of the lithosphere along the Hawaiian‐ Emperor seamount chain. Geophysic al Journal International, 38(1), 119-141.

Watts, A. B., Bodine, J. H., Steckler, M. S. 1980. Observations of flexure and the state of stress in the oceanic lithosphere. Journal of Geophysical Research: Solid Earth, 85(B11), 6369-6376.

Watts, A. B., 2001. Isostasy and Flexure of the Lithosphere. Cambridge University Press.

Zuber, M.T., Bechtel, T.D. and Forsyth, D.A., 1989. Effective elastic thickness of the lithosphere and mechanisms of isostatic compensation in Australia. J. Geophys. Res., 94, 9353-9367.

Figure Captions

Figure 1. a) Major tectonic structures in the Central Anatolian region and fault traces (modified after Koçyiğit, 2003). b) Surface elevation data from ETOPO1 (Amante and Eakis, 2009). c) Bouguer anomaly data compiled from the Earth Gravitational Model (EGM08) (Pavlis et al., 2008).

Figure 2. The estimation of an average EET from Bouguer-topography coherence for different NW values. (a-c) The misfit curves that have well-marked minimum, are used to estimate the EET. (d-f) Observed and predicted isotropic Bouguer-topography coherence curves. Note that the best fitting is chosen by error minimization and obtained for NW=4 and EET=21.3 km.

Figure 3. Moho deflection under the effect of the load (L) of the topography. EET: Effective elastic thickness, T0: Compensation depth, wo: Maximum deflection, w: Moho deflection at a distance r (modified after Abd-Elmotaal, 2004).

Figure 4. a) A detailed topography from ETOPO1 digital elevation model in central Anatolia region. b) Regional isostatic Moho depths to central Anatolian region.

Figure 5. A plot of isostatic Moho depths versus surface elevation data of the central Anatolia region. Note that scatter plot shows strong positive linear correlation.

Figure 6. a) Residual gravity anomaly map based on removing a third-order trend surface from the EGM08 Bouguer anomalies in Figure 1b. b) AB profile from residual gravity anomalies. c) Interpretative block model used in the inverse solution. The subsurface was subdivided into 882 cells (9×98) by setting minimum and maximum density contrasts of ±0.5 g/cm3. RMS of error is found as 0.42 in 100 iterations.

Graphical abstract

Highlights

 Spectral coherence estimates provide an average effective elastic thickness value of 21.3 km.  Lithospheric strength may be controlled primarily by young geological provinces and partial melting of the lithospheric mantle.  The Moho undulation is generally smooth with an average depth of 40 km.

 A strong positive linear relationship between regional isostatic Moho and topography shows that regional isostatic mechanism is generally achieved.  The isostatic conditions may be potentially affected by uppermost crustal density contrast

variations.