Crustal Density Model of the Sea of Marmara: Geophysical Data Integration and 3D Gravity Modelling

The Sea of Marmara, in Northwest Turkey, is a transition zone where the dextral North Anatolian Fault Zone (NAFZ) propagates westward from the Anatolian plate to the Aegean plate. The area is of interest in the context of seismic 10 hazard in the vicinity of Istanbul, a metropolitan area with about 15 million inhabitants. Geophysical observations indicate that the crust is heterogeneous beneath the Marmara Basin, but a detailed characterization of the crustal heterogeneities is still missing. To assess if and how crustal heterogeneities are related to the NAFZ segmentation below the Marmara Sea, we develop a new crustal-scale 3D density model which integrates geological and seismological data and is additionally constrained by 3D gravity modelling. This model indicates that the observed gravitational anomalies originate from significant 15 density heterogeneities within the crust. Two layers of sediments, one syn-kinematic and one pre-kinematic with respect to the Marmara Sea formation are underlain by a heterogeneous crystalline crust. A felsic upper crystalline crust (average density of 2720 kg.m-3) and an intermediate to mafic lower crystalline crust (average density of 2890 kg.m-3) appear to be crosscut by two large, dome-shaped mafic high-density bodies (average density of 3050 kg.m-3) of considerable thickness above a rather uniform lithospheric mantle (3300 kg.m-3). The spatial correlation between the bent segments of the fault and the location of 20 the high-density bodies suggests that the distribution of lithological heterogeneities within the crust controls the rheological behaviour along the NAFZ, and consequently, influences fault segmentation and propagation dynamics.


Introduction
The Sea of Marmara in NW Turkey is an extensional basin associated with a right-stepping jog in the orientation of the North Anatolian Fault Zone (NAFZ; Fig. 1), a westward-propagating right-lateral strike-slip fault that constitutes the plate boundary between the Anatolian and the Eurasian plates (Fig. 1a;McKenzie, 1972;Şengör et al., 2005).As one of the most active platebounding strike-slip faults in the world, and being located in the Istanbul metropolitan area with a population of approximately al., 2017a: https://www.gonaf-network.org;MARsite: http://www.marsite.eu)have been embarked on to improve the observational basis for hazard and risk assessments in the area of the Marmara Sea.
The Marmara section of the NAFZ, considered to be a 150-km-long seismic gap between the ruptures of two big events in 1912 (M 7.3) and 1999a (M 7.4), is a zone of strong earthquakes (M ~ 7.4) with a recurrence time of approximately 250 years (Fig. 1b and 1c); this section experienced the last ground-rupturing events in 1509 and 1766, suggesting that the fault is mature and a ground-rupturing event might be expected (Ambraseys, 2002;Barka et al., 2002;Parsons, 2004;Janssen et al. 2009;Murru et al., 2016;Bohnhoff et al., 2016a;2016b;2017b).There, the potential for a large seismic event is regarded as being high (Bohnhoff et al., 2013;2017b).A key question is if this 150-km-long seismic gap will rupture in the future in one event or in several separate events due to the segmentation, an issue that will depend a lot on the stress evolution along strike among other forcing factors.In this regard, three-dimensional (3D) geological models are the fundament of geomechanical models, and the distribution of density is of key importance as density controls body forces.Density modelling is generally done by integrating geological information, seismic observations, and gravity data.In the Marmara region, however, 3D gravity modelling has not been at the focus compared with the use of other geophysical methods.Furthermore, gravity models can also help clarify the density distribution at greater depth where borehole observations and/or seismic surveys have limitations.
Our study aims at assessing the deep crustal configuration of the Marmara Sea and surrounding areas.To address the question if there is a spatial relationship between fault activity and the distribution of certain physical properties in the crust, we developed a 3D density model that integrates available seismological observations and is consistent with observed gravity measurements.In a previous gravity modelling effort (Kende et al., 2015), after wavelength-depth related crustal correction, an inversion method was applied to calculated the Moho depth below the Marmara region.This study (Kende et al., 2015) revealed that assuming laterally homogeneous crustal layers is insufficient to integrate density heterogeneities.Building on an earlier 3D structural model (sediments; crust; Moho; Fig. 2 and 3) developed to evaluate the stress-strain state in this region (Hergert and Heidbach, 2010;2011;Hergert et al., 2011), we used crustal and regional-scale forward 3D gravity modelling and seismic data as additional constraints.Here we show that significant density heterogeneities are laterally present within the crust below the Sea of Marmara.
Our results from forward 3D gravity modelling point to the existence of lateral density heterogeneities within the crust in form of two local high-density bodies.This supports the notion of rheologically strong crustal domains that may influence the kinematics of the NAFZ below the Sea of Marmara.In addition, considering the results of the earlier geomechanical-numerical models (Hergert and Heidbach, 2010;2011;Hergert et al., 2011), we suggest that implementing the high-density bodies in the structural model may change earlier results.

