Crustal thickness variations and isostatic disequilibrium across the North Anatolian Fault, western Turkey

We use teleseismic recordings from a dense array of seismometers straddling both strands of the North Anatolian Fault Zone to determine crustal thickness, P/S velocity ratio and sedimentary layer thickness. To do this, we implement a new grid search inversion scheme based on the use of transfer functions, removing the need for deconvolution for source normalization and therefore eliminating common problems associated with crustal‐scale receiver function analysis. We achieve a good fit to the data except at several stations located in Quaternary sedimentary basins, where our two‐layer crustal model is likely to be inaccurate. We find two zones of thick sedimentary material: one north of the northern fault branch, and one straddling the southern branch. The crustal thickness increases sharply north of the northern strand of the North Anatolian Fault Zone (NAFZ), where the fault nearly coincides with the trace of the Intra‐Pontide Suture; the velocity ratio changes across the southern fault strand, indicating a change in basement composition. We interpret these changes to indicate that both strands of the NAFZ follow preexisting geological boundaries rather than being ideally aligned with the stress field. The thick crust north of the northern NAFZ strand is associated with low topography and so is inconsistent with simple models of isostatic equilibrium, requiring a contribution from mantle density variations, such as possible loading from underthrust Black Sea oceanic lithosphere.


Introduction
The formation and evolution of active interplate faults is a balance between accommodation of the required strain field and exploitation of preexisting weak zones from past geological processes. How this balance is resolved in particular cases and in three dimensions is the source of considerable debate. A good example of this type of setting is the North Anatolian Fault Zone (NAFZ), a ≈1500 km long right-lateral strike-slip system which forms the northern boundary of the Anatolian Plate [Salah et al., 2007]. The NAFZ is a young geologic feature that formed by westward propagation, reaching western Turkey circa 200 ka [Şengör et al., 2005]; its seismicity also displays a westward progression and represents a major source of seismic hazard. In western Turkey, the NAFZ divides into northern and southern strands [Karimi et al., 2014] whose relationship is not entirely clear. The northern strand is currently the more seismically active of the two and was the location of two devastating earthquakes nearİzmit and Düzce in 1999 [Barka et al., 2002;Gülen et al., 2002].
Turkey consists, for the most part, of an amalgam of continental fragments and oceanic remnants that assembled in the late Tertiary, due to the closing of the Tethyan Ocean and the Alpide orogenic system [Şengör and Yilmaz, 1981;Okay, 2008]. The NAFZ approximately parallels the resulting east-west tectonic fabric; in western Turkey, the NAFZ follows the trend of Tethyside accretionary complexes [Şengör et al., 2005], with the northern branch lying close to the Intra-Pontide Suture (IPS) between the accretionary Sakarya Terrane and theİstanbul Zone (a continental fragment with Late Proterozoic basement) [Okay, 2008]. Sedimentary cover of variable age and thickness is prominent in western Turkey. Notably, extension along a number of east-west trending grabens resulted in thick deposits of Neogene sediments [Sari anḑ Salk, 2006], in places up to 4 km thick.
Researchers from the University of Leeds, the Kandilli Observatory, and Sakarya University installed a dense array of broadband seismometers in western Turkey from May 2012 through September 2013, FREDERIKSEN ET AL.

Figure 1.
Study area and instrumentation. The SEIS-UK Dense Array for North Anatolia instrumentation (circles) is supplemented by temporary stations KO01 through KO07 (diamonds) and permanent stations GULT, SAUV, and SPNC, operated by Kandilli Observatory, Bogaziçi University (triangles). Background shading is topography from the NASA Shuttle Radar Topographic Mission [Jarvis et al., 2008]. Thick black lines are fault traces; dashed black lines are major sutures [Okay, 2008], and red lines indicate the rupture zone of the 1999 Izmit earthquake [Barka et al., 2002]. Thin black lines outline Quaternary basins [Oguz and Sasatani, 2004;Dogan et al., 2014]. in order to examine crustal and mantle structures associated with the NAFZ. The array consisted of 66 stations with a nominal station spacing of 7 km (DA01 through DF11; Figure 1) and included three stations of the Kandilli Observatory and Earthquake Research Institute network (GULT, SAUV, and SPNC). The stations were deployed on six north-south trending lines, each containing 11 stations. This dense array was supplemented by seven stations forming a semicircle around the east side of the dense grid (KO01 through KO07). The array was deployed across the rupture of the 1999İzmit (Kocaeli) earthquake [Barka et al., 2002] and crosses both surface strands of the NAFZ in this region. The array straddles the IPS and samples three main crustal blocks: theİstanbul Zone in the north and the Sakarya Terrane in the south, separated by the Armutlu-Almacık Block in the center of the array.

