Secondary Fault Activity of the North Anatolian Fault near Avcilar, Southwest of Istanbul: Evidence from SAR Interferometry Observations

Abstract

Strike-slip faults may be traced along thousands of kilometers, e.g., the San Andreas Fault (USA) or the North Anatolian Fault (Turkey). A closer look at such continental-scale strike faults reveals localized complexities in fault geometry, associated with fault segmentation, secondary faults and a change of related hazards. The North Anatolian Fault displays such complexities nearby the mega city Istanbul, which is a place where earthquake risks are high, but secondary processes are not well understood. In this paper, long-term persistent scatterer interferometry (PSI) analysis of synthetic aperture radar (SAR) data time series was used to precisely identify the surface deformation pattern associated with the faulting complexity at the prominent bend of the North Anatolian Fault near Istanbul city. We elaborate the relevance of local faulting activity and estimate the fault status (slip rate and locking depth) for the first time using satellite SAR interferometry (InSAR) technology. The studied NW-SE-oriented fault on land is subject to strike-slip movement at a mean slip rate of ~5.0 mm/year and a shallow locking depth of <1.0 km and thought to be directly interacting with the main fault branch, with important implications for tectonic coupling. Our results provide the first geodetic evidence on the segmentation of a major crustal fault with a structural complexity and associated multi-hazards near the inhabited regions of Istanbul, with similarities also to other major strike-slip faults that display changes in fault traces and mechanisms.

1. Introduction

Crustal strike-slip faults are found to be highly segmented, with patches of high or low slip, and parasitic or associated faults that may sporadically become seismically active [1,2,3] and often localize at geometric bends of the main fault trace [4]. The North Anatolian Fault (NAF) is one of the most active strike-slip faults in the world and can be traced from eastern Turkey to the northern Aegean Sea with a scale of ~1200 km [5]. Near the Almacık block, the NAF separates into two branches, and the northern branch of the NAF (NNAF) continues submarine into the Marmara Sea, aligning with the north and south slopes of the Çınarcık Basin [6] (Figure 1). The strike of the northern branch then bends from the northwest to the west [6]. As we shall discuss in this work, the localization of this geometric bend is also the region of fault complexity near Küçükçekmece Lake, associated with secondary faulting processes that can be observed using synthetic aperture radar (SAR) persistent scatterer interferometry (PSI), a modern space geodetic technique.

The 1999 İzmit/Düzce earthquake disaster highlighted the high seismic hazard and revealed a westward propagation of earthquakes towards the Marmara Sea [10]. The tectonic evolution of the NNAF beneath the Marmara Sea has been widely discussed in past decades; some authors suggest that the NNAF deforms as a single shear zone [11], whereas others suggest that the crust beneath the Marmara Sea is separated by an E-W trending graben [12]. Besides, the evolution of the NNAF zone has also been interpreted as a combination of en-echelon strike-slip faults and pull-apart basins [6,13,14] associated with numerous subparallel fault scarps of unknown activity [7,15]. These investigations mostly focused on the fault activity beneath the Marmara Sea, but the possible activity of secondary faults north of the Marmara Sea remained to be studied. As the faults are shaping the land and directly affect the populated regions of western Istanbul, monitoring the fault activity in this region is of great importance.

Previous studies on seismic and bathymetric observations show evidence of the existence of secondary faults near the two lakes (Büyükçekmece Lake (BL) and Küçükçekmece Lake (KL); Figure 1), north of the Marmara Sea [8,16,17,18,19]. These secondary faults were found to be NW-SE oriented following the geographic alignment of the lake, which also show general consistency with the inland drainage pattern and the scarps of the surrounding terrain [19]. Furthermore, the agreement between the strike and location of the secondary faults and off-shore active faults indicates that the secondary faults near the lakes are not local fault systems, but related to the NW extension of the NNAF [16].