Geological setting
In the large-scale plate-tectonic framework of Asia Minor, the NAFZ accommodates the westward escape of the Anatolian plate in response to the northward motion and indentation of the Arabian plate into Eurasia (Fig. 1a); this has resulted in numerous deformation features along the well-defined trace of the fault and regionally, along the northern flanks of the Anatolian Plateau (Barka and Hancock, 1984;Barka, and Reilinger, 1997;Pucci et al., 2006;Yildirim et al. 2011;2013).
In the westernmost sector, the NAFZ bifurcates into several strands of locally variable strikes, which has resulted in a mosaic of pull-apart basins flanked by steep mountain fronts and intervening push-ups; all of these morphotectonic features of the greater Marmara region are characterized by active Quaternary deformation process (Yildirim and Tüysüz, 2017).To the west of the Almacik Block, a transpressional push-up ridge, the NAFZ splits into three main strands (Fig. 1b; Armijo et al., 1999Armijo et al., , 2002)): the northern, middle, and southern branches.The northern branch traverses the Sea of Marmara and forms the N70°E striking Main Marmara Fault (MMF; Le Pichon et al., 2001;2003).The approximately E-W-striking middle branch passes through the Armutlu Peninsula and continues along the southern coast of the Marmara Sea; this branch changes strike to NE-SW in the southern part of the Kapıdağ Peninsula (Yaltırak and Alpar, 2002;Kurtuluş and Canbay, 2007).The southern branch traverses the Biga Peninsula, the region to the south of the southern margin of the Marmara Sea.
The Marmara Sea is an E-W elongated transtensional basin with up to 1300 m water depth along its axial part; it is surrounded by onshore domains at about 600 m average elevation (Fig. 1c).The deepest part of the basin is the North Marmara Trough (NMT; Bécel et al. 2009), which hosts three main sedimentary basins along the NAFZ.These include the Çınarcık, Central, and Tekirdağ basins.These depocentres are separated from each other by the shallower Central High (East) and the Western High (West), respectively.In the deep parts of the basin protracted subsidence has resulted in the accumulation of more than 5 km of Pliocene-Holocene sediments (Le Pichon et al. 2001;2003;2015;Armijo et al. 2002;Parke et al. 2002;Carton et al., 2007;Laigle et al., 2008;Bécel et al. 2009;2010;Bayrakci et al., 2013).
The region of the Marmara Sea is an integral part of the NAFZ, which began its activity in the east approximately 13 to 11 Ma ago (Şengör et al., 2005).Although different models and timing constraints for the onset of basin formation in the Marmara Sea have been presented in the context of the evolution of the NAFZ and the Aegean region (e.g., Armijo et al., 1999;Ünay et al., 2001;Yaltırak, 2002;Şengör et al., 2005;Le Pichon et al., 2014;2015), offset geological marker horizons, displace structures, and paleontological data point to a transtensional origin during the propagation and sustained movement of the NAFZ with displacement and block rotations after the Zanclean transgression.Such a geodynamic scenario of transtensional dextral strike-slip faulting is compatible with space-geodetic data, the pattern of seismicity and geomorphic indicators in the landscape (Reilinger et al., 1997;2006;Barka and Kadinsky-Cade, 1988;Bürgmann et al., 2002;Pucci et al., 2006;Akbayram et al., 2016;Yildirim and Tüysüz, 2017).

Method and model setup
Like for the earlier 3D model (Hergert and Heidbach, 2010), our study area extends from 40.25° N-27.25°E to 41.15° N-30.20°E and is projected as a rectangular shape in WGS84 UTM Zone 35N with a dimension of 250-by-100 km (Fig. 1c).It covers the Sea of Marmara and the adjacent onshore areas, as well as the city of Istanbul and the Bosporus.The principal approach used for this study is crustal-scale 3D gravity forward modelling to assess the density configuration of different structural units.In this methodology, the gravity response of a model is calculated and compared with the observed gravity field.The model is iteratively modified to find the best fit with observations.Since the solution is not unique in gravity modelling, it is required to reduce the number of free parameters by integrating other available geophysical and/or geological data as additional constraints.In the spirit of this philosophy, the workflow adopted in this study consists in: (1) setting up an initial density model (Fig. 2 and 3) -in our case based on the previous studies (Hergert and Heidbach, 2010;2011;Hergert et al., 2011); (2) calculating the gravity response of this initial model and analysing the misfit (gravity residual) between modelled and observed gravity; (3) modifying the initial model by introducing additional density variations while integrating additional constraining data to obtain the density-geometry configuration that reproduces the observed gravity field best.In general, positive residual anomalies indicate that more mass is required in the model to fit the observed gravity field, whereas negative residuals imply that the mass in the model is too large in the domain of the misfit.3D forward gravity modelling has been performed using the Interactive Gravity and Magnetic Application System-IGMAS+ (Transinsight GmbH©;Schmidt et al, 2011).In IGMAS+, the gravity response of a 3D structural and density model is calculated and compared with the observed gravity field over the model area.Therefore, the model has to be defined in terms of the geometric configuration of its individual structural units.In addition to geometry, information on the densities needs to be assigned to the different units of the model to calculate the gravity response.The chosen parameter combinations for the different models studied are detailed in Sect. 4. IGMAS+ provides the density-geometry configuration in the form of triangulated polyhedrons over the 3D model domain.These polyhedrons are spanned between 2D vertical working sections where the model can be interactively modified (Schmidt et al., 2011).For this study, a lateral resolution of 2500 meter is considered that results in 100 North-South oriented working sections.Downward the model extends to a constant depth of 50 km b.s.l and the unit comprised between the Moho and the lower model boundary is considered as the uniform lithospheric mantle.
Key horizons where major contrasts in density are expected are the air-water interface, the sediment-water interface, the interface separating sediments and crystalline crust and the crust-mantle boundary (Moho).These interfaces also are well imaged with seismic methods and can therefore easily be integrated.Internal heterogeneities within the crust, may not be identified by seismic methods or only locally along individual profiles.This is where 3D gravity modelling can be used in addition to translate velocities to densities first along the seismic section and use density modelling to close the gaps in between.
This strategy together with the three-dimensionality of the calculation strongly reduces the non-uniqueness of gravity modelling as densities need to be in certain ranges for different rock types and density anomalies at different depths produce gravity effects of different wavelengths (e.g.Schmidt et al., 2011;Maystrenko et al., 2013;Sippel et al., 2013;Maystrenko and Scheck-Wenderoth, 2013).
To assess the density variations in the deeper crust of the Marmara Sea region, we calculate the gravity response for models of increasing complexity concerning their 3D structural and density configuration: (1) the initial model with homogeneous crust below the sediments, (2) a more differentiated model integrating additional seismic observations for the different crustal Throughout the modelling procedure, the uppermost surface, the bathymetry (Fig. 1c), the top-basement depth (Fig. 2a) and the depth to the Moho discontinuity (Fig. 2b) are kept fixed as defined in the initial model since the geometries of these interfaces are well-constrained by geological and geophysical data.In all tested models an average density of 1025 kg.m -3 was assumed for seawater, and a homogeneous density of 3300 kg.m -3 is assigned to the mantle below the Moho.For all gravity models presented, we define the uppermost surface of the model as the onshore topography and as the sea level offshore.
Accordingly, the thickness between sea level and bathymetry (Fig. 1c) corresponds to the column of seawater (Fig. 3a) which attains largest values in the Tekirdağ, Central, and Çınarcık basins.

