Interpretation of gravity anomalies to delineate some structural features of Biga and Gelibolu peninsulas, and their surroundings (north-west Turkey)

Abstract

On the basis of gravity data, derived anomaly traces were presented for the interpretation of some shallow structural features of Biga and Gelibolu peninsulas. Since building general understanding of subtle details about subsurface geology is of great importance considering that the study area is tectonically important, some advanced data processing techniques were implemented to gravity anomalies in a detailed manner. The procedures were performed using a MATLAB-based software package (Gravity and Magnetic Interpretation – GMINTERP). First, a finite element method was utilised to produce the residual data-set which is expected to reflect short wavelength anomalies arising from shallower geological structures, and thereafter some derivative-based algorithms were executed to analyse the residual data. The general anomaly patterns obtained from the applications clearly corresponded to the well-known surface geology map of the study area. Derivative-based anomaly maps put forward some findings about the existence of an old caldera structure in the western part of the Biga Peninsula. Additionally, abrupt lateral changes in anomaly amplitudes indicated the presence of some major structural discontinuities. Thus, findings yielded to make significant geological interpretations that might be important for further investigations. This study also showed that GMINTERP software package proved useful in assisting geological interpretation using geophysical potential field data-sets.

Keywords:

• gravity anomalies
• structural features
• data processing
• GMINTERP
• Biga and Gelibolu peninsulas

1. Introduction

Herein, we aim to contribute to the understanding of the nature of some shallow structural features of the Biga and Gelibolu peninsulas, and their surroundings by analysing Bouguer gravity anomalies. It is well known that the complex geological structure of the western Anatolia has attracted intense attention so far, and numerous geophysical investigations have been performed using potential field data-sets. Some of these studies have involved entire or some parts of the area that is presented here, and they have mainly focused on mapping the geologic contacts, faults and the seismogenic zones, and determining the crustal features and crustal thickness of the region (e.g. Albora, Uçan & Aydoğan, 2007; Ateş, Bilim, Büyüksaraç, & Bektaş, 2008; Ateş, Kayıran, & Sincer, 2003; Aydoğan, 2007; Aydoğan, Elmas, Albora, & Uçan, 2005; Hisarlı & Orbay, 2000; Oruç & Keskinsezer, 2008a; Ulugergerli, Seyitoğlu, Başokur, Kaya, & Dikmen, 2007). However, except for some recently reported works (Ekinci & Yiğitbaş, 2012a, 2012b) which have studied the airborne magnetic anomalies of the igneous rocks that occupy a considerable amount of area in the Biga Peninsula and its surrounding, no special geophysical investigation has been yet performed in the entirety of the Biga and Gelibolu peninsulas and their surroundings which focuses on the shallow geological structure of the region. Hence, to fill that gap, we attempted to investigate the shallow subsurface geology and structural features in some detailed framework by interpreting and correlating the gravity anomalies with the geology of the region which is better known now by the recent regional geological projects of Duru, Dönmez, Ilgar, Pehlivan, and Akçay (2012) and Yiğitbaş, Şengün, and Tunç (2014).

Gravity anomaly interpretation procedures consist of many special techniques depending on the quality of the data-set and the goals of the analysis (Hinze, von Frese, & Saad, 2013). These data processing techniques may facilitate the interpretation of the potential field anomaly maps, and may also aid geological mapping (Blakely, 1995). Gravity anomaly maps are generally analysed by employing some linear transformations, directional derivative-based techniques and inverse modelling techniques. Among them, linear transformations may assist geological interpretation, and provide some insights about the nature of the causative geological sources (Blakely, 1995). Analytic continuations (upward and downward) are the most commonly used linear transformations in gravity data processing (e.g. Ateş et al., 2012; Aydemir, Ateş, Bilim, Büyüksaraç, & Bektaş, 2014; Büyüksaraç, 2007; Gibert & Galdeano, 1985; Gunn, 1975). Directional derivative-based techniques commonly use vertical or horizontal gradients (derivatives) of degree first or second, analytic signal (AS) amplitudes or their various combinations in order to improve the edges, and also to outline the possible boundaries of the gravity and pseudogravity sources (e.g. Balkaya, Göktürkler, Erhan, & Levent, Ekinci, 2012; Blakely & Simpson, 1986; Debeglia & Corpel, 1997; Oruç, 2011, 2014; Roest, Verhoef, & Pilkington, 1992). Inversion techniques are frequently used to estimate the model parameters whose responses are similar to the measured data, but the well-known under-determined and non-uniqueness nature of the potential field data makes the processing and interpretation difficult. Thus, the inverse modelling problem intensely requires prior information and some constraints in order to recover interpretable and realistic model solutions (e.g. Camacho, Montesinos, & Vieira, 2000; Ekinci, 2008a, 2008b; Feng, Meng, Chen, & Zhang, 2014; Last & Kubik, 1983; Li & Oldenburg, 1998; Roy, Sen, McIntosh, Stoffa, & Nakamura, 2005; Zhdanov, 2009; Zhdanov, Ellis, & Mukherjee, 2004; Sun, Liu, Shou & Zhang (2012)).

The success of the interpretation of potential field data-sets undoubtedly depends on the selection of the data processing techniques that are proper to the nature of the problem. Consequently, potential field data-sets generally need many processing techniques in order to produce meaningful results which are in accordance with the nature of the geology of the area of interest (Ekinci & Yiğitbaş, 2012b, 2012c). Numerous commercial software packages (e.g. Geosoft Oasis Montaj, Geosystem WinGlink, MagPick and IGMAS) are commonly used for analysing the potential field data-sets by employing some of the techniques specified above. Additionally, many freeware and open-source codes which have focused on particular issues and problems can be found in the literature (e.g. Arısoy & Dikmen, 2011; Cooper, 2009; Cooper & Cowan, 2005, 2006, 2008; Durrheim & Cooper, 1998). However, freeware and commercial software packages may not meet the requirements in some cases, since they do not include all the required data processing techniques in their own libraries. Therefore, an efficient toolkit named as Gravity and Magnetic Interpretation (GMINTERP) that performs numerous analysing techniques was developed and used for processing the gravity data-set of the Biga and Gelibolu peninsulas, and their surroundings. A brief overview of the developed toolkit was previously given by Ekinci and Yiğitbaş (2012b). Before the real gravity data analyses, a large number of synthetic simulations were implemented with success, indicating the reliability of the GMINTERP software. The MATLAB-based interactive toolkit which consist of a set of linked functions in conjunction with a graphical user interface (Ekinci, 2010a; Ekinci & Yiğitbaş, 2012b, 2012c) produced some noteworthy anomaly maps from the gravity data-set of the study area and we discussed the geological implications of those derived gravity signatures.