Near-field geodetic observations can provide an independent constraint on the present activity of these secondary faults. Recent GPS observations show local crust deformation (~5 mm/year) near KL [17,20], however, the inferred results may suffer from the spatial limitation of the survey-mode GPS observation. Previous studies [21,22] already described the pronounced subsidence of the Avcilar area (Figure 1), but the horizontal movement related to the secondary fault activity remained to be studied.

In this paper, well-covered long-term PSI measurements were carried out by exploiting the ESA ENVISAT ascending (Track 429) and descending (Track 336) data archive, in order to study the stable near-field crustal deformation (Figure 1). After geometric decomposition and model-based correction of the post-seismic viscoelastic relaxation associated with the 1999 İzmit/Düzce earthquakes, the slip rate and locking depth of the secondary fault beneath KL were determined for the first time based on the assumption of a single dextral fault that was inferred from geological studies [16,17]. We compare our findings to independent information and discuss the implications of these secondary faults, as well as the relationship of the fault slip and fault trend to the geomorphology and landslide activity.

2. Data Handling and Method

Spaceborne SAR interferometry (InSAR) allows investigating the study area at high resolution and identifying subtle deformation occurrences. The basic principles of SAR interferometry have been extensively described in the literature [23,24], and the techniques were then significantly improved by the application of time series procedures. The small baseline subset technique [25] allows such time series analysis and is based on interferograms with small baselines to limit geometric decorrelation and was used to map subsidence zones in the Avcilar area nearby Istanbul [22]. The PSI technique, with the idea of selecting distinct pixels, named persistent scatterers (PS), characterized by a coherent signal over a long series of acquisitions, irrespective of the baseline, brought important advances in the field, allowing millimetric precision of displacement detection [26]. Since then, different PSI and, in general, multi-interferogram InSAR techniques have been developed [25,27,28,29,30,31]. By means of PSI techniques, very long time series, such as available over the study area west of Istanbul, can be processed, and very dense PS measurements can be obtained in spite of temporal and geometrical decorrelation effects.

In this work, an advanced PSI technique named persistent scatterer pair (PSP) [30] was used to process the full ENVISAT archive of SAR data, in order to study aseismic deformations to the west of Istanbul, a region that is located in the vicinity of an important bend of the NNAF.

2.1. SAR Interferometry Processing

The PSP approach was recently proposed [29,30,32] for the identification of PSs in series of full resolution SAR images and the retrieval of the corresponding 3D positions and displacements. Differently from standard PSI techniques, the PSP method relies only on the relative signal in pairs of nearby points, both for identification and analysis of PSs. For this reason, the PSP technique is intrinsically robust to spatially-correlated disturbances (such as atmospheric and orbital artifacts) and does not need a preliminary processing aimed at compensating these contributions to the signal. In addition, deviations from the models of displacement evolution typically used in PSI techniques (for example, due to strong nonlinear evolutions) are mitigated by the PSP approach because close points likely have similar deviations. Against the mentioned advantages, working with arcs requires many more computations than analyzing single points: in fact, N points can form (N − 1) × N/2 different arcs. The PSP approach is based on a complex algorithm designed to find a set of valid (in the sense of the multi-acquisition coherence or analogous parameters) arcs to form a graph with given connectivity and cardinality. The PS points are identified as the nodes of this graph.

The high number of arcs and the high connectivity of the graph guarantee reliable PS identification and accurate PS measurements. Moreover, it makes it possible to fully extract coherent information even from low-intensity SAR signals, even when strongly scattering structures are not present, the sensed signal is weak and the scattering distributed, as in the case of rather smooth surfaces or natural terrains. In this sense, PS are intended in the general sense of scatterers whose backscattering properties persist unaltered over time (at least for the period of interest), as the name suggests, but that can be also fully distributed. In fact, thanks to the higher range resolution and the stricter orbital control, in the recent satellite SAR systems, the interferometric baseline is below the critical limit, and there are no strong spatial decorrelation phenomena.