Input data
The database for this study includes topography-bathymetry data, geometrical and density information from a previous 3D structural model, seismic observations, and satellite gravity data.

Topography and bathymetry
The topography-bathymetry (Fig. 1c) was exported from 1 Arc-Minute Global Relief Model (ETOPO1; Amante and Eakins, 2009).This dataset, over the study area, integrates the 30 arc-second grid obtained from NASA's Shuttle Radar Topography Mission (SRTM) and a bathymetry dataset (MediMap Group, 2005) with 1 km resolution.Figure 1c illustrates that the presentday Marmara Sea is surrounded by up to 1500 m high regions.The configuration of the present-day sea floor shows that the Marmara Sea is structured into the three main depocentres of the Tekirdağ Basin, the Central Basin, and the Çınarcık Basin where the water depth reaches up to 1300 m.While the axis of the Central Basin is aligned along the MMF, the Çınarcık Basin and the Tekirdağ Basin extend only south and mostly north of the MMF, respectively.The MMF bends along the northern boundary of the Çınarcık Basin, at the Tuzla Bend, from an E-W directed strike (East of the Marmara Sea) to an ESE-WNW strike direction at the north-western margin of the Çınarcık Basin before it resumes the E-W strike direction at the Istanbul Bend.The segment of the MMF between the two bends is the Princes Islands Segment.Farther in the West of the Marmara Sea, at the Ganos Bend, the MMF once more changes strike direction from E-W to ENE-WSW.There, the MMF exits the Sea of Marmara and creates the Ganos Fault segment of the NAFZ.

Initial model
The 3D structural model (Hergert and Heidbach, 2010;Hergert et al., 2011), considered to be the initial model for our study, differentiates three main horizons: (1) the topography-bathymetry (Fig. 1c), (2) a top-basement surface (Fig. 2a), and (3) the Moho discontinuity (Fig. 2b).In their study, Hergert and Heidbach (2010), modeled the top-basement geometry based on seismic observations (Parke et al., 2002;Carton et al., 2007;Laigle et al., 2008;Bécel et al., 2009;2010) and other geophysical and geological data such as 3D seismic tomography (Bayrakci, 2009), well data (Ergün and Özel, 1995;Elmas, 2003) and geological maps (Elmas and Yigitbas, 2001).This surface, however, has been interpreted by others as the top of a Cretaceous limestone that is pre-kinematic with respect to the opening of the Sea of Marmara (Ergün and Özel, 1995;Parke et al., 2002;Le Pichon et al., 2014).Hergert and Heidbach (2010) derived the thickness of the sediments of the Marmara Sea as the difference between bathymetry-topography and top-basement.Accordingly, their "basement" delineates the base of the sediments and not the crystalline basement.We, therefore, regard these sediments as the "syn-kinematic sediments" and refer to top-basement as the "base syn-kinematic sediments" in the following.
The thickness between the topography-bathymetry and the base syn-kinematic sediments represents the syn-kinematic sediment fill (Fig. 3b).This thickness is on average about 2.5 km over the Marmara Sea area and the unit is missing outside the Marmara Sea in response to its syn-kinematic origin.Two thickness maxima indicate localized subsidence and sediment accumulation, one aligned along the MMF where the syn-kinematic sediments are more than 5.2 km thick below the presentday Central Basin and southeastern part of the Tekirdağ Basin, and the second maximum of up to 5 km below the Çınarcık Basin limited northward by the MMF.
The depth to the Moho interface in the initial model (Fig. 2b) has been obtained by interpolating between various seismic data covering a larger area than the model area (Hergert et al., 2011, Supporting Information, Fig. S1).To constrain the Moho depth to the model area, Hergert et al. (2011) applied a Gauss filter to adjust the local variation of the Moho depth.The Moho is distinctly shallower below the Marmara Sea than below the surrounding onshore areas and shows updoming to a depth of 27 km below the basin.Along the basin margins, the Moho is about 30 km deep and descends eastwards to more than 35 km depth beneath Anatolia.