2. Study area

Turkey lies within a highly active tectonic belt which is characterised by significant neotectonic structures (Figure 1), and consists of several continental fragments which were joined together into a single landmass finally in the late Tertiary (Okay, 2008). These continental terranes were separated by oceans whose relics are widely distributed throughout the Anatolia during most of the Phanerozoic (Okay, 2008). Strandja, İstanbul and Sakarya are the well-known terranes in northern Turkey (Figure 2), and they are called as the Pontides (Okay, 2008). Armutlu–Ovacık Zone (Figure 2) is a tectonic mix of İstanbul and Sakarya zones (Elmas & Yiğitbaş, 2001, 2005; Yiğitbaş, Elmas, & Yılmaz, 1999). The Pontides, the Anatolide–Tauride block and the northern part of Arabian Plate are pronounced as the main tectonic units of Turkey (Ketin, 1966), and they are separated by sutures that represent the boundaries of former lithospheric plates i.e. tectonic lines (Okay & Tuysuz, 1999). The İzmir–Ankara–Erzincan suture (Figure 2) that was formed by the closure of the northern branch of the Neo-Tethys Ocean separates the Pontides and the Anatolide–Tauride block from each other (Şengör & Yılmaz, 1981). The northern part of the Arabian Plate was separated from the Anatolide–Tauride block by the southern branch of the Neo-Tethys Ocean (Şengör & Yılmaz, 1981). The continuing convergence of the African, Arabian and Eurasian plates since the Middle Miocene, and the westward migration of the Anatolian Plate predominantly effects the neotectonic regime of Turkey (Figure 1) (Dewey & Şengör, 1979; Jackson & McKenzie, 1984, 1988; Ketin, 1948; McKenzie, 1970, 1972; Şengör & Yılmaz, 1981).

Figure 1. Simplified neotectonic and relief maps of Turkey and the surrounding area (compiled with Yiğitbas, Elmas, Sefunç, and Özer, 2004) and CGIAR-CSI GeoPortal, 2012). Abbreviations: SBT, Southern Black Sea Thrust; NAFZ, North Anatolian Fault Zone; NEAFZ, Northeast Anatolian Fault Zone; WAEP, West Anatolia Extensional Province; EAFZ, East Anatolian Fault Zone; DSFZ, Dead Sea Fault Zone.

Figure 2. Simplified tectonic map of Turkey and surrounding regions, showing the main sutures and continental fragments (modified after Tunç et al., 2012) and the location of the study area (see rectangle)

The study area, Biga and Gelibolu peninsulas and their surroundings, is located in north-west Turkey, and covers an area of about 21600 km² (black rectangle in Figure 2). It is bounded in the west by the Aegean Sea, in the north by Thracian Basin and Strandja Massif along the Intra-Pontide suture (Şengör & Yılmaz, 1981), and in the south by the Edremit Gulf. Tectonically, the study area is situated at the western end of the Sakarya Zone (Okay, Siyako, & Bürkan, 1991; Şengör & Yılmaz, 1981) which is separated by the İzmir–Ankara–Erzincan suture zone from the Anatolide–Tauride block (Şengör & Yılmaz, 1981).

2.1. Geology of the Gelibolu Peninsula

The upper Cretaceous Lört formation is the oldest unit in the Gelibolu Peninsula (Figure 3). This formation composed mainly of pelagic limestone is unconformably overlain by the Ypresian–Lower Lutetian (Siyako, Bürkan, & Okay, 1989) regressive sequence (Karaağaç and Fıçıdere formations). These units are overlain by an upper Lutetian–upper Eocene transgressive sequence (Soğucak and Ceylan formations). The late Eocene–late Oligocene era is represented by strata displaying a regressive character (Mezardere and Armutburun formations) (Elmas & Meriç, 1998; Önal, 1986; Sümengen et al., 1987). The basement of the middle–upper Miocene/Quaternary sequence (Gazhanedere and younger formations in Figure 3) in the Gelibolu Peninsula is formed by these entire assembly (Elmas & Meriç, 1998). In the Gelibolu Peninsula, the Anafartalar Fault formed as a positive flower structure and elevated the western parts of the Peninsula (Elmas & Meric, 1998; Yaltırak, 1995). The Ganos Fault bounds the Gelibolu Peninsula along its north-west area (Tüysüz, Barka, & Yiğitbaş, 1998). The Anafartalar and Ganos faults are the main structural features that have induced the evolution of the Çanakkale Strait (Yaltırak, Alpar, Sakınç, & Yüce, 2000). The Anafartalar thrust fault continues throughout the Gazhanedere formation (Albora et al., 2007; Yaltırak, 1995; Yaltırak, Alpar, & Yüce, 1998).

Figure 3. Geological map of the Gelibolu (modified after Dönmez, Akçay, Duru, Ilgar, and Pehlivan 2008) and Ilgar et al., 2008) and Biga peninsulas (modified after Duru et al., 2012).

2.2. Geology of the Biga Peninsula