The redundant information provided by the highly connected PSP graph is conveniently exploited to guarantee more accurate measurements by a recent method for robust phase unwrapping and finite difference integration [33], which are key processing steps in the reconstruction of PS elevations and displacements from InSAR data. The use of a multi-scale decomposition method [34] allows one to perform a global processing approach also in large areas with high PS density, despite the high computational demand required by the exploitation of very redundant information.

The advantages of globally exploiting redundant arcs connecting different pairs of PS can be particularly relevant in order to obtain reliable and accurate measurements when the PS are very sparse, as in the considered region of Istanbul, which is characterized by the presence of water bodies, mountains and large vegetated areas.

For this study, an area of about 600 km² surrounding KL was analyzed exploiting the ESA ENVISAT ascending (Track 429) and descending (Track 336) data archive. The ascending dataset consisted of 33 images covering the period December 2002–April 2010 and the descending dataset consisted of 32 images covering the period April 2003–September 2010. The images acquired from the ENVISAT ASAR sensor have the following characteristics: strip map image mode, a resolution of 5 m × 25 m per pixel, VV polarization, a nominal incidence angle of 23 degrees and a maximum baseline of 1800 m. The 3-arcsec SRTM DEM was used as an initial guess for the removal of the topographic contribution. During the PSP processing, the elevation information was further refined from the SAR data stack, by estimating in correspondence of each PS the correction to be applied to the input DEM. The PSP (and, in general, PSI) processing makes it possible to reach a metric elevation precision for each PS with ENVISAT data. It is to be noted that the use of a higher resolution DEM as the initial guess could reduce the computational cost of the algorithm, but would not improve the final elevation accuracy of each PS.

For each dataset, the interferograms with respect to the first image of the series were generated (single master interferogram configuration). In fact, when working at full resolution with PSI techniques, there is no added information in the other interferograms (all the possible pairs can be exactly obtained as a linear combination of the single master interferograms). Persistent scatterers are by definition quite insensitive to spatial and spectral decorrelation effects, and large temporal and spatial baseline interferograms can be used, allowing the full exploitation of long time series of acquisitions. The diagrams of the temporal and spatial distribution of SAR images corresponding to the two analyzed datasets are shown in Figure 2.

By using the PSP method, about 240,000 and 200,000 PS points were extracted for the study area in ascending and descending geometries, respectively. Based on the statistical analysis that we performed, the standard deviation of the deformation measurement error is in the order of 3 mm, typically ranging between 1 and 5 mm, depending on the quality of the specific PS and on the type of ground cover, whereas the refined height measurements have a metric accuracy. The average measurement density is 100 PS/km² over the whole area and reaches 1000 PS/km² in the most urbanized areas.

2.2. Signal Decomposition

Due to the radar viewing geometry, ENVISAT observations are ~3-times more sensitive to vertical motion than to fault-parallel motion in this area, assuming a fault trend of N30°W and a strike-slip motion. Previous investigations show that the Avcilar district underwent clear subsidence [21,22]; this signal can overwhelm the horizontal motions induced by fault movement. In order to investigate the particularities of these deformations, we apply a signal decomposition approach that we designed for investigations related to the N30°W striking fault slip. By combining ascending and descending PS-InSAR observations, the LOS (line of sight) velocities can be projected onto any desired set of two orthogonal basis vectors under the assumption that the third component is zero [35]. Such an approximation is particularly suitable for areas around the middle of a straight fault segment as in the present case.

The time-series data of the study region suggest vertical, but also significant horizontal displacement components. In order to improve the signal-to-noise ratio, we estimate the mean LOS velocity between 2002 and 2010. Then, we projected the obtained ascending and descending LOS velocities (Figure 3) onto motions in the fault-parallel and vertical direction using a least square method [35,36]. A uniform fault-parallel direction of horizontal movement was assumed, whereas the fault-perpendicular horizontal component was not considered in the modelling. This horizontal strike-slip motion was already independently inferred from the GPS observation in this area [37], agreeing with the orientation of the KL fault that was shown by seismic investigation [18], as well as by structural studies [16].