Geophysical data
The gravity observations considered in this study are based on EIGEN-6C4 (Förste et al., 2014).This dataset is a combined global gravity field model up to degree and order 2190 correlating satellite observations (LAGEOS, GRACE, GOCE) and surface data (DTU 2'x2' global gravity anomaly grid).We used the free-air gravity anomaly, downloaded with the resolution of ETOPO1 (1 Arc-Minute), from the International Centre for Global Earth Models-ICGEM (Barthelmes et al., 2016).The free-air anomaly map of the study area (Fig. 4a) displays generally low gravity values (±20 mGal) over the basin area indicating that the basin is largely isostatically compensated.An exception is a pronounced negative anomaly with values as low as -80 mGal in the northwestern segment of the Marmara Sea around the NAFZ.Comparing the bathymetry (Fig. 1c) with the freeair gravity anomaly map, it is evident that this negative anomaly is not related to a larger basin depth as bathymetry is rather uniform along the entire axial part of the basin.Likewise, the basement of the syn-kinematic sediments (Fig. 2a) is in the same range in both sub-basins.Accordingly, the negative anomaly is not due to thickness variations of the young sediments or water depth.Apart from the onshore area next to this negative anomaly, the Marmara Sea basin is surrounded by a chain of positive free-air gravity anomalies in a range of +70 to +120 mGal that largely correlate with high topographic elevations.In summary, it can be stated that apart from the local negative anomaly domain, the syn-kinematic sediments need to be isostatically balanced in the crust, given that the Moho topography is varying on a far longer wavelength below the basin.
The seismic observations considered for this study, in addition to those taken into account in the initial model, include P-wave velocity profiles from an offshore-onshore reflection-refraction survey (Bécel et al., 2009) and from a 3D seismic tomography study focused on the sediment-basement configuration of the North Marmara trough (Bayrakci et al., 2013).Both studies are based on the SEISMARMARA-Leg1 seismic survey (Hirn and Singh, 2001), and the location of the related profiles in the model area are shown in Fig. 4b.Based on 3D seismic tomography modelling in the North Marmara trough (Bayrakci et al., 2013) indicates that the P-wave velocities vary between 1.8 and 4.2 km.s -1 within the syn-kinematic sediments, and derives the top of the crystalline basement as an iso-velocity surface with a P-wave velocity of 5.2 km.s -1 .In addition, relying on wideangle reflection-refraction modelling, Bécel et al. (2009) interpreted a refractor below the base syn-kinematic sediments with a P-wave velocity close to 5.7 km.s -1 as the top of the crystalline basement.Furthermore, they interpreted a reflective horizon with a P-wave velocity of 6.7 km.s -1 and largely parallel to the Moho topography as the top lower crystalline crust.Moreover, multichannel seismic reflection data collected in the southwestern part of the Central Basin and in the northeastern part of Marmara Island, documented a 43 km long low-angle dipping reflector interpreted as a normal detachment fault cutting through the upper crystalline crust down to the lower crust (Bécel et al., 2009).In brief, these seismic studies suggest that the crust beneath the syn-kinematic sediments is not homogeneous as assumed in the initial model, but that there is a unit of pre-kinematic sediments beneath the syn-kinematic sediments with an average Pwave velocity of 4.7 km.s -1 above the crystalline crust.Within the upper crystalline crust, the P-wave velocity varies between 5.7 km.s -1 at the top of the crystalline basement to 6.3 km.s -1 above the top of the lower crystalline crust.Lateral velocity variations (~ 0.3 km.s -1 ) are also observed surrounding the detachment fault in the upper crystalline crust.

Results
In addition to the initial structural model with a homogeneous crustal layer below the syn-kinematic sediments (Fig. 3), relying on seismic profiles (Fig. 4b), we modified the structural model that differentiate crustal layers (Fig. 5).The gravity response of these two different 3D structural density models and their corresponding residual gravity anomaly are shown in Fig. 6 and   7, respectively.Here, we show these models and present a best-fit model obtained from the 3D forward gravity modelling (Fig. 8, 9, 10).These results are summarized in Table 1.

Initial model
The initial model (Hergert and Heidbach, 2010;Hergert et al., 2011) resolves only the three structural units: water, synkinematic sediments, and a homogeneous crust (Fig. 3).Hergert et al. (2011) considered a depth-dependent density gradient based on seismic velocities for the sediments and crust.The gradient profile varies between 1700 to 2300 kg.m -3 within the syn-kinematic sediments, between 2500 to 2700 kg.m -3 for the first 20 % of the crust, and from 2700 to 3000 kg.m -3 for the for the syn-kinematic sediments and the crust, respectively.
The calculated gravity response of the initial model (Fig. 6a) indicates a significant misfit with respect to the observed gravity (Fig. 4a).In the onshore domain and the eastern part of the model, the misfit between observed and modelled gravity is rather small and ranges between ±20 mGal (Fig. 7a).Furthermore, within the offshore domain, along the MMF, there are two local positive residual gravity anomalies with more than +90 mGal (A and B in Fig. 7a).These positive anomalies indicate mass deficits in the model and spatially correlate with the bends along the MMF: one occurs in the southern part of the Princes Islands Segment, between the Tuzla Bend and the Istanbul Bend, and the other one is present south of the Ganos Bend.There is also a local short-wavelength positive residual anomaly, reaching values higher than +60 mGal at the location of the Imralı Basin (C in Fig. 7a).In addition, a pronounced West-East oriented continuous negative residual anomaly of around -50 mGal is detected adjacent to the southern coastline.Overall, these results indicate that the long-wavelength gravity field is reproduced by the initial model and that the Moho topography (Fig. 2b) is consistent with observed gravity.However, the large residual anomalies of a few tens of km in diameter indicate the presence of crustal density heterogeneities causing gravity anomalies of smaller wavelengths, i.e. shallower depth.

Differentiated crust
In addition to this indication of density heterogeneities in the crust from gravity, also seismic observations (e.g.Laigle et al., 2008;Bécel et al., 2009;2010;Bayrakci, 2009;Bayrakci et al., 2013) point to crustal heterogeneity expressed as distinct lateral and vertical variations in seismic velocity.To integrate the outcomes of the seismic studies, we vertically differentiate the crust in a next step into three units: (1) a unit of pre-kinematic sediments, (2) a unit of upper crystalline crust, and (3) a lower crystalline crustal unit.