Biga Peninsula comprises different types of tectonic units that comprises the Sakarya Zone as an Alpine continental fragment in the south and Armutlu–Ovacık Zone which is a late Cretaceous–Eocene shear zone in the north (Elmas & Yiğitbaş, 2001; Şengün, Yigitbaş, & Tunç, 2011; Yiğitbaş, Şengün, & Tunç, 2009). Figure 3 shows the geological map of the Biga Peninsula (Duru et al., 2012). The Biga Peninsula in north-west Turkey is characterised mostly by plutonic and associated volcanic rocks, ophiolites and amphibolite to granulite facies (Şengün, Davis, Tunç, & Yiğitbaş, 2012) metamorphic rocks. Gneiss, amphibolites and marbles form the basement of the Biga Peninsula and high- to medium-grade basement rocks occur in several parts (Okay & Satır, 2000a, 2000b) such as Kazdağ Massif, Çamlıca and Kemer metamorphics (Figure 3). A geological cross-section from Gelibolu to Biga peninsulas (oriented in north–south and north-north-west–south-south-east as shown in Figure 3) showing the major tectonic units is illustrated in Figure 4. In addition to its complex geology, the Peninsula is also situated in a neotectonically active region (Figure 1) that has been affected by three main splay branches of North Anatolian Fault Zone (NAFZ) since the Middle Miocene (Aldanmaz, Pearce, Thirlwall, & Mitchell, 2000; Okay et al., 1991). The southern branch consists of the Bursa, Ulubat, Manyas, Yenice-Gönen, Evciler and Edremit faults (Kürçer, Yalçın, Gülen, & Kalafat, in press), the middle branch continues across the Peninsula through the town of Ezine while the northern branch passes underneath the Sea of Marmara (Aldanmaz et al., 2000). In addition to these branches of NAFZ, the Peninsula has been deformed by the West Anatolian Extensional Province (WAEP) which consists of east–west trending horsts, sediment filled grabens and their basin-bounding active faults (Bozkurt, 2001; Gürer, Sangu, & Özburan, 2006). These deformations has caused dextral strike-slip faults with normal component (Barka, 1992; Barka & Kadinsky-Cade, 1988; Herece, 1990; Şaroğlu, Emre, & Kuşçu, 1992; Siyako et al., 1989), tectonic basins related to these faults and other deformation structures with obvious transcurrent component (Kürçer et al., 2008). Çan-Biga fault zone (ÇBFZ), Sarıköy Fault (SAF) and Yenice-Gönen Fault (YGF) are the active members of NAFZ in the Biga Peninsula (Figure 5). Among them, the YGF is one of the most important active faults, and it has produced a very destructive earthquake Mw = 7.2 on 18 March 1953 (Caputo, 2005) with a focal depth estimation of 10–12 km (Herece, 1985). About 70-km-long dextral surface rupture occurred along the fault, and the fault can be followed through the alluvial area of Gönen and the Neogene deposits near Muratlar village, where the fault scarp is still evident (Kürçer et al., 2008). Active faults and the distributions of the epicentres of the earthquakes in the study area with magnitudes of M ≥ 4 are shown in Figure 5.

Figure 4. A simplified geological cross-section along A–B and B–C lines shown in Figure 3. Note that vertical axis of the cross-section was exaggerated and an approximate topography was used for the illustration.

Figure 5. Topographic map of the study area showing the distributions of the epicentre of the earthquakes with magnitudes of M ≥ 4. The sizes of the stars are proportional to the earthquake magnitudes. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Black lines show coastline. Abbreviations: KF, Kestanbol Fault; EVF, Evciler Fault; ÇBFZ, Çan-Biga fault zone; SF, Sinekli Fault; SAF, Sarıköy Fault; YGF, Yenice-Gönen Fault; PF, Pazarköy Fault; EF, Edincik Fault; HBFZ, Havran–Balıkesir fault zone; OF, Ovacık Fault; BF, Bekten Fault; MFZ, Manyas fault zone; EFZ, Edremit fault zone; ATF, Anafartalar thrust fault; GF, Ganos Fault.

Biga Peninsula is generally dominated by the widespread occurrence of Tertiary plutonic and associated volcanic rocks of variable age and composition which are characteristic of middle–late Tertiary events (Beccaletto, 2003). Ezine–Kestanbol, Eybek, Ilıca–Şamlı, Karabiga and Evciler are the high-level granitoid intrusions which are related to magmatic activity (Demirel, Yıldırım, & Burçak, 2004). The exposed granitoid intrusions trend either north-east, following the major tectonic grain of the Peninsula, or east-north-east, cutting the major tectonic grain (Yigit, 2012). Active faults developed by neotectonic activities and also the young granitic and volcanic rocks have key importance on the geothermal sources in the Biga Peninsula. The Peninsula comprises more than 50 hot springs with temperatures ranging from 30 to 100 °C (Demirel et al., 2004), and therefore it is considered as a significant geothermal region in Turkey. Additionally, high heat flow values clearly dominate the extensional tectonic region in the southern part of the Peninsula (Pfister, 1995).

3. Gravity data and data processing

Gravity surveys of the entire Turkey were carried out by the General Directorate of Mineral Research and Exploration (MTA) of Turkey in the years between 1973 and 1988. MTA has accomplished 62,000 gravity observations in Turkey at station intervals of 2–5 km. The Geodesy Department of MTA has provided the locations and elevations of the gravity stations from 1:25,000 scale topographic maps (Ateş, Kearey, & Tufan, 1999). Gravity values were tied to MTA and General Command of Mapping base stations related to the Potsdam 981,260.00 mGal (1 Gal = 1 cm s⁻²) absolute gravity value accepted by the International Union of Geodesy and Geophysics in 1971 (Ateş et al., 1999). The gravity values are related to the International Gravity Standardization Network (Erden, 1979). All the corrections of the observed gravity data were carried out by the MTA. The corrections can be given as follows (Ateş et al., 1999): latitude corrections were carried out using the Gravity Formula of 1967; Bouguer corrections were performed assuming a density of 2.40 Mg m⁻³; topographic corrections were completed to a distance of 167 km assuming a terrain density of 2.40 Mg m⁻³; and lastly free-air and tidal corrections were performed. Complete Bouguer gravity anomaly map of entire Turkey can be seen in some publications (e.g. Ateş et al., 1999; MTA, 2006).

Herein, the gravity anomaly map of the study area was produced with a grid spacing of 2 km, which covers an area of 180 km × 120 km (Figure 6). At the data processing stages, some potential field techniques (i.e. vertical derivative, downward continuation) were performed in wavenumber domain using the fast Fourier transform based on the method given by Blakely (1995). As is well known, the nature of the fast Fourier transformation assumes that a grid is periodic, in effect repeating itself infinitely many times in all horizontal directions (Blakely, 1995). However, if the edges of the data grid do not meet smoothly with their repetitive neighbours the procedure may cause some edge effects (Blakely, 1995). Thus, to remove the possible discontinuities at the grid edges of the data-sets, the size of data grid was augmented to next higher power of 2 by adding artificially produced data bands to the east and north edges of data grid prior to the fast Fourier transform. The extra artificial data bands were removed from the data grid at the end of the operation.