Data and Analysis
From the 16 months of recorded data, teleseismic events with magnitudes greater than 5.5 were visually examined for P wave signal-to-noise ratio on the vertical and horizontal components. A total of 21 high-quality events were retained for further analysis, from which we retained a time window stretching from 10 s prior to the IASP91-predicted P wave arrival to 30 s after the predicted arrival. A few traces were removed at individual stations due to noise problems or were absent due to data loss. The retained time windows are dominated by the P arrival and its associated coda (see, e.g., Figure 2b).
Crustal thickness is commonly measured from the teleseismic P coda using the popular H-k stacking method [Zhu and Kanamori, 2000]. This is a two-stage process in which receiver functions are generated by deconvolution of (most commonly) the vertical component from the radial component, followed by stacking along predicted arrival-time curves for Ps conversions and reverberations for a grid of assumed Moho depths (H) and crustal P/S velocity ratios (k). The total stacked amplitude is expected to reach a maximum when the predicted arrival-time curves are most accurate, i.e., when the Moho depth (H) and P/S velocity ratio (k) values are closest to being correct.
The assumptions of H-k stacking are violated in the presence of complex crustal structures, particularly in cases where strong intracrustal velocity contrasts are present, as is often the case in sedimentary basins. Sedimentary basins can introduce a strongly oscillatory character to receiver functions [Zelt and Ellis, 1999], which can complicate H-k analysis and yield inaccurate or ambiguous results [Yeck et al., 2013]. In lieu of H-k stacking, we apply a recently developed transfer function technique which avoids deconvolution and allows for the inclusion of a sedimentary layer in the analysis.
The transfer function technique involves using ray theory to predict the transfer function T zr , which relates the vertical to the radial component of teleseismic data, for a given assumed Earth model. If the radial and vertical Green's functions of the receiver-side structure are represented by G r and G z , and the source waveform is S, the recorded teleseismic traces are, in the frequency domain, R = SG r and Z = SG z for the FREDERIKSEN ET AL. radial and vertical components, respectively. The transfer function is then simply the spectral ratio of the Green's functions: and the relation between the radial (R) and vertical (Z) components is given by (2) where the source contribution S is canceled out. For a given Earth structure, we can use ray theory to generate a transfer function and then use it to predict the radial component from the vertical. This approach makes it possible to determine the consistency of a data set with a range of models without needing to deconvolve, in similar fashion to the cross-convolution approach of Bodin et al. [2014]. It also becomes possible to include a multilayered crust in the analysis.
We used the transfer function technique to perform a grid search over three parameters: the total crustal thickness (H), the P/S velocity ratio of the basement (k), and the thickness of an overlying sedimentary layer. We held the sedimentary P and S velocities fixed, in this case at 5.0 km/s and 2.9 km/s, respectively, in line with the values used by Salah et al. [2007] in the same area. Crustal thickness was varied from 30 to 45 km at 0.5 km intervals, the range being selected to be in line with past studies of crustal thickness in western Turkey [Mutlu and Karabulut, 2011;Tezel et al., 2013;Vanacore et al., 2013]. The P/S velocity ratio was allowed to vary from 1.6 to 1.9 at intervals of 0.02; the sediment thickness range used was 0 to 5.5 km at 0.5 km intervals [Sari anḑ Salk, 2006;Tigli et al., 2012]. The misfit we used was the variance of the difference between the predicted and measured radial traces: allowing us to calculate error bars using the F test approach commonly used in shear wave splitting analysis [Walsh et al., 2013]. The square root of this value (root-mean-square misfit) is given in the table and plots.
A sample grid search result for a good quality station (DA01) is shown in Figure 2; the misfit (a root-meansquare measure) obtained in this case was 0.0902, with a recovered sedimentary thickness of 3.0 km.
FREDERIKSEN ET AL.
©2015. The Authors.  Figure 3a shows the variation in misfit across the instrument array; stations in black were removed from the data set due to misfits exceeding a threshold of 0.18, which were often accompanied by an obvious visual mismatch between the observed and predicted radial components (generally in the form of a major difference in amplitude). Figures 3c and 3d show the error bars on crustal thickness and velocity ratio, which are generally larger for noisier stations. The majority of these high-misfit stations lies within the Quaternary Adapazarı and Pamukova basins, where the near-surface velocity is very low [Komazawa et al., 2002] and our two-layer model is likely inadequate, particularly given the considerable basement topography under these basins. Our assumed sedimentary velocities are more appropriate for older basins containing sedimentary rock than for unconsolidated Quaternary material.
In order to test the sensitivity of the results to variations in sediment velocity, we ran three additional analyses with different P velocities (3.5, 4, and 4.5 km/s) in the uppermost layer, with the P/S ratio fixed at 1.73. The final misfit was insensitive to the sediment velocity, varying on average by 4% between the best and worst values. We found that the sediment thickness trades off with velocity (for example, DA01's recovered sediment thickness ranges from 2 km at low sediment velocities to 3 km at the 5 km/s velocity used in the final results) and is therefore only interpretable in a relative sense. The overall crustal thickness is similarly affected, though by a proportionally smaller amount (32.5 to 34.5 km at DA01), and its overall pattern is not changed. The pattern of basement P/S ratio shows changes at low sediment velocities (3.5 or 4 km/s) and so should be interpreted with some caution. We conclude that the teleseismic data are sensitive to aspects of the near-surface structure that they cannot constrain; our preferred approach is therefore to retain a fixed value of 5 km/s for the near-surface velocity and to bear in mind that it may be a possible source of bias.

Results
Though sediment thickness is not a parameter that we can recover with great accuracy, due to both our assumption of constant sediment velocities and the frequency-content limitations of the data set, the pattern we observe (Figure 3b) is nonetheless of some interest. Only two stations were best fit by a one-layer (no sediment) model, with all other stations requiring sedimentary basins of 1.5 to 5.5 km in thickness. The basin thickness pattern is spatially coherent, with two zones of thick (≥ 3.5 km) sediment: one north of the north branch of the NAFZ, and one straddling the southern branch. These zones do not correspond exactly to the surface locations of the Quaternary basins, though there is some overlap. The northern zone also matches a zone of low velocity in the shallow crust imaged by local earthquake tomography [Koulakov et al., 2010]. Given that our assumed sediment velocity is higher than that expected in the Quaternary basins and that we are unable to achieve an acceptable fit at many of the basinal stations, we believe that the sedimentary layer included in our models reflects older, more consolidated material.
The crustal thicknesses we observe (Figure 4a) range from 30 to 45 km. The overall pattern shows a thickening of the crust from south to north, with a sharp increase (≈ 7 km) corresponding to the northern strand of the NAFZ, though the thick crust is also truncated on its eastern edge at about 30 • 25'E by thinner crust. A small patch of thick crust is present along the southern fault strand at about 30 • 15'E, apart from which all stations south of the northern strand exhibit crustal thicknesses of less than 40 km.
The pattern of basement P/S velocity ratio variation (Figure 4b) shows more scatter than the crustal thickness, perhaps in part because the velocity ratio is less well resolved (note the shape of the low-misfit region in Figure 2a). The simplified two-layer model we are using does not completely capture the complexity of the sedimentary deposits in our study area, potentially leading to sedimentary contributions in some of these values. Despite this scatter, there is a significant increase in the velocity ratio south of the southern fault strand. North of the southern strand, the velocity ratio is generally low (1.6-1.7) except for isolated stations and for a localized zone at approximately 40 • 50'N, 30 • 10 ′ W.

Discussion
The primary observation of this study is that the crustal thickness increases sharply north of the northern strand of the NAFZ, while the P/S velocity increases sharply south of the southern strand. Although previous receiver function studies of crustal thickness and P/S ratio [Tezel et al., 2013;Vanacore et al., 2013] did not have sufficient station coverage to detect sharp changes across the fault strands, they detected northward crustal thickening in our study area. A crustal thickness study based on gravity measurements [Arslan et al., 2010], however, interpreted the crustal thickness to decrease northward across the NAFZ assuming Airy FREDERIKSEN ET AL.

10.1002/2014GL062401
isostasy. The large thickness jump we observe at the northern strand of the NAFZ also corresponds to the location of the Intra-Pontide Suture (IPS) and so may predate the NAFZ; it is also possible that the fault's trajectory was influenced by the presence of the suture.
Though our P/S velocity ratio measurements are more scattered (and likely have greater uncertainties) than our thickness measurements, they indicate a change in velocity ratio (and so presumably crustal composition) across the southern strand of the NAFZ. The tomographic models of Salah et al. [2007] and Koulakov et al. [2010] also show significant P/S velocity ratio variations in the area, though the variations they observe are less systematic. The velocity ratio generally decreases with increasing felsic content, due to the very low ratio for quartz [Christensen, 1996], indicating more mafic material south of the southern fault strand. Given that we see a compositional change, which is unlikely to be related to the fault itself, it seems likely that both strands of the NAFZ follow preexisting geological boundaries. If this is the case, the resulting imperfect alignment with the stress field would potentially produce localized uplift and pop-up structures and so affect the pattern of topography. The coincidence between the surface location of the fault strands and the changes in crustal properties that we observe suggests that both strands of the NAFZ are vertical or near vertical through much of the crust. Models of fault zone deformation [Platt and Behr, 2011] suggest that plate-boundary strike-slip faults become broad shear zones in the midcrust, which would imply that the NAFZ strands merge at depth; our results suggest that this is not the case.
A peculiarity of the thickness change that we observe is that the thickest crust is associated with low rather than high topography (Figure 4a), indicating that under the assumption of constant-density crust and no mantle variations, the crust is not in Airy isostatic equilibrium. This observation is the presumed cause of differences between our crustal thickness pattern and that determined by Arslan et al. [2010]. A likely explanation is that both of these assumptions are incorrect: the crustal density is likely to be nonuniform, given that the northern strand of the NAFZ corresponds to a suture and the southern strand corresponds to a change in velocity ratio, and the mantle is likely to contain significant density variations. Notably, Fichtner et al. [2013] used waveform tomography to resolve a linear low-velocity zone in the lithosphere below the eastern and central portion of the NAFZ, extending to a depth of 100 km, which they interpret as the lithospheric expression of Tethyan sutures (including the IPS) as well as a zone of weakness that may have played a role in determining the trajectory of the NAFZ. If this low-velocity zone represents a zone of higher temperature than the surrounding lithosphere and continues westward beneath our target area, then it may provide isostatic support for the thinner crust between the fault strands and so account for the association of thinner crust with higher topography [Flament et al., 2013]. There is also evidence that oceanic lithosphere from the Black Sea underlies the crust north of the NAFZ [Gülen et al., 2002] and may be a cause of isostatic loading and low topography.

Conclusions
We have successfully measured crustal thickness and P/S velocity ratio over a densely sampled section of the North Anatolian Fault Zone, using a new transfer function approach that avoids deconvolution artifacts and accounts for a near-surface sediment layer. Our results show changes in crustal properties at both strands of the NAFZ: an increase in crustal thickness north of the northern strand, and an increase in velocity ratio south of the southern strand. Given that the northern strand lies close to the Intra-Pontide Suture, we interpret both strands of the NAFZ to follow preexisting geological features related to Tethyan accretion and to continue vertically through the full thickness of the crust. The association of thick crust with low topography at the north end of the study area indicates that a simple Airy isostatic model of the region is insufficient and that the isostatic situation is likely to be complex, with a potentially significant influence from mantle density variations.