Pre-kinematic sediments
In the initial model (Hergert and Heidbach, 2010;Hergert et al., 2011), the upper limit of the crust below the syn-kinematic sediments (their "top-basement") was mainly defined as pre-kinematic Cretaceous limestone (Ergün and Özel, 1995;Parke et al., 2002;Le Pichon et al., 2014): a surface corresponding to an increase of P-wave velocity to values larger than 4.5 km.s -1 .Furthermore, Bécel et al. (2009) interpreted a top crystalline basement as a surface where P-wave velocity increases to values above 5.7 km.s -1 based on seismic imaging.In addition, Bayrakci et al. (2013) derived the top of the crystalline crust at an isovelocity surface of 5.2 km.s -1 based on a 3D P-wave tomography model beneath the North Marmara Trough.These seismic observations justify the differentiation of an additional unit of pre-kinematic sediments.Accordingly, we implement a unit the upper limit of which corresponds to the top of the pre-kinematic Cretaceous limestone (=base syn-kinematic sediments in the initial model) while its base corresponds to the top crystalline basement (Fig. 5).The top crystalline crust topography proposed by Bécel et al. (2009) and by Bayrakci et al. (2013) is similar, and the depth difference between the surfaces presented in the two studies is mainly less than 2 km.Therefore, we derive the geometry of the top crystalline basement for the gravity test applying a convergent interpolation between the seismic profiles (Fig. 4b) of Bayrakci et al. (2013) and of Bécel et al. (2009).
As the newly implemented pre-kinematic sedimentary unit represents Pre-Marmara Sea deposits, it is mostly absent in the realm of the present-day Marmara Sea (Fig. 5a).Its thickness displays maxima of up to 6.5 km along the north-western and southern margins of the present-day Marmara Sea and significantly decreases eastwards to less than 1.5 km.Bayrakci et al. (2013) showed that the average velocity of the pre-kinematic sediments is around 4.7 km.s -1 .To convert the velocity information for this unit into density, we use an empirical equation (Eq. 1) which is a polynomial regression to the Nafe-Drake Curve valid for P-wave velocities between 1.5 to 8.5 km.s -1 (Brocher, 2005after Ludwig et al., 1970).

Crystalline crust
Apart from the unit of pre-kinematic sediments, the P-wave velocity model of Bécel et al. (2009) differentiates an additional crustal interface across which P-wave velocities increase from values of around 6.2 km.s -1 above the interface to values higher than 6.7 km.s -1 below.They interpreted this interface as the top of the lower crystalline crust.Consequently, we applied a convergent interpolation between the seismic profiles (Fig. 4b) of Bécel et al. (2009) to derive the top lower crystalline crust implemented into the next model.Eventually, we considered the thickness between the top crystalline basement and the top of the lower crystalline crust as the upper crystalline crustal unit.Its thickness distribution (Fig. 5b) shows pronounced thickness minima below the thickness maxima of the syn-kinematic sediments, where the upper crystalline crust is less than 12 km thick.
In contrast, the upper crystalline crust is up to 23 km thick below the south-western margin of the present-day Marmara Sea and reaches more than 26 km in thickness along the eastern margin.
Below the upper crystalline crust, a lower crystalline crustal unit follows, bounded on its top by the top lower crystalline crust and at its base by the Moho discontinuity.It is characterized by an almost uniform thickness distribution (Fig. 5c) of around 10 km across the model area.In the north-western corner of the model area, where the Moho surface (Fig. 2b) descends, the thickness of the lower crystalline crust reaches its maximum of up to 15 km.In contrast, the lower crystalline crust thins to less than 6 km below the south-western margin of the present-day Marmara Sea, where the upper crystalline crust thickens to 23 km.Two other local areas where the lower crystalline crust has a reduced thickness (less than 7 km) are aligned in a N-S direction in the eastern central part of the study area: one adjacent to the Bosporus in the North, and a second one south of Armutlu Peninsula.
Throughout the upper crystalline crustal unit, seismic velocities increase with depth from 5.7 km.s -1 to 6.3 km.s -1 (Bécel et al., 2009).Therefore, we considered 6 km.s -1 as the average P-wave velocity of the upper crystalline crust.P-wave velocities for the lower crystalline crust show less variation, thus, 6.7 km.s -1 has been adopted as the average P-wave velocity within the lower crystalline crust.The density for both crystalline crustal layers, are calculated respecting the P-wave velocities (Eq. 1) as 2720 kg.m -3 and 2890 kg.m -3 for the upper and lower crystalline crust, respectively.
The gravity calculated for this refined model shows a better fit with the observed gravity field in comparison to the initial model (Fig. 7a and 7b).Nevertheless, the three local large positive residual gravity anomalies observed for the initial model (A, B, and C in Fig. 7a) are still evidently indicating that the implemented subdivision of the crust alone is insufficient.The short-wavelength positive anomaly at location C could be interpreted as a local lack of mass within the modelled sedimentary fill of the Imralı Basin.However, the wavelength of the two other positive residual anomalies at A and B is too large to be caused by a high-density feature at the sedimentary fill level but too small to be a result of density heterogeneities in the mantle.
Thus, we concluded that these misfits are most likely related to high-density bodies within the crystalline crust.