Figure 6. Complete Bouguer gravity anomaly map of the study area. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Abbreviations used for the faults are given in Figure 5.

3.1. Gravity anomaly separation

It is a well-known fact that the potential field anomaly maps present the effects of both shallow and deeper geological structures. These two components are named as regional trends and local (residual) anomalies, and can be expressed as follows (Pawlowski, 1994):(1)

where is the observed gravity field, s and d refer to the gravity response of shallow and deeper structures, respectively. Thus, one of the most important issues in potential field data interpretation is the removal of regional trends when dealing with relatively shallower local geological structures (Dobrin, 1976). Therefore, some mathematical methods are required to separate the map data into two components which are the regional nature and the local fluctuations (Davis, 2002). Since the study presented herein deals with the shallow geological structures of the Biga and Gelibolu peninsulas, regional/residual separation process was applied to gravity data-set in order to estimate the amplitude of the regional background.

Polynomial fitting procedure is the most commonly used technique for potential field anomaly separation. The procedure is based on the minimisation of the sum of the squares of the differences between the data and a polynomial of degree m evaluated at corresponding given observations (Hedgley, 1978). The function to be minimised is expressed as follows (Hedgley, 1978):(2)

where d is the degree of the model, x_i and y_i are the discrete points, a_j is the coefficient to be determined, is the weighting factor and NP is the number of data points. Although there are some statistical approaches such as determining the correlation coefficients or goodness of fit, there have always been some difficulties of how to choose the best degree of the polynomial and what the functional form of the trend should be to obtain the optimum solution without specifying any prior information about nature of the data. When the higher order of the polynomial is preferred the better fit can be achieved and the curve may pass precisely through all the data points, but the regional anomaly may have no reasonable geological meaning (Lowrie, 2007). On the other hand, a finite element method (FEM) based on the use of element shape functions may yield optimum solutions without introducing a prior information (Mallick & Sharma, 1997, 1999). Additionally some other researchers have gained benefit from applying FEM to the potential field data-sets so far (e.g. Agarwal & Srivastava, 2010; Ekinci, 2010b; Ekinci & Yiğitbaş, 2012b; Kaftan, Salk, & Sari, 2005).

In the procedure of FEM, the amplitudes of the regional background are estimated by a weighted sum of discrete gravity values at eight stations coinciding with the eight nodes of a second-order isoparametric element superimposing the map space (Mallick & Sharma, 1999). Firstly, an eight node element is defined in x–y space (real-map space). x_c, y_c is the location of the centre and the sides parallel to the x- and y-axes are defined as 2a and 2b. The transformation of the real-map space to a corresponding reference space is performed by the following substitutions (Mallick & Sharma, 1999):(3)

A reference space is defined by the non-dimensional coordinates ξ and η, ranging between the values of −1 and 1 (Mallick & Sharma, 1999). The regional gravity at any point in the reference element is expressed by a series which is a weighted sum of the nodal gravity values (Mallick & Sharma, 1997):(4)

where the weights, , are the shape functions of the element. The shape functions at the nodes 1 through 8 in reference space are computed by using following formulations (Cheung & Yeo, 1979):(5)

Equation (5(5) ) satisfies that the shape function at the ith node and zero at the other nodes. Further, (Mallick & Sharma, 1999). At the last step, the regional field computed on the reference ξ − η space is transformed to the real-map space by using following coordinate transformation (Mallick & Sharma, 1999):(6)

where are the shape functions and x_i and y_i are the nodal coordinates in the real-map space (Mallick & Sharma, 1997, 1999).

Before the application of regional/residual separation procedure on real gravity data-set, some synthetic simulations were performed to test the success of FEM and the polynomial fitting technique. Considering a rectangular prism is one of the simple ways to approximate a volume of mass. The gravitational attraction of a single rectangular prism is easily computed by integration over the limits of the prism (Blakely, 1995), and the following derived equation can be used for the procedure (Plouff, 1976):(7)

where ρ is the density contrast and γ is the gravitational constant. The other terms are given by:(8)

The gravity response of three vertical-sided prisms was generated for the test simulations by running GMINTERP based on Equation (7(7) ). We used a 0.25 km × 0.25 km grid interval for three-dimensional forward modelling. Model parameters of the causative prismatic bodies are given in Table 1. Figure 7(a) and (b) demonstrates the plan view of three-dimensional prismatic model bodies and the synthetically generated gravity anomaly map, respectively. Note that the gravity response of the model 3 was assumed to represent the regional effect. Figure 8(a)–(f) illustrates the computed regional and residual gravity anomaly maps from the application of FEM and polynomial fitting (second- and third-order) procedures. It is obviously seen that the gravity response of the deeper body (model 3) dramatically enhanced by the use of FEM (Figure 8(a)). However, the maximum gravity amplitudes could not be located over the causative deeper body by the use of polynomial fitting procedure (Figure 8(d) and (f)). It is worth mentioning here that these synthetic simulations indicated the superiority of the FEM over the polynomial fitting procedure. Thus, FEM was performed on the real gravity data-set shown in Figure 6, and regional (Figure 9) and residual gravity anomalies (Figure 10) were computed.

Table 1. Model parameters of three prismatic bodies used for the test simulations.

	Model parameters
Model number	x1 (km)	x2 (km)	y1 (km)	y2 (km)	z1 (km)	z2 (km)	ρ (kg/m^3)
1	12	14	12	14	2	4	100
2	7	10	7	10	5	8	110
3	2	6	2	6	10	15	120

Figure 7. (a) General plan view of three prismatic theoretical model bodies used for the synthetic simulations. (b) Synthetically generated gravity data-set consisting of the anomalies caused by three-dimensional prismatic source bodies shown in (a), Contour interval: 0.4 mGal.

Figure 8. The results of the synthetic simulations obtained by regional/residual separation algorithms. Contour interval: 0.04 mGal. (a) Regional background computed by using FEM. (b) Residual anomaly map computed by using FEM. (c) Regional background computed by using second-order polynomial fitting procedure. (d) Residual anomaly map computed by using second-order polynomial fitting procedure. (e) Regional background computed by using third-order polynomial fitting procedure. (f) Residual anomaly map computed by using third-order polynomial fitting procedure.