The decomposed horizontal and vertical movement are shown in Figure 4a,b. Clear relative movement in the fault-parallel direction was observed on two sides of the KL fault, which agrees well with the tectonic background (dextral character) of this fault. We found that relative horizontal movements striking N30°W locally exceed 5.0 mm/year. This deformation pattern is in good agreement with horizontal velocities from GPS data (see Figure 4a). The quantitative comparison between the horizontal GPS velocity and the inverted horizontal InSAR velocity is shown in Figure A1, from which we can see that the InSAR velocities are generally comparable with that of the GPS observation in terms of magnitude and direction, though a difference is found on a few stations. Note that the inverted horizontal InSAR velocities on two continuous stations (AVCT and KCEK in Figure A1) fit well with that of the GPS observation. Furthermore, clear subsidence was observed in two areas (Figure 4b): the large one is located at the west of KL aligning along a pronounced valley; the small one is located at the east of the lake fault. Subsidence along KL thus is associated with the steep topography and clearly associated with the steeper slopes. Subsidence east of KL, in turn, is found on a relatively flat plain.

3. Effect of Post-Seismic Viscoelastic Relaxation

The post-seismic viscoelastic relaxation effect (PVRE) may affect crustal deformation near large earthquakes for decades [38]. Studies on post-seismic deformation of the 1999 İzmit/Düzce earthquakes suggest that the PVRE induced by these events is clear and should not be ignored in the PS-InSAR observation period [39]. In this paper, we estimated and removed the PVRE using a model-based method following [40]. The detailed model construction and PVRE estimation have been given in [40], and some key steps are described in this section.

Firstly, a three-layer viscoelastic model was set up to forward simulate this viscoelastic effect on each PS-InSAR point within the observation period (Figure 5). The model has one elastic upper crust layer and two viscoelastic layers representing the lower crust and upper mantle, for which the Maxwell rheology was applied (Figure 5).

Secondly, the lower crust and upper mantle viscosity of the Maxwell rheology play an important role in estimating the PVRE effect [38]. These parameters were inferred from the post-seismic GPS observation in a similar time period by employing a grid-search method [40].

Thirdly, we calculated the PVRE at each PS-InSAR data point using the derived optimal viscosity parameters and the numerical code provided by [41]. The results show that the mean PVRE effect on ground velocity is almost constant in the observation period. Using the optimal viscosities, the calculated PVRE is estimated to be ~2 mm/year horizontally in the study area and is therefore considerable for inverting the fault movement (Figure 4c,d). However, the PVRE reveals a smooth spatial pattern in the area around the KL fault, so that its influence on the estimation of the fault slip is not so significant.

At last, we removed the PVRE from the InSAR velocities, in order to obtain the deformation signal that is largely free of the post-seismic effect (Figure 4e,f).

To avoid the effects induced by other local deformation processes (landsliding, groundwater extraction and others), we low-pass filtered the data using a simple median filtering method. We got the velocity of each InSAR point by calculating the mean value of neighboring points within a distance (filtering length). After filtering the small-scale signal, the relative fault movement becomes more visible (Figure 6), which highlights the relative motions of the KL fault. Other small-scale processes may slightly affect the deformation pattern and play a secondary role in controlling the crustal deformation in this area.

4. Modelling and Results

After removing the PVRE, a velocity profile perpendicular to the fault orientation was constructed (Figure 7), which was used to estimate the fault slip rate and locking depth based on a widely-used 2D screw dislocation model [42]. In this model, the fault-parallel velocity ($v$) as a function of its perpendicular distance to the fault ($x$) can be drawn by:

$$v\left(\mathrm{x}\right)=\left(V/\pi \right)\times \arctan \left(x/D\right)$$

(1)