Best-fit Model
To overcome these remaining misfits between modelled and observed gravity (A and B in Fig. 7a), we applied forward gravity modelling to test different structural density configurations within the crust.The gravity response of the best fit model is shown in Fig. 6c and the residual in Fig. 7c.Over most of the model area, the residual gravity anomaly (Fig. 7c) shows differences between modelled and observed gravity of ±20 mGal.Achieving this fit required the implementation of two dome-shaped high-density bodies of considerable dimension in the crystalline crust.In the best-fit model, these bodies have an average density of 3050 kg.m -3 , being thus denser than the lower crystalline crust (average density 2890 kg.m -3 ), but less dense than the Mantle (3300 kg.m -3 ).They extend from the Moho upward, cutting through the lower crystalline crust, and reaching into the upper crystalline crust as shallow as the top crystalline basement at about 4 km depth.Accordingly, the high-density bodies attain thicknesses of up to 26 km (Fig. 8 and 9).
The position of these high-density bodies spatially correlates with the domains where the MMF bends (Fig. 9 and 10).At the western margin of the Marmara Sea and below the Ganos Bend, the high-density body cuts the lower crystalline crust at a depth of around 21 km b.s.l and continues through the upper crystalline crust.The shallower part of this body (less than 5 km b.s.l) is located directly south of the Ganos Bend, where the MMF changes its strike direction from E-W to ENE-WSW (Fig. 10).Likewise, the second high-density body is modelled beneath the Princes Islands Segment at the eastern margin of the Marmara Sea, and the top of the body is located at a depth of less than 5 km b.s.l (Fig. 10).
By introducing the two high-density bodies into the structural model, eventually, the thickness distribution of the upper and lower crystalline crust has changed (Fig. 9) below the Çınarcık and Tekirdağ basins, where the high-density bodies largely replace the crystalline crustal units.Over the rest of the model area, the thickness distribution of the crystalline crustal units is similar to the one in the model explained in Sect.4.2.2.Remarkably, the long axis of the high-density bodies follows the strike directions of the bent segments of the MMF (Fig. 9 and 10).In addition, a spatial correlation is evident between the location of the two high-density bodies with the position of the young depo-centres of the Çınarcık and Tekirdağ basins as indicated by deepest present-day bathymetry and by thickness maxima of the syn-kinematic sediments (Fig. 1c and 3).

Interpretation and discussion of the best-fit model
The response of the best-fit gravity model (Fig. 6c) and its corresponding misfit (Fig. 7c) confirmed that the crust below the Marmara Sea is characterized by significant density heterogeneities.In summary, the model predicting gravity best resolves six structural units with different densities that indicates different lithological settings (Table 1).
The uppermost and youngest layer is the present-day water column (Fig. 3a) that is largest in the present-day sub-basins of the Marmara Sea and underlain by the unit of syn-kinematic sediments of the Marmara Sea (Fig. 3b).These syn-kinematic The third modelled unit is characterized by an average density (2490 kg.m -3 ) and by observed seismic velocities (4200-5200 m.s -1 ) representative for sediments.At the same time, the unit is largely missing below the present-day Marmara Sea.We, therefore, interpret this unit as a Pre-Marmara Sea sedimentary unit above the top crystalline basement.The areas where the maximum thickness of more than 6 km are modelled for the pre-kinematic sediments (NW and S of the Sea of Marmara) coincide spatially with the location where Pre-Neogene rocks are present according to surface geology (Yaltırak, 2002).Other surface geological observations (Ergün and Özel, 1995;Genç, 1998;Turgut and Eseller, 2000;Le Pichon et al., 2014) also report the presence of Eocene-Oligocene sediments at the location where the maximum thickness of the pre-kinematic sediments unit is modelled.
The sedimentary units are underlain by the upper crystalline crust, which is thinned below both the Marmara Sea and the prekinematic sediments.This indicates that upper crustal thinning accompanied both phases of basin evolution.Both, the modelled average density and observed seismic velocities for the upper crystalline crust indicate that this unit is dominantly composed of felsic metamorphic rocks.A comparison of the average density (2720 kg.m -3 ) and average P-wave velocity (6 km.s -1 ) of the upper crystalline crust with a compilation velocity-density pairs derived from laboratory measurements (Christensen and Mooney, 1995) indicates a composition corresponding to phyllites and/or biotite gneisses.
Below the upper crystalline crust, the lower crystalline crust follows, the top of which is largely parallel to the Moho topography.The thickness of this unit (Fig. 5c 6.7 km.s -1 ) and the property compilations of Christensen and Mooney (1995), the lithology of the lower crustal unit could be interpreted as diorite and/or granulite.
The sixth unit resolved in the best-fit gravity model encompasses two high-density bodies rising from the Moho in a domeshaped manner through both crystalline crustal layers, partly up to the top crystalline basement (Fig. 8 and 9).For these bodies, a rather high density (3050 kg.m -3 ) has to be assumed which indicates that they are of mafic composition.Considering the seismic velocity and density relationship (Eq.1), a corresponding average P-wave velocity for such a high-density body would be around 7.27-7.28km.s -1 .This combination of physical properties would indicate gabbroic intrusive rock (Christensen and Mooney, 1995) as a possible lithological interpretation for the high-density bodies.Their locations correlate with the bent segments of the MMF (Fig. 9 and 10) indicating that such a mafic composition in concert with their considerable thickness could result in greater strength compared to the surrounding felsic upper crust or the intermediate-mafic lower crust.
The modelled upper and lower crystalline crustal units are in consistency with seismic observations and velocity modelling (Laigle et al., 2008;Bécel et al., 2009;2010;Bayrakci et al., 2013).In contrast, seismic studies did not report the presence of large high-velocity bodies that would coincide spatially with the modelled high-density bodies.There are only a few indications from seismic tomography (Bayrakci et al., 2013) discriminating a zone of high P-wave velocity (Vp > 6.5 km.s -1 ) beneath the top crystalline basement beneath the Çınarcık Basin.This high-velocity zone approximately correlates with the top of the highdensity body in this area.In addition, other tomography results (Yamamoto et al., 2017) indicate a zone of higher S-wave velocity and slightly higher P-wave velocity at about 20 km depth b.s.l, in the area where the western high-density body cuts the boundary between the upper and the lower crystalline crust.