Figure 9. Regional gravity background of the study area computed by FEM based on the use of element shape functions. Contour interval: 4 mGal.

Figure 10. Residual gravity anomaly map of the study area. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Black lines show coastline. Abbreviations used for the faults are given in Figure 5.

3.2. Approximating the locations of the causative geological structures

In order to better understand the nature of the residual gravity anomalies and to trace the abrupt lateral changes in densities, firstly, we performed the technique of Blakely and Simpson (1986) to the residual gravity data-set shown in Figure 10. Briefly, the technique is based on a two-step procedure when dealing with gravity data-sets. The shallow source bodies tend to produce gravity or pseudogravity anomalies with maximum horizontal gradients near their edges or boundaries (e.g. Balkaya et al., 2012; Blakely, 1995; Blakely & Simpson, 1986; Ekinci & Yiğitbaş, 2012b), therefore a total horizontal derivative (THD) anomaly map (Figure 11) of the residual gravity data was produced by using the following expression (Blakely, 1995):(9)

Figure 11. THD anomaly map of the data shown in Figure 10. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Black lines show coastline. Abbreviations used for the faults are given in Figure 5.

where Δg is the gravity anomaly. Horizontal derivatives in Equation (9(9) ) were estimated using simple finite-difference method and discrete values of as follows (Blakely, 1995):(10)

where i and j represent the discrete values of on the observation plane at uniform sample intervals Δx and Δy. At the second stage, we attempted to search for the maximum horizontal gradient amplitudes (Blakely, 1995; Blakely & Simpson, 1986) on the grid plane. Computed maxspot map, showing the maxima locations on the grid plane and their amplitudes is demonstrated in Figure 12.

Figure 12. Maxspot map showing the locations of the maximum horizontal gradient amplitudes. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Black lines show coastline. Abbreviations used for the faults are given in Figure 5.

Secondly, a profile curvature anomaly map (Figure 13) was produced to improve the possible geological boundaries and lineaments and also to compare the results with those obtained by analysing the THD magnitudes (Figure 11) of the residual gravity data. Recently, the effectiveness of the technique on gravity and pseudogravity data-sets for contact mapping has been reported by some researchers (Ekinci, Balkaya, Şeren, Kaya, & Lightfoot, 2014; Ekinci, Ertekin, & Yiğitbaş, 2013a, 2013b). Briefly, profile curvature operator measures the rate of change of the slope in the direction of the steepest gradient at each point of the data (Cooper, 2009), and produces a zero contour value at the source edges and either a local maxima or minima over the causative sources (Lee, Morris, Leblanc, & Harris, 2013). The definition of the profile curvature amplitude (P) for a gridded data-set is given by the following equation (Mitášová & Hofierka, 1993):(11)

Figure 13. Profile curvature anomaly map of the data shown in Figure 10. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Black lines show coastline. Abbreviations used for the faults are given in Figure 5.

where g is the gridded data, ∂x and ∂y are the horizontal derivatives while ∂x² and ∂y² are the second-order horizontal derivatives of the gridded data. Here, horizontal derivatives were computed using the terms in Equation (10(10) ), while second-order horizontal derivatives were obtained using the following approximate definitions (Carnahan, Lither, & Wilkes, 1969):(12)

Thirdly, 3D normalized full gradient (NFG) technique was applied to the residual gravity data-set to improve the subtle details of the geological sources, and to follow their downward extensions. NFG technique is based on the downward continuation of the potential field data and the NFG amplitude is computed by dividing the AS amplitude of the field by the average of AS amplitude using the following equation:(13)

where Δz is the downward continuation level and M is the number of data on the grid plane. Downward continuation of a potential field data can be obtained using the following substitution (Blakely, 1995):(14)

where g is the gravity field, F is the Fourier transform, F⁻¹ is the inverse Fourier transform, Δz is the downward continuation level and e is the exponential function. is the radial wavenumber at grid intersections throughout the k_x and k_y plane, and is defined by following equation (Blakely, 1995):(15)

The last component in 3D NFG, the AS amplitude, is a complex quantity and formed through a combination of the vertical and horizontal derivatives of the potential field data (Roest et al., 1992). The maxima locations in AS anomaly map are usually used to locate source bodies and also to improve the edges of them on the grid plane. The AS amplitude is obtained by using the equation as given below (Roest et al., 1992):(16)

where Δg is the measured gravity anomaly and ∂x, ∂y and ∂z are the directional derivatives. In order to produce the AS amplitudes, the horizontal derivatives (∂x and ∂y) were computed in space domain using the simple finite-difference method as expressed previously in Equation (10(10) ) while vertical derivatives (∂z) were computed in wavenumber domain using the fast Fourier transform using the definition below (Blakely, 1995):(17)

where g is the gravity field and the other terms were defined previously. The NFG technique was first used for exploring oil reservoirs (Berezkin, 1967; Berezkin & Buketov, 1965) and then used for various investigations (e.g. Aghajani, Moradzadeh, & Zeng, 2011; Ekinci & Yiğitbaş, 2012a, 2012b; Oruç & Keskinsezer, 2008b; Zeng et al., 2002). The normalisation procedure in NFG computation helps to minimise the strong field perturbation associated with downward continuation of the potential field data (Aghajani et al., 2011). Application of NFG technique typically enhances the anomalies by computing the anomaly to a level close to the source bodies and (Oruç & Keskinsezer, 2008b). Therefore, we produced the NFG map (Figure 14) using the residual gravity data-set in order to monitor downward extent of the high amplitude gravity anomalies. The depth level of more than 1000 m was not computed due to the oscillations and undesired effects which distorted the real anomalies.

Figure 14. NFG depth slice map of the data shown in Figure 10. The overlapped active faults in Biga Peninsula are modified after Duru et al. (2012). The overlapped faults in Gelibolu Peninsula are modified after Yaltırak et al. (2000). Black lines show coastline. Abbreviations used for the faults are given in Figure 5.

4. Results and discussion