where $V$ and $D$ are the slip rate and the locking depth of the KL fault, respectively. A grid-search method was applied to search for the slip rate and locking depth over ranges of 0.5–20 mm/year and 0.5–20 km with 0.5-mm/year and 0.5-km intervals. The RMS (root mean square) misfit between the model and observation for each combination of the parameters was calculated, from which the optimal values of the slip rate and locking depth were fixed by minimizing the RMS misfit.

The uncertainty of the obtained optimal slip rate and locking depth were estimated statistically using the Monte Carlo method [43]. We firstly constructed the 1D variance-covariance matrix of the ascending and descending InSAR data using the method of [44] and generated spatially-correlated noise. Secondly, we obtained the perturbed dataset by adding the noise to the observation. Finally, we searched for optimal parameters based on each disturbed dataset and could analyze the best-fitting model. We run these steps 500 times in order to obtain a robust estimation about the uncertainty of the parameters. The distribution of the slip rate and locking depth are inferred from a Monte Carlo simulation and show a general consistence with the RMS misfit map (Figure 8). The result statistically reveals a mean slip rate of 5.1 ± 0.4 mm/year and a locking depth of 1.0 ± 0.5 km (the mean value and the one-sigma uncertainty; Figure 8), suggesting that the KL fault is active with a shallow locking depth. Therefore, a tectonically and potentially seismically active fault is located in the immediate vicinity of western Istanbul, with important implications as discussed in the following section.

In order to test the effect of possible small-scale deformations (such as landsliding), we also inverted the low-pass filtered data using the same method mentioned above. As shown in Figure 8, we found that similar results as those inverted from the unfiltered data could be obtained, although the data misfits became smaller as the filtering length increasing. The estimated locking depths tend systematically to a very small value (<1 km), suggesting possible creep of this secondary fault. Similar fault creep and shallow locking depth have been found on other strike-slip fault, such as the central section of the San Andreas fault [45] and western sections of the NAF [46].

5. Discussion and Implications

5.1. Limitations

A number of limitations of the present study need to be discussed that may affect the interpretation and correct understanding of our results. Although the herein used PS-InSAR technique provides very robust results at a pixel resolution that is higher than hitherto presented results in the region, the signal investigated is still rather small. Horizontal motions in the sub-centimeter scale are very challenging to observe, which is why we have decided to exploit the average velocity rather than certain time-dependent increments of the data. The advantage of this approach is that the total displacement values become significant. The disadvantage is that temporal variations may not be interpreted. The strong seasonality of the subsidence region [22] or the temporal variations of the landslides [17] are not discussed, as is the nature of the strike-slip fault motion. In this, we cannot decipher if this inferred fault slip is a short-term behavior or a long-term process, warranting further studies. Furthermore, future studies may use the temporal variability of the different processes to investigate the different types of deformation (faulting, landsliding and subsidence).

In addition, compared with the descending data, the ascending data are less sensitive to the horizontal movement of the KL fault due to the geometric setting. This limitation may highly affect the extraction of horizontal motion if only the descending data are available. In this paper, owing to the combined analysis of the ascending and descending data, such an effect would be decreased to a certain degree. Synthetic tests were carried out in order to investigate how the horizontal and vertical deformation signal depend on the geometric setting of the ascending, descending data and the fault orientation. As the results show in Figure A2, Figure A3, Figure A4 and Figure A5, we find that the input horizontal and vertical signals can generally be well recovered with the geometric setting of the InSAR data and the fault orientation in this area. Nevertheless, more close investigations, such as continuous cross-fault GPS and strain-meter observation, are strongly required in order to obtain a more robust estimation of the fault movement.

This study is based on assumptions made from earlier studies on the KL fault mechanism. In addition, it is hard to further distinguish whether the KL fault is a simple fault or a narrow zone including a few small active faults [18] using the geodetic data shown in this paper. We assumed it to be a single dextral fault as a first order approximation.

5.2. Tectonic Implication