Model limitations
The gravity response of the best-fit model presents a good fit (±20 mGal) over most of the model area.Nevertheless, there are still some negative residual gravity anomalies where observed and modelled gravity differ by more than 70 mGal at Marmara Island (D in Fig. 7c) and by more than 50 mGal at the Armutlu Peninsula (E in Fig. 7c).These short-wavelength negative residual anomalies indicate that shallow low-density features remain unresolved in the model.
The thickness distribution maps (Fig. 3, 5, and 9) show that Marmara Island is dominantly exposing rocks of the upper crystalline crust.More precisely, geological surface observations in this area (Aksoy 1995;1996;Attanasio et al., 2008;Karacık et al., 2008;Ustaömer et al., 2009) differentiate three main rock types in outcrops: A Permian Marble unit in the North, an Eocene Granodiorite unit in the centre, and a Permian Metabasite in the South of Marmara Island.Considering the residual anomalies (Fig. 7c), these three units have densities that are different from the average density assumed for the upper crystalline crust (2720 kg.m -3 ).This could explain the obtained misfit with the observed gravity in this region (D in Fig. 7c).
Our result of obtaining a negative residual indicates that the subsurface extent of rocks with densities lower than the assumed average for the upper crystalline crust is larger than that of the units with higher densities.In other words, the marbles would make a larger portion of the island's subsurface than the metabasites or granodiorites.The negative residual anomaly at Armutlu Peninsula (E in Fig. 7c) is found where the syn-kinematic sedimentary unit is absent (Fig. 3b), whereas a thickening of the pre-kinematic sediments is modelled there (Fig. 5a).Geological maps (Genç, 1998;Yaltırak, 2002;Akbayram et al., 2016) show that this area is mainly covered by Pre-Neogene basement, Miocene acidintermediate volcanic rocks, and some Pliocene-Holocene clastic sediments.However, the model does not account for these locally documented occurrence of syn-kinematic sediments (Pliocene-Holocene clastics) and of Miocene volcanic rocks in this domain, which overall could explain the negative residual anomaly.

Implications
The gravity modelling demonstrates that considering a homogenous crystalline crust beneath the Sea of Marmara is not a valid assumption, but that rather a two-layered crystalline crust crosscut by two large local high-density bodies is plausible.The mechanisms and timing of the emplacement of the high-density bodies are, however, difficult to determine.The modelled density indicates that the high-density bodies represent magmatic additions to the Marmara crust, potentially originating from larger depths and have buoyantly risen into domains of local extension.The spatial correlation between the position of the high-density bodies and the position of thickness maxima in the syn-kinematic sediments indicates that subsidence in the synkinematic basins at least partly took place in response to cooling of the previously emplaced (magmatic) high-density bodies.
This would imply that the emplacement of the high-density bodies predates the formation of the Marmara Sea sub-basins and the propagation of the MMF.To assess the possible contribution of thermal cooling to the subsidence history of the Marmara Sea, a detailed subsidence analysis with determination of the tectonic subsidence would be required.
A second interesting finding is the spatial correlation between the position of the high-density bodies and the bent segments of the MMF.If the high-density bodies represent high-strength domains of the Marmara Sea crust, it would cause local stress deviations influencing the fault propagation direction.The 3D view on the fault plane in relation to the high-density bodies position illustrate how the fault bends in these high-strength domains (Fig. 10).This would imply that the emplacement of the high-density bodies also predates the propagation of the MMF.Such an interpretation would support the previously proposed hypothesis that the NAFZ reached the eastern part of the present-day Marmara Sea (Izmit) around 4 Ma before present, when the area was a domain of distributed (trans)tensional deformation, and started to propagate beneath the present-day Sea of Marmara as the MMF about 2.5 Ma ago (Le Pichon et al., 2014;2015).
Another implication from density modelling is that the compositional and therefore also rheological heterogeneity of the Marmara crust may result in a differential response of the area to present-day far-field stresses.Accordingly, conclusions drawn from earlier studies investigating the stress-strain state in the region of the Marmara Sea (Hergert andHeidbach, 2010, 2011;Hergert et al., 2011) need to be checked if there are still valid.
One of the important discussions in the area of the Marmara region is on aspects that govern the dynamics of the MMF, where a 250-year lasting seismic gap in the southern vicinity of Istanbul is observed.The western segment of the MMF is considered as a partially creeping segment (Schmittbuhl et al., 2016;Bohnhoff et al., 2017b), whereas the eastern-central segment of the MMF is thought to be locked down to 10 km depth (Bohnhoff et al., 2013;2017b;Ergintav et al., 2014;Sakic et al., 2016).hypothesis could have implications for hazard and risk assessment in this area, but need to be tested by geodynamic models considering thermo-mechanical principles.