The complete Bouguer gravity anomaly map of the study area shown in Figure 6 clearly displays both high and low amplitude anomalies. The north–south-trending high amplitude gravity anomaly zone which extends massively along the western coasts of the Biga and Gelibolu peninsulas is observed from the complete Bouguer gravity anomaly map. Some other high amplitude gravity anomalies located at the surrounding of Bandırma and at the east of Lapseki are also detected (Figure 6). Additionally a high amplitude anomaly zone exists between the ÇBFZ and SAF. The southern and eastern edges of that anomaly zone are clearly bordered by SAF and a north-west–south-east-trending fault, respectively. The north-east–south-west-oriented lineament located in the Gelibolu Peninsula dramatically cut through the Peninsula and also indicates a sharp gravity contrast between the high and low amplitude anomalies (Figure 6). Except for the low amplitude anomalies traced on the middle part of the Çanakkale Strait, north of Sinekli Fault (SF) and the vicinity of Etili district, complete Bouguer gravity anomaly characterisation seems to be mostly low (negative) which stretches from the Evciler Fault (EVF), SAF and SF to eastwards (Figure 6). The lowest gravity anomaly pattern in the complete Bouguer gravity anomaly map is observed at the south of Havran–Balıkesir fault zone (HBFZ) which is near the province of Balıkesir. On the basis of those anomaly zones the general gravity amplitude variation observed from the complete Bouguer gravity anomaly map (Figure 6) could be simply characterised by a decreasing trend towards the eastern parts of the study area i.e. inland. As is well known, gravity anomalies observed over the continental areas show generally negative gravity signatures due to the effect of thick crust, and oceanic regions where the crust is very thin are mostly characterised by positive gravity anomalies (e.g. Lowrie, 2007). In the same manner, the complete Bouguer gravity anomaly map of the study area (Figure 6) demonstrably indicates that most of the land areas are represented by low amplitude (negative) gravity anomalies while high amplitude (positive) gravity anomalies mostly dominate the coastal and marine environment. Furthermore, as also explained previously by Ateş et al. (1999), the gravity amplitude variation mentioned above is quite likely due to the effect of the isostatically thickened continental crust which gradationally increases in thickness towards the eastern regions.

Since the complete Bouguer gravity anomalies intensely represent the effects of the longer wavelength fields caused by the deeper and broader geological causative sources (Mallick, Vasanthi, & Sharma, 2011), the regional background (Figure 9) was computed by using element shape functions of FEM. In order to amplify the responses of the relatively shallower causative geological units (short wavelength anomalies) and to attenuate the crustal effects (long wavelength anomalies), a residual gravity anomaly map (Figure 10) was produced by removing the regional background (Figure 9) from the complete Bouguer gravity field (Figure 6) as mentioned in previous section. The computed regional background (Figure 9) indicates a remarkable gravity discrepancy ranging between −15 and 57 mGals, and also exhibits an amplitude variation which can be characterised by an eastwardly decreasing trend in gravity values as expected. Figure 10 displays that the computed residual gravity anomalies indicates a clear gravity contrast in the study area which could be regarded as sufficiently clear indicator to distinguish the effects of causative geological sources. High amplitude anomalies on the study area are located between Biga–Çan–Yenice–Gönen, Manyas–Balya–İvrindi–Balıkesir, Çanakkale-south of Ayvacık and the surroundings of Bandırma region (Figure 10). On the contrary to the complete Bouguer gravity anomaly map shown in Figure 6, the western coasts of Biga Peninsula demonstrate a low amplitude gravity anomaly characterisation in the residual anomaly map (Figure 10). When considering the geological map of the study area (Figure 3), it is seen that alluvium and Bayramiç formation (conglomerate, sandstone, claystone, siltstone, etc.) are mostly prevalent at the Aegean coasts of the Biga Peninsula. Thus, it can be mentioned that those sedimentary units having relatively lower densities evidently contribute to the decreasing of the amplitudes of gravity anomalies. Additionally, the residual anomaly map (Figure 10) points out that the aforementioned low amplitude gravity responses become wider in the east–west direction towards northern parts of the Aegean coasts of the Biga Peninsula. The behaviour of those anomaly signatures fits quite well with the surface geology map of the Peninsula (Figure 3), because the geological map displays that Pliocene deposits (Bayramiç formation) and upper Miocene deposits (Çamrakdere member: sandstone, calcarenite, marl and siltstone) are wider in the east–west direction at the south of the entrance to the Çanakkale Strait in comparison to Ezine coasts. Similarly, the presence of alluvium and Bayramiç formation units give rise to the lowering of the gravity amplitudes at the coasts of Biga district. On the other hand, although the sedimentary cover (alluvium and Bayramiç formation) exists at the east of the Biga coasts (towards Bandırma district), gravity anomalies exhibit relatively higher amplitudes. One possible interpretation of this case is that the sedimentary covers overlying the Şapçı volcanics (andesite and pyroclastic rocks) and its correlative units in that region is relatively thinner in comparison to those of western parts.

The lowest amplitude residual gravity anomalies in the study area are observed in the eastern and southern parts of the Gelibolu Peninsula. A north-east–south-west-trending nearly longitudinally sharp gravity contrast is apparently recognised in the Peninsula (Figure 10). When bearing in mind the surface geology map of the Gelibolu Peninsula (Figure 3) and the residual gravity anomaly map of the study area (Figure 10) together, it can be stated that the sharp gravity contrast with a north-east–south-west-trending lineament is associated with the response of the Anafartalar thrust fault which continues throughout the Gazhanedere formation (continental conglomerate, sandstone, mudstone, etc.). The anomaly signature of the Anafartalar thrust fault can be easily followed from the THD map (Figure 11), maxspot map (Figure 12). Profile curvature anomaly map (Figure 13) also demonstrates the position of the thrust fault by producing zero contour values. Additionally, the anomaly trace of the westward fold of the Anafartalar thrust fault at the south of the Salt Lake is dramatically enhanced by computing the maximum horizontal gradient magnitudes (Figure 12). NFG depth slice anomaly map (Figure 14) displays a wider strong anomaly having high amplitudes from the surrounding of the Anafartalar thrust fault towards inland of the Biga Peninsula. Moreover, as clearly seen from Figure 10 that the lowest gravity anomaly amplitudes dominate the vicinity of the Çanakkale Strait. Thus, it may be suggested that the thickest sedimentary deposits having low densities are located at the vicinity of the Çanakkale Strait in the study area.