Geodetic results in this paper show clear dextral movement of the fault near KL, revealing the movement of an active secondary fault north of the Marmara Sea. Our results of relative motion (~5 mm/year) on the two sides of KL provide an important indication of the existence of the secondary faults in this area. Combined with previous geological and seismic results [8,16,17,18,19], this active fault perhaps can be interpreted as an extension of the northern branch of the NAF (NNAF) [16,18]. This interpretation suggests two possibilities about the evolution of the north branch of the NAF: one is that the NNAF does not bend from NW to WSW at KL and connect to the Ganos fault at the west, but extends to KL northwest [16]; the other one is that the NNAF might separate into two branches, one running northwest into the KL and the other one bending to the WSW direction.

For detailed analysis of the two possibilities, more close observations, such as more dense near-field GPS [36] and numerical simulations [47], are required. However, no matter which possibility reveals the development of the north branch of the NAF, the active fault beneath KL would play an important role in regional tectonic evolution and related seismic hazard assessment.

Fault complexities at fault bends are observed worldwide. These are commonly explained by perturbations induced by geometrical complexities in a main fault trace [4]. During multiple earthquake cycles residual stresses may accumulate and lead to the formation and dislocation of secondary faults. Whether the herein investigated N30°W trending fault on land is seismogenically active is not known. However, due to the shallow locking depth shown by our modelling result, the KL fault might be subjected to aseismic creep. This perhaps can explain that only a few earthquakes were observed in this area [48].

5.3. Local Landslides and Subsidence

Faulting processes are often associated with topography reshaping, steep erosion and landslide processes. Many small-scale landslides have been observed in the Avcilar area (upper left panel in Figure 1), which may affect the local deformation process. Landslide inventory maps of that region suggest steeper slopes to be vulnerable for flank instability. For this reason, one may suspect that the observed horizontal PS-InSAR deformation might be induced by slow landslides. Many factors can affect the development of landslides, where the terrain gradient has a dominant role in controlling the susceptibility and the final movement [9,49]. In other words, if the observed motions in this area mainly result from the landslides, they should show general spatial correlation with the terrain gradient.

With this assumption, we decomposed the LOS velocities into N120°W and vertical directions using the method mentioned in Section 2.2, where the N120°W direction is close to the horizontal component of the maximum terrain gradient in this area. We then draw the decomposed horizontal velocities along the selected profile and compared them with the corresponding terrain. As shown in Figure 9, no clear spatial correlation was found between the horizontal motion and the slope of the terrain, even when filtered data were applied. We thus conjecture that the modelled long-term N30°W PS-InSAR velocity is distinct from that of the local landslides, but mainly induced by stable tectonic movement.

The subsidence area at the west of the KL fault agrees well with the results of [21] in terms of location and subsidence rate; however, our result suggests that this subsidence shows also a high spatial correlation with the Haramidere valley, which perhaps was induced by ground water redistribution [21], soli liquefaction [50] or other unknown mechanisms, which deserve future investigations. The subsidence area east of the KL fault is thought to be mainly related to substratum material consolidation processes, as well as the human activity of this rapidly growing urban area [22].

6. Conclusions

In this paper, we use long-term PS-InSAR observation (2002–2010) to investigate crustal deformation and related secondary fault activity southwest of Istanbul. The mean LOS velocity in ascending and descending geometries was decomposed into fault-parallel and vertical direction motions, in which the post-seismic viscoelastic relaxation effect of the 1999 İzmit/Düzce earthquakes was removed based on model-based estimations. Based on the assumption of the fault behavior (i.e., fixed fault orientation and strike-slip character) inferred from geological studies, our results suggest that a secondary fault near KL is now active with a dextral slip rate of ~5 mm/year and a shallow locking depth of <1.0 km. The results show the first geodetic evidence about the existence and present activity of secondary faults north of the Marmara Sea, which has valuable implications for regional tectonic evolution and close seismic hazard assessment of the mega city Istanbul.