Conclusions
In this study, a 3D density model for the Sea of Marmara was presented that integrates available seismological observations and is consistent with observed gravity.Testing successively models of increasing complexity, a best-fit model is derived that resolves six structural units with different densities (Table 1).From our results we conclude: (1) The present-day seafloor of the Marmara Sea has a more complex structure than during the phase of its initiation and is structured into the three main depocentres of the Tekirdağ Basin, the Central Basin, and the Çınarcık Basin.This structural pattern is consistent with a pull-apart setting in a releasing zone of the NAFZ.
(2) Below the present-day seafloor, the unit of syn-kinematic sediments of the Marmara Sea indicates that two main depocentres were subsiding during the early phase of basin formation.A lower sedimentary unit is interpreted as pre-kinematic sediments of the Marmara Sea.The sedimentary units are underlain by the felsic metamorphic upper crystalline crust that is significantly thinned below the basin.The lowest crustal layer of regional extent is the intermediate to mafic lower crystalline crust.Both crystalline crustal layers are cut by two up-doming high-density mafic bodies that rise from the Moho to the base of the syn-kinematic basins.The Moho surface used in this study is consistent with the satellite gravity observation.
(3) The emplacement of the high-density bodies within the crystalline crust could have a causal relationship with the basinforming mechanism.
(4) The spatial correlation between the high-density bodies with the bent segments along the MMF indicates that rheological contrasts in the crust may control the propagation and movement of the MMF; these high-density bodies are a possible explanation for the bends of the MMF, and support the hypothesis that the MMF is geomechanically segmented.
(5) The high-density bodies may have an impact on the stress variability along the MMF.Therefore, an update of the existing geomechanical model is now essential to account for this more detailed density distribution.In addition, the 3D density model could be considered for further thermo-mechanical modelling of the system dynamics.

Solid
Earth Discuss., https://doi.org/10.5194/se-2018-113Manuscript under review for journal Solid Earth Discussion started: 22 October 2018 c Author(s) 2018.CC BY 4.0 License.levels below the sediments, and (3) a final best fit model in which the remaining residual anomaly is minimized by implementing additional density-geometry changes in the crust but respecting the seismic data.
are present only inside the Marmara Sea domain and their thickness distribution indicates a subsidence regime similar to the present-day one.The relationship between the individual sub-basins of the Marmara Sea and the course of the MMF are however different: The shape of the present-day Tekirdağ Basin is not evident in the thickness distribution of the syn-kinematic sediments, whereas the Central Basin along the MMF and the Çınarcık Basin are largely following their presentday counterparts.This indicates that the differentiation into the present-day Central and Çınarcık basins postdates the synkinematic phase of the Marmara Sea.The average density and the observed seismic velocities indicate that this unit is mainly composed of poorly consolidated clastic deposits.There is, however, little information on their precise ages; suggested time intervals for the deposition of this unit range from Late Miocene to Holocene with a longer deposition portion of the unit assigned to Pliocene to Holocene times(Le Pichon et al., 2014;2015).
) indicates no clear relationship with the formation of both generations of pre and syn-kinematic basins.Here, the modelled average density and observed seismic velocities are indicative for an intermediate to mafic composition.Combining the physical properties of the lower crystalline crust (ρ = 2890 kg.m -3 & V p = Solid Earth Discuss., https://doi.org/10.5194/se-2018-113Manuscript under review for journal Solid Earth Discussion started: 22 October 2018 c Author(s) 2018.CC BY 4.0 License.

Solid
Earth Discuss., https://doi.org/10.5194/se-2018-113Manuscript under review for journal Solid Earth Discussion started: 22 October 2018 c Author(s) 2018.CC BY 4.0 License.The reasons why this seismic gap of the MMF has not ruptured over the past 250 years are debated.The felsic to intermediate crustal composition deduced from our gravity model would favour creep between the two crustal high-density bodies, whereas the two domains of the high-density bodies could represent locked segments that would require high-stress levels to fail.In case of failure, however, the energy would probably be released in a strong earthquake.These high-density bodies are interpreted as mafic and therefore represent stronger material than the surrounding felsic to intermediate crustal material of the same depth.Such rheological heterogeneities would explain the distribution of different deformation modes with creeping segments in the felsic to intermediate crustal domains and locked to critically stressed segments in the mafic domains.This

Figure 1 :
Figure 1: The location and the tectonic setting of the model area; (a) Plate tectonics map of the Anatolian plate and its relation to the Arabian, African, and Eurasian Plates; (b) The NAFZ propagation and its branches in NW Anatolian plate; (c) The model area (WGS84 UTM Zone 35N) including the relief map of the Sea of Marmara and its surrounding onshore domain, the seismic gap since 1766 (thick red line; Bohnhoff et.al., 2017b), and the faults system identified in the investigations by Armijo et al. (2002; 2005), and 5

Figure 2 :
Figure 2: Main horizons within the initial model (WGS84 UTM Zone 35N); (a) Depth to top-basement; (b) Depth to Moho.The corresponding thickness maps are illustrated in Fig. 3.

Figure 6 :
Figure 6: Calculated gravity over the model area (WGS84 UTM Zone 35N): (a) Initial model gravity response; (b) Gravity response of a model with differentiated crust based on the seismic observations (Fig. 5); (c) Gravity response of the best-fit model based on the forward gravity modelling.The corresponding residual gravity anomaly of each model is shown in Fig. 7.

Figure 7 :
Figure 7: Residual gravity anomaly maps show the misfit between the observed (Fig. 4a) and calculated gravity (Fig. 6) of different structural model across the study area (WGS84 UTM Zone 35N): (a) Initial model; (b) Model with a differentiated crustal unit; (c) Best-fit model based on the 3D forward gravity modelling.AA´, BB´, and CC´ show the location of the cross-sections in Fig. 8.

Figure 8 :
Figure 8: Cross-sections of the observed and calculated gravity based on the best-fit gravity model and the corresponding structural settings.The profile locations are shown in Fig. 7c: (a) AA´; (b) BB´; (c) CC´.

Figure 9 :
Figure 9: Thickness map of different units of the crystalline crust based on the best-fit gravity model (WGS84 UTM Zone 35N): (a) Upper crystalline crust; (b) Lower crystalline crust; (c) High-density bodies.

Figure 10 :
Figure 10: 3D view of the Moho, the high-density bodies, and the MMF plane across the model area (WGS84 UTM Zone 35N).The high-density bodies location spatially correlates with the bent segments of the MMF.The Moho depth and the 3D fault plane from Hergert and Heidbach (2010).