A low amplitude gravity anomaly characterisation is clearly seen at the east and south-east of Balıkesir province (Figure 10). Geological map of the study area (Figure 3) exhibits that alluvium and Soma formation (conglomerate, sandstone, claystone, limestone, etc.) exist at the east and south-east of Balıkesir province. Thus, those sedimentary deposits having lower densities in comparison to bedrock cause a decrease in gravity values. On the other hand, Bornova flysch (ophiolitic rock and sandstone with limestone blocks, claystone and conglomerate), Yayla melange (serpentinite, spilitic basalt, gabbro, diabase, sandstone, claystone, etc.) and Hallaçlar volcanics (andesite, dasite, pyroclastic rocks, etc.) produce higher amplitude gravity anomalies at the surrounding of Balıkesir province due to their relatively higher densities (Figure 10). THD anomaly map (Figure 11) and, in particular, maxspot map (Figure 12) clearly indicates the trace of the boundaries of the sedimentary deposits in the south of HBFZ around Balıkesir province. Maxspot map (Figure 12) also indicates that the northern and western edges of the sedimentary deposits are characterised by larger maximum horizontal gradient amplitudes. This finding is most likely due to the abrupt lateral changes in densities between the sedimentary deposits and the high density Neo-Tethyan ophiolites and melange (Bornova flysch and Yayla melange) that form the northern edges of the İzmir–Ankara suture zone. Moreover, the spatial locations of the edges of sedimentary deposits on the grid plane correspond to zero contour value in profile curvature anomaly map (Figure 13) which is in well agreement with the maximum horizontal gradient locations displayed in Figure 12. The existence of downward extend of this anomaly zone is also shown in NFG anomaly map (Figure 14) and it must be noted that the anomaly zone appears to be wider in the image. Since the residual gravity anomaly map (Figure 10) demonstrates that the anomaly of that sedimentary deposits are not characterised by very low values, it may be mentioned that the deposits are not quite thick. Thus, the wider anomaly zone seen in NFG image map (Figure 14) might be caused by the high-density Neo-Tethyan ophiolites and melange which cover larger area in the subsurface than as exposed in the ground surface. Additionally, it may be stated that the sedimentary deposits overlie the Neo-Tethyan ophiolites and melange in that region.

One of the most remarkable high amplitude gravity anomaly zone is located between the districts of Biga–Çan–Yenice–Gönen (Figure 10). The boundaries of that high amplitude gravity anomaly zone can be clearly observed from the maxspot map (Figure 12). Zero contour values of profile curvature anomaly map (Figure 13) are also compatible with the locations of maximum horizontal gradients. The northern boundary of that zone located nearly on the eastern parts of ÇBFZ produces the largest maximum horizontal gradient magnitudes (Figures 11 and 12) and clearly exhibits the vertical tectonic contact between Işıkeli rhyolite (rhyolite, perlite, ignimbrite and tuff) and Pliocene deposits (Bayramiç formation). This high amplitude anomaly elongated from north-east to south-west (Figures 11 and 12) appears to be important and could be the continuation of the Main Fault Zone (MFZ). MFZ has been previously defined in the south shore of the Marmara Sea by Ateş et al. (2009) based on the analyses of aeromagnetic anomalies. Yiğitbaş, Yılmaz, and Şengör (1995) and Yiğitbaş et al. (1999) have also suggested a fault zone in that region. Additionally, this lineament is consistent with ‘Western Pontide Fault’ which has been defined by Elmas and Yiğitbaş (2001). In the south of this anomaly zone, some low and moderate values of maximum horizontal gradients (Figure 12) indicate the contact between the Oligo–Miocene granitoid and Hallaçlar volcanics (andesite, dacite, pyroclastic rocks, etc.) (Figure 3) due to detectable density discrepancy. However, the contact between the Oligo–Miocene granitoid (Ezine–Kestanbol granitoid) and Hallaçlar volcanics located between the Ezine and Ayvacık districts could not be detected. This case may be explained by the alteration on the Hallaçlar volcanics which might have reduced the density difference between those igneous rocks. The surface geology map (Figure 3) shows that Triassic units (Karakaya complex-undifferentiated, Camialan limestone, Çal formation, Mehmetalan formation, Arkosic sandstone and Orhanlar greywacke) are prevalent at that high amplitude gravity anomaly zone between the ÇBFZ and SAF. However, these Triassic units do not seem to produce such a high amplitude gravity response as shown in Figure 10 due to their relatively lower densities. Therefore, that high amplitude anomaly zone are most likely caused by the relatively higher densities of the Işıkeli rhyolite, Şapçı volcanics and Şahinli formation (basalt, pyroclastic and volcanoclastic rocks) which are overlaid by Triassic units between the ÇBFZ and SAF. Maxspot (Figure 12) and NFG anomaly (Figure 14) maps indicate that the high anomaly zone are bordered by the north-east–south-west-trending faults (ÇBFZ) at the south of the Biga district. Additionally, Figures 12 and 14 demonstrate that the location of the eastern boundary of that anomalous zone correspond to a north-west–south-east-trending fault formed at high-angle to the main fault systems (ÇBFZ and SAF). The existence of another north-west–south-east-trending fault formed at high-angle to the fault system is clearly seen at the western part of the anomalous zone. Thus, the high amplitude gravity anomalies may be supposed to be the trace of a transtensional or transpressional fault zone in the strike-slip fault systems. Additionally, the high amplitude gravity anomaly zone observed through the districts of Bandırma–Biga–Çan (Figure 14) probably indicates a significant shear zone at that region. However, more detailed geophysical and geological studies should be carried out in that region to strengthen those interpretations.

Another noteworthy feature observed on the residual gravity anomaly map (Figure 10) is located at the western parts of the Biga Peninsula. That significant high amplitude gravity response extends from the south of the Çanakkale province to the Ayvacık coasts. Interestingly, low amplitude gravity signatures arising from the effects of alluvium and Bayramiç formation are located at about the centre of that high amplitude gravity anomaly zone. THD (Figure 11), maxspot (Figure 12) and profile curvature anomaly (Figure 13) maps dramatically exhibit the horizontal locations of the edges of that nearly circular-shaped anomaly zone. Additionally, NFG depth slice (Figure 14) also displays the possible edges of that causative geological source at the subsurface. The highest amplitude gravity anomalies in that zone are observed at the south-west of the Bayramiç district (Figure 10). The derivative-based anomaly maps (Figures 11–13) also indicate the north-east–south-west-trending highest amplitude anomaly at the south-west of the Bayramiç district. That nearly circular-shaped gravity anomaly zone might be the trace of a large caldera existed during the active magmatism period in the Biga Peninsula. If this assumption is true, then it is possible to expect that the north-east–south-west–trending high anomaly trace located at the centre of the caldera indicates a volcanic eruption centre or a volcanic neck where volcanic rocks younger than the surrounding rocks exist. Indeed, a volcanic neck representing a volcanic eruption centre can be seen at the east of the Ezine–Ayvacık highway, and this instance strongly corroborates the assumption of the existence of a large and old caldera. Additionally, it must be noted that this old caldera structure might have been significantly altered or covered by neotectonic units and sedimentation in the region. Considering the geology of the region (Figures 3 and 4), it may be mentioned that the Pliocene deposits (Bayramiç formation) might be preserved at about the centre of the caldera structure. Moreover, the north–south-trending incongruous end of the western boundary of Pliocene deposits in proportion to the general structure of the area between Ezine and Bilaller (a red square at the north of Ayvacık in Figure 3) espouses that assumption. As is known, instead of flowing into the Aegean Sea by passing through Ezine–Kestanbol (a blue square at the west of Ezine in Figure 3)-Dalyan (a green square at the west of Ezine in Figure 3), Karamenderes River has changed its route northward at the vicinity of the Ezine district. One of the possible interpretations that can be offered for that case is that the Bayramiç basin (Figure 4) might have been in the same form as it now appears namely a depression surrounded by ridges. This depression may be attributed to a caldera over a shallow intrusive magma chamber. Moreover, provided that a coeval volcanism was developed at the surface and the supposed volcanic eruption centre is represented by the present topographic Bayramiç depression area, it can be stated that the nearly circular-shaped gravity anomaly zone (Figures 10–14) indicates the remains of an old caldera structure. However, a detailed study of the volcanic rocks distributed from the volcanic centre should be carried out to verify the existence of the old caldera. To the best knowledge of the authors, no special investigation has been yet performed on those volcanic rocks in such detailed manner which may provide essential data for the study presented here. Another possible interpretation is that the Bayramiç basin (Figures 3 and 4) was a north-east–south-west-trending depression area that extend over the sea but later on it was blocked by north–south-trending faults. Therefore, Kestanbol has been eroded by uplifting, and Karamenderes River has changed its route northward. Above-stated scenarios seem possible for explaining the meaning of nearly circular gravity anomaly. Additionally, if the aforementioned gravity anomaly zone (Figures 10–14) does not represent the trace of the remains of a large caldera structure, the existence of a magma chamber (may be granitic) under the Bayramiç depression gives rise to think about a cauldron subsidence representing this depression, which should be supported by detailed geological data in the region.

5. Conclusions

We analysed the Bouguer gravity anomalies in a detailed manner to delineate shallow subsurface geology and some structural features of Biga and Gelibolu peninsulas and also their surrounding area. Gravity data analyses were carried out by executing a developed MATLAB-based toolkit named as GMINTERP which is composed of a set of linked functions in conjunction with a graphical user interface. The coded algorithms in GMINTERP software that are incorporation into a large mathematical library yielded significant geophysical anomaly maps and also brought out some subtle details which shed light on the geological interpretation of the studied region.

A finite element algorithm based on the use of element shape functions was used to produce the residual gravity anomaly map which is mostly compatible with the surface geology map of the study area. In order to better understand the nature and geological meanings of the gravity anomaly zones, some derivative-based edge detection techniques were performed to residual gravity data-set. Implementation of these techniques yielded to locate the abrupt lateral changes in densities on the grid plane and the findings obtained lead to a sound some geological interpretations about tectonic and magmatic structures in Biga and Gelibolu peninsulas. The most noticeable feature detected in the anomaly maps derived from different data processing techniques is the existence of a nearly circular-shaped gravity anomaly signature located at the western part of the Biga Peninsula. In consideration of some geological knowledge, this high amplitude gravity anomaly traces were attributed to represent a large and old caldera structure which might have existed during the active magmatism period in the Biga Peninsula. The high amplitude anomaly zone located at the centre of that anomaly zone strengthened the interpretation, and was supposed to indicate a volcanic eruption centre or a volcanic neck. Additionally, the existence of a volcanic neck-like structure at the east of the Ezine–Ayvacık highway may be accepted as an evidence of a caldera structure. On the other hand, it was thought that the presence of a cauldron subsidence representing the Bayramiç depression is within the bounds of possibility, but it was concluded that more detailed geophysical and geological studies would be able to show whether that gravity anomaly traces represent an old caldera or a cauldron subsidence or any other geological structure.

This study clearly highlighted the key role of the potential field methods such as gravity in better understanding the structural features of the complex geological areas and to reveal some significant geological implications. Additionally, the findings obtained indicated that the developed interactive data processing toolkit proved useful in analysing the potential field data-sets and it may effectively aid geological interpretation.

Acknowledgments

This paper is a part of PhD dissertation undertaken by the corresponding author (YLE) at Graduate School of Natural and Applied Sciences in Çanakkale Onsekiz Mart University, under the supervision of co-author (EY). Can Ertekin is warmly thanked for his great effort in digitising the 1/250,000 scaled surface geological map of the Biga Peninsula and its surrounding area. Thanks are due to Dr Abdullah Ateş, Dr Aydın Büyüksaraç, Dr Ömer Feyzi Gürer and Dr Mustafa Bozcu for their contributions. The location map in Figure 5 was generated by using Generic Mapping Tools (GMT) (Wessel & Smith, 1995) and Alper Demirci is thanked for his help in preparing this map. Two anonymous reviewers provided many helpful comments that clearly improved our paper.

Disclosure statement

No potential conflict of interest was reported by the authors.