Interseismic strain accumulation across the central North Anatolian Fault from iteratively unwrapped InSAR measurements

The North Anatolian Fault (NAF) is a major tectonic feature in the Middle East and is the most active fault in Turkey. The central portion of the NAF is a region of Global Navigation Satellite Systems (GNSS) scarcity. Previous studies of interseismic deformation have focused on the aseismic creep near the town of Ismetpasa using radar data acquired in a single line‐of‐sight direction, requiring several modeling assumptions. We have measured interseismic deformation across the NAF using both ascending and descending data from the Envisat satellite mission acquired between 2003 and 2010. Rather than rejecting incorrectly unwrapped areas in the interferograms, we develop a new iterative unwrapping procedure for small baseline interferometric synthetic aperture radar (InSAR) processing that expands the spatial coverage. Our method corrects unwrapping errors iteratively and increases the robustness of the unwrapping procedure. We remove long wavelength trends from the InSAR data using GNSS observations and deconvolve the InSAR velocities into fault‐parallel motion. Profiles of fault‐parallel velocity reveal a systematic eastward decrease in fault slip rate from 30 mm/yr (25–34, 95% confidence interval (CI)) to 21 mm/yr (14–27, 95% CI) over a distance of ∼200 km. Direct offset measurements across the fault reveal fault creep along a ∼130 km section of the central NAF, with an average creep rate of 8 ± 2 mm/yr and a maximum creep rate of 14 ± 2 mm/yr located ∼30 km east of Ismetpasa. As fault creep is releasing only 30–40% of the long‐term strain in the shallow crust, the fault is still capable of producing large, damaging earthquakes in this region.


Introduction
The North Anatolian Fault (NAF) is a major continental right-lateral transform fault located in northern Turkey. Together with the East Anatolian Fault, it facilitates the westward motion of Anatolia, caught in the convergence zone of the Eurasian plate with the Arabian plate [McKenzie, 1972]. Since the 1939 M w 7.9 Erzincan earthquake in eastern Turkey, the NAF has ruptured in a sequence of large (M w > 6.7) earthquakes with a dominant westward progression in seismicity [Barka, 1996;Stein et al., 1997]. Stein et al. [1997] and Hubert-Ferrari et al. [2000] have interpreted this sequence to result from stress transfer along strike, where one earthquake brings the adjacent segment closer to failure.
In order to understand the role that the NAF plays in regional tectonics and seismic hazard, there have been numerous estimates of the fault slip rate for the NAF using present-day deformation measured with Global Navigation Satellite Systems (GNSS) [e.g., Straub et al., 1997;Reilinger et al., 2006;Ergintav et al., 2009] or offset geological features [e.g., Hubert-Ferrari et al., 2002;Pucci et al., 2008;Kozaci et al., 2009]. There have also been several interferometric synthetic aperture radar (InSAR)-derived estimates of the fault slip rate, which have focused on the western or eastern regions of the NAF where the InSAR coherence is better [e.g., Wright et al., 2001a;Cakir et al., 2005;Walters et al., 2011;Kaneko et al., 2013;Cakir et al., 2014;Cetin et al., 2014;Walters et al., 2014;Cavalié and Jónsson, 2014;Hussain et al., 2016].
However, slip rate estimates for the central NAF are relatively poorly constrained, with sparse GNSS data north of this portion of the fault ( Figure 1) and wide ranging geological and geodetic estimates. Geological fault slip rate ranges from as low as 5 mm/yr to as high as 44 mm/yr [e.g., Barka and Hancock, 1984;Barka, 1992;  [Kreemer et al., 2014]. The colored sections indicate previous ruptures along this section of the fault. (b) The Envisat satellite data tracks used in this study. Descending tracks are colored in red and ascending tracks in blue. Hubert-Ferrari et al., 2002;Kozaci et al., 2007;Kozaci et al., 2009], while GNSS studies estimate the slip rate for the region to a range of 17-34 mm/yr [e.g., Oral et al., 1993;Noomen et al., 1996;Ayhan et al., 2002;Reilinger et al., 2006]. Shallow aseismic slip on the fault plane, i.e., fault creep, on the central portion of the NAF was first documented by Ambraseys [1970], who observed increasing displacements of a wall that was built across the fault near the town of Ismetpasa, over multiple years. Ambraseys [1970] estimated a fault creep rate of ∼20 mm/yr for the time period [1955][1956][1957][1958][1959][1960][1961][1962][1963][1964][1965][1966][1967][1968][1969]. Since this original investigation, the fault creep has been the focus of numerous geodetic studies [e.g., Cakir et al., 2005;Kutoglu et al., 2010;Karabacak et al., 2011;Ozener et al., 2013;Cetin et al., 2014]. Cetin et al. [2014] suggested that the fault creep rate has been decaying since the first measurements in 1970 to a current steady state value of ∼6-8 mm/yr. Most previous InSAR studies in this region have only used satellite data from a single-look direction, e.g., the use of descending Envisat data by Cakir et al. [2005] and Cetin et al. [2014]. Kaneko et al. [2013] used a combination of ascending tracks from the ALOS satellite and one descending frame from Envisat track 207, limiting their observational period to 2007-2011. They suggested that aseismic creep at a rate of ∼9 mm/yr is limited to the upper 5.5-7 km of the crust, which exhibits velocity strengthening frictional behavior.
Recently, Rousset et al. [2016] used high-resolution COSMO-SkyMed satellite data spanning the time window between July 2013 and May 2014 to show evidence of periods of elevated fault creep spanning a month with total slip of 20 mm, indicating that episodic creep events may be an important mechanism producing aseismic slip.
In this study we use a more complete data set covering the entire central NAF in both ascending and descending geometries and spanning the ∼8 year time window between 2003 and 2010. We remove long wavelength trends from the InSAR data using published GNSS velocities [Kreemer et al., 2014] and deconvolve the InSAR line-of-sight velocities into fault-parallel and vertical motion.
We use simple elastic dislocation models to estimate geodetic fault slip rates and locking depths and investigate the spatial variation of fault creep along the central NAF. We also develop and apply a new iterative unwrapping algorithm that minimizes unwrapping errors during the InSAR processing.

InSAR Processing
Our data set consists of 191 Envisat images from four descending tracks (250, 479, 207, and 436) and three ascending tracks (28, 71, and 343) ( Figure 1b). Together these cover the central NAF between 31.5 ∘ E and 35 ∘ E, and span the time interval 2003-2010. Details of the processed data for each track are given in Table 1.
We focus the Envisat images using ROI_PAC [Rosen et al., 2004] and use the DORIS software [Kampes et al., 2003] to construct 494 interferograms. For each track we produce a redundant connected network of interferograms while minimizing the temporal separation between acquisitions and the spatial separation of the satellite (the perpendicular baseline) ( Figure S1 in the supporting information). We correct topographic contributions to the radar phase using the 90 m SRTM Digital Elevation Model [Farr et al., 2007] and account for the known oscillator drift for Envisat according to Marinkovic and Larsen [2013]. We unwrap the interferometric phase using a new iterative unwrapping process described in section 3.
We apply the StaMPS (Stanford Method for Persistent Scatterers) small baseline time series technique [Hooper, 2008;Hooper et al., 2012] to remove incoherent pixels and reduce the noise contribution to the deformation signal, by selecting only those pixels that have low phase noise on average in the small baseline interferograms used in the analysis.
The atmospheric contribution is often the largest source of error in radar interferograms [e.g., Doin et al., 2009;Walters et al., 2013;Jolivet et al., 2014;Bekaert et al., 2015a]. To mitigate this we estimated a troposphere correction using auxiliary data from the ERA-Interim global atmospheric model reanalysis product [Dee et al., 2011]. We use the TRAIN (Toolbox for Reducing Atmospheric InSAR Noise) software package [Bekaert et al., 2015c] to correct each individual interferogram for tropospheric noise. After removing a planar phase ramp from each interferogram, the ERA-I correction reduces the standard deviation of our tracks by 8% on average. The average reduction in standard deviation is small after correction, implying that some residual atmospheric signals remain in the interferograms after the ERA-I correction. Our final redundant small baseline networks consist of a total of 297 interferograms over the seven tracks ( Figure S1). We use these networks to calculate the average line-of-sight (LOS) velocity map for each track.
Any nontectonic long wavelength signals (>100 km), including those due to orbital errors, are effectively removed from each track when the InSAR line-of-sight (LOS) velocities are transformed into a Eurasia-fixed GNSS reference frame (details in section 4). The uncertainties on the final velocity for each pixel are calculated using bootstrap resampling [Efron and Tibshirani, 1986] and are presented at the 1 sigma level in the following work.
We calculate the LOS variance-covariance matrix of the noise for each InSAR track by computing the average radial covariance versus distance (autocorrelation) using the velocities in a 50 km by 50 km region ∼250 km to the south of the fault. This region is assumed to have no tectonic deformation and contains only atmospheric noise. We fit an exponential covariance function [e.g., Lohman and Simons, 2005;Parsons et al., 2006], C(r), as where we estimate the variance ( 2 ) and the characteristic length ( ), which give the spatial correlation of noise as a function of distance between pixels (r). Our values for each track and the center of the region used to calculate the covariance function are shown in Table 2. These covariances are used in section 5 when modeling the horizontal velocities and fault creep rates.

Method Description
Phase unwrapping is the process of recovering continuous phase values from phase data that are measured modulo 2 radians (wrapped data) [Ghiglia and Pritt, 1998]. Original 2-D phase-unwrapping algorithms unwrapped the phase of each individual interferogram independently [e.g., Goldstein et al., 1988;Costantini, 1998;Zebker and Lu, 1998]. However, a time series of selected interferogram pixels can be considered a 3-D data set, the third dimension being that of time. Hooper and Zebker [2007] showed that treating the unwrapping problem as one 3-D problem as opposed to a series of 2-D problems leads to an improvement in the accuracy of the solution in a similar way to which 2-D unwrapping provides an improvement over one-dimensional spatial methods.
Fully 3-D phase-unwrapping algorithms commonly assume that the phase difference between neighboring pixels is generally less than half a phase cycle (2 radians) in all dimensions [Hooper and Zebker, 2007]. However, due to atmospheric delays, InSAR signals are effectively uncorrelated in time, violating this assumption. Other unwrapping algorithms require the assumption of a temporal parametric function, such as a linear phase evolution in time [Ferretti et al., 2001], to unwrap the phase signals.
The standard unwrapping algorithm used in the Stanford Method for Persistent Scatterers (StaMPS) software [Hooper, 2010] uses the actual phase evolution in time to guide unwrapping in the spatial dimension without assuming a particular temporal evolution model. The phase difference between nearby pixels (double-difference phase) is filtered in time to give an estimate of the unwrapped displacement phase for each satellite acquisition and an estimate of the phase noise. This is used to construct probability density functions for each unwrapped double-difference phase in every interferogram. An efficient algorithm (SNAPHU) Zebker, 2000, 2001] then searches for the solution in space that maximizes the total joint probability, i.e., minimizes the total 'cost'.
For a connected network of small baseline interferograms, the phase unwrapping of individual interferograms can be checked for network consistency by summing the phase around closed interferometric loops [e.g., Pepe and Lanari, 2006;Biggs et al., 2007;Cavalié et al., 2007;Jolivet et al., 2011] (Figure 2). In the standard unwrapping approach used in StaMPS, any interferograms identified to have large unwrapping errors are removed from the small baseline network, which can result in loss of information and/or reduction in network redundancy. Note that some other InSAR practitioners [e.g., Biggs et al., 2007;Wang et al., 2009;Walters et al., 2011] generally do not drop badly unwrapped interferograms, but attempt to correct unwrapping errors by manually adding integer multiples of 2 to badly unwrapped regions of pixels. However, this is a time-consuming process.
In our method, we iterate the standard StaMPS unwrapping procedure while calculating the sum of the unwrapped phase around closed loops for every pixel in every interferogram, using the following equation: where UW is the StaMPS unwrapping operator, n is the number of interferograms on the path around an interferometric loop, ( i+1 − i ) are the interferometric phase values of a pixel in the interferograms created by calculating the phase difference between image i+1 and i relative to a reference point, and is the error term. The reference point is chosen to be north of the fault for all tracks. Any pixels satisfying the requirement of | | < 1 rad are defined as "error-free pixels" and are assumed to be correctly unwrapped. An error term is needed because the interferograms are multilooked before unwrapping, and so we do not expect to have perfect closure around each interferometric loop. Using = 1 is reasonable as it is well below the 2 radians required to produce unwrapping errors and allows for a small amount of closure error introduced by the nonlinear nature of multilooking. In our tests setting to 0.5 made no significant impact on the acceptance rates.
In each iteration, we keep all unwrapping parameters fixed (such as the number of interferograms and filtering) but assume that pixels identified as error-free in the previous iteration are likely unwrapped correctly and apply a high cost to changing the phase difference between these pixels in the next iteration. The StaMPS unwrapping algorithm uses the double-difference phase evolution in time to calculate a probability density function of unwrapped phase for each pixel pair in each interferogram. For interferograms where both pixels in a pair are identified as unwrapped correctly, we set the weighting to 100 times those of the other interferograms, to effectively ensure that the evolution in time is fixed. In this way, the iterative unwrapping method uses the error-free pixels as a guide to unwrapping the regions that contained unwrapping errors in previous iterations.

López-Quiroz et al.
[2009] describe a process where unwrapping is iterated on the residual interferogram after the removal of an estimate of the deformation signal while our technique iterates the StaMPS unwrapping procedure on the actual interferometric phase.

Testing the Iterative Unwrapping Procedure
We tested the new algorithm on data from Envisat descending track 207, which covers a region roughly 100 km by 400 km in central Turkey ( Figure 1b). Each iteration consists of the following steps: running the StaMPS unwrapping algorithm, determining the pixels unwrapped correctly in each interferogram using the method described above and in Appendix A, applying a high cost to unwrapping across these pixels, and rerunning the unwrapping algorithm again. We iterate this procedure 30 times. The results from standard unwrapping does not change as no modifications are made to its inputs and is represented by the straight line indicating no change in the number of error-free pixels per iteration. Figure 3 shows that the percentage of error-free pixels in the entire small baseline network increases sharply with the first eight iterations from 70% to 83%, reaching a maximum of 84% after 30 iterations; meaning that there are some unwrapping errors the method is unable to fix. This is also evident from the individual interferograms ( Figure 4), which show this same rapid increase in the percentage of error-free pixels followed by a plateau. It is clear that there are some unwrapping errors that cannot be corrected (blue colors in Figure 5) using the iterative method. However, the iterative procedure greatly reduces the total number of unwrapping errors and thus increases the InSAR coverage while minimizing errors.
After eight iterations the percentage of error-free pixels increased from 90% to 94% for track 250, from 65% to 80% for track 436, from 92% to 95% for track 479, from 83% to 87% for track 343, from 71% to 77% for track 28, and from 91% to 93% for track 71. . Total percentage of pixels in the small baseline network for descending track 207 that were identified as closed, i.e., correctly unwrapped, using our iterative unwrapping procedure (blue) and the standard unwrapping (red) algorithm. There is a rapid increase in the number of error-free pixels for the first eight iterations after which it reaches a plateau. As no modification is made to the input of the unwrapping algorithm, there is no change for each iteration of the standard unwrapping algorithm.

Interseismic Velocity Field Across the Central NAF
To investigate the pattern of interseismic strain accumulation along the fault, we decompose our full InSAR velocity field into the fault-parallel and fault-perpendicular components of motion. Following the method described in Hussain et al. [2016], we do this first by resampling our InSAR LOS velocities ( Figure 6) onto a 1 km by 1 km grid encompassing the spatial extent of all our tracks. We use a nearest neighbor resampling technique including only those persistent scatterer pixels with a nearest neighbor within 1 km of the center of each grid point. We reference each track to a Eurasia-fixed GNSS reference frame by first averaging the InSAR velocities that fall in a 1 km radius around every GNSS station within the boundaries of each InSAR track. We project the GNSS velocities into the local satellite line of sight and calculate the difference from the  Error-free pixels are identified in red, while pixels that did not close, i.e., have unwrapping errors, are in blue. The unwrapped phase for each iteration is shown in Figure S7 in the supporting information.
InSAR velocities. The vertical component of the GNSS velocities are not available on the Global Strain Rate Model website. Ergintav et al. [2009] showed that the vertical GNSS component is small and very noisy over western Turkey; therefore, we only use the horizontal velocities in our analysis. We determine the best fit plane through the residual velocities and remove this from the InSAR velocities to transform the LOS velocities into a Eurasia-fixed GNSS reference frame. This procedure is done separately for each track.
To estimate the uncertainties in the data, we calculate the RMS residual in horizontal velocities in the overlapping areas between neighboring tracks assuming negligible vertical motion ( Figure S4). The residuals are approximately Gaussian with mean values close to zero. The average RMS misfit is 5 mm/yr, which gives an empirical uncertainty of ∼4 mm/yr for the individual tracks.
For every pixel where information from both ascending and descending geometries are available, we use equation (3) to invert for the east-west and vertical components of motion following the method described by Wright et al. [2004] and Hussain et al. [2016] while taking into account the local incidence angles: where D LOS is the LOS velocity, is the local radar incidence angle, is the azimuth of the satellite heading vector, and [D E , D N , D U ] T is a vector with the east, north, and vertical components of motion, respectively. Equation (3) contains three unknowns (D E , D N and D U ), but we only have two input velocities with large differences in satellite look angle in the inversion (the ascending and descending InSAR LOS velocities). Therefore, it is impossible to calculate the full 3-D velocity field without a prior assumption. The common assumption made in previous studies is that there is no vertical motion across the region of interest [e.g., Walters et al., 2014;Hussain et al., 2016]. In our case we note that both the ascending and descending tracks are equally insensitive to motion in the north-south direction. We therefore use the smooth interpolated north component of the GNSS velocities ( Figure S5) to constrain the north-south component (D N ) in the inversion, and solve for the east-west and vertical components of motion using the InSAR LOS velocities. We calculate the fault-parallel component of the horizontal velocity by assuming that motion occurs on a strike-slip fault trending at N81 ∘ E.
Our fault-parallel velocities (Figure 7a) show the expected right-lateral interseismic motion across the NAF, with red colors representing motion to the northeast and blue to the southwest. Our estimated vertical component shows that there is little vertical motion across the NAF in this region (Figure 7b).
There is a relatively sharp change in fault-parallel velocity south of the NAF (Figure 7) that coincides with the B-B ′ profile line. We believe that this is due to a combination of postseismic deformation from the 2000 Orta earthquake (M w 6) [Taymaz et al., 2007], residual atmosphere introduced mainly from ascending track 71, and postseismic deformation from the 1999 Izmit and Düzce earthquakes.

Modeling Profile Velocities
We analyze three profiles across the fault where velocities from within 20 km are projected onto the profiles shown in Figure 7a. Walters et al. [2014] noted that there is a variation in the fault-parallel velocity away from the fault that is not due to interseismic loading but due to the proximity to the Euler pole of rotation. For example, GNSS velocities presented by Nocquet [2012] show fault-parallel velocity vectors with magnitude ∼25 mm/yr close to the NAF but ∼8 mm/yr in Cyprus roughly 800 km away from the fault. This variation is mostly due to the proximity of the Cyprus GNSS stations to the pole of rotation of Anatolia with respect to Eurasia. We use the pole of rotation calculated for Anatolia with respect to Eurasia by Reilinger et al. [2006], who estimated a rotation rate of 1.23 ∘ /Myr about a pole located at 32.1 ∘ E, 30.8 ∘ N near the Nile delta. In a Eurasia-fixed reference frame this rotation effect only applies to the region south of the NAF and corresponds to a value of rot = 0.0215 mm/yr/km or 2.15 mm/yr at a distance of 100 km from the fault.
Assuming the fault-parallel velocities far to south of the fault (>200 km) are mostly due to atmospheric noise and contain no tectonic deformation, we calculate the variance-covariance matrix of the noise using the method described in section 2, using velocities from a 50 km by 50 km region centered on 32.5 ∘ E, 39 ∘ N. The estimated variance ( 2 ) and characteristic length ( ) for the covariance function (equation (1)) is 6.35 (mm/yr) 2 and 35.8 km, respectively.
Profiles A-A ′ and C-C ′ do not cross the creeping section of the fault. For these profiles we fit a 1-D model [Savage and Burford, 1973] through the profiles where the fault-parallel velocity, v par , at a fault normal where a is a static offset.
However, profile B-B ′ crosses the creeping section of the fault. For this profile we model the fault-parallel velocity as a combination of two signals: a long wavelength signal that represents interseismic loading at rate S and locking depth d 1 , and a short wavelength signal that represents the fault creep at a rate C from the surface down to depth d 2 [e.g., Wright et al., 2001a;Elliott et al., 2008;Hussain et al., 2016].
where (x) is the Heaviside function.
We find best fit values for each model parameter (S, d 1 , C, and d 2 ) and an offset a, using a Bayesian approach, implementing the Goodman and Weare Our MCMC sampler uses 600 walkers to explore the parameter space constrained by 0 < S(mm/yr) < 60, 0 < d 1 (km), < 60, 0 < C(mm/yr), < 30, 0 < d 2 (km), < 40, −40 < a(mm/yr) < 40, assuming a uniform prior probability distribution over each range. An important constraint we impose is that the maximum creep depth cannot be greater than the locking depth, i.e., d 2 ≤ d 1 . Our MCMC model runs over 300,000 iterations and produces 48,000 random samples from which we estimate both the maximum a posteriori probability (MAP) solution and corresponding parameter uncertainties.
The results of our analysis are shown in Figure 8, with the observed profile velocity in red and the MAP solution in bold dashed line. The sampled marginal probability distributions for the fault slip rate, the locking depth, creep rate, and the static offset are approximately normally distributed ( Figure 9). As expected of elastic dislocation models there is a strong trade-off between the fault slip rate and the locking depth (top left box for each profile in Figure 9) where a slower slip rate can be compensated by a shallower locking depth.
The average slip rate for the whole region from the three profiles is 26 mm/yr, which is slightly faster than the GNSS-derived block model slip rate for the same region of 24.2 mm/yr [Reilinger et al., 2006]. We find that only 10% of our models for profile A-A ′ show similar slip rates to the GNSS block model constant rate to within 2 mm/yr. Sixteen percent of the models for profile B-B ′ and 28% for profile C-C ′ fall in the same range, implying that there is a systematic eastward decrease in the probability density functions for the slip rates.
To test whether the difference in MAP slip rate between profiles A-A ′ and C-C ′ is significant, we consider the null hypothesis that each of the estimated slip rates are one draw from a Gaussian distribution with the same expected value (but with different standard deviations).
If the hypothesis is true, the distribution of the difference in MAP slip rates will be Gaussian with a mean of zero and standard deviation = √ 2 A + 2 C , where 2 A and 2 C are the variance of the estimator for slip rate between profiles A-A ′ and C-C ′ , respectively. The ratio of (S A −S C )∕ √ ( 2 A + 2 C ), where S A and S C are the MAP slip rates for A-A ′ and C-C ′ , respectively, can therefore be used to test the null hypothesis. A value of 1.96 or more should only occur 5% of the time if the null hypothesis is true. In our case we find the ratio to be equal to 2.28, so we reject the null hypothesis at the 5% level meaning our results indicate that the rates are different with >95% confidence.
Our map of fault-parallel velocity (Figure 7a) shows a lateral variation in far-field velocities. For example, at 40 ∘ N the fault-parallel velocity decreases from 28-30 mm/yr on profile A-A ′ to 15-20 mm/yr on profile C-C ′ . Assuming the far field to the north is pinned to zero, as would be the case in a Eurasia-fixed reference frame, the fault-parallel velocities show an eastward decrease in relative velocity between the region north of the fault and the region to the south, which would result in decreasing fault slip rate.
The GNSS study of Yavaşoglu et al. [2011], which overlaps with the eastern edge of our fault-parallel InSAR velocities estimated a fault slip rate of 20.5 ± 1.8 mm/yr, which is consistent with our estimate of 21 mm/yr (14-27, 95% CI) for the eastern profile (C-C ′ ). In general, our estimates are comparable with the slip rate estimates from GNSS studies in this region, which range between 17 and 34 mm/yr [e.g., Oral et al., 1993;Noomen et al., 1996;Ayhan et al., 2002;Reilinger et al., 2006]. However, our rate of 30 mm/yr to the west is at the higher edge of the range of published estimates.
An important limitation of the simple dislocation models used in this study is that they assume the elastic properties of the crust do not vary along the fault, which is not always the case for faults. These differences may arise due to changes in fault zone geometry and elastic properties due to permanent damage [e.g., Perrin et al., 2016], or to specific rock geology [e.g., Ben-Zion, 2008] and the presence of fluids. Variations in crustal rheology could change the strain accumulation on the fault, which would result in different slip rates. However, the simple elastic dislocation model matches the data well and is able to give a first-order estimate of the fault slip rate and locking depth.

Fault Creep Along the Central NAF
To investigate the pattern of aseismic creep along the central NAF, we plot short profiles extending 5 km either side of the fault at regular locations (every ∼5 km) along the central NAF (Figure 10b), projecting the LOS velocities from within 2.5 km onto each profile. We fit two straight lines through the velocities on either side of the fault, taking into account the covariance and determine the offset at the fault trace, which corresponds to the LOS creep rate.
Our results (Figure 10a) clearly show that a ∼130 km section of the central NAF is undergoing aseismic creep at average rates of ∼4 mm/yr in the LOS for descending and ∼3 mm/yr for ascending. The extent of creep is in agreement with the ∼125 km estimated by Cetin et al. [2014] but larger than the ∼70-80 km estimated by Cakir et al. [2005] and Kaneko et al. [2013]. We find no fault creep above our noise level (∼1 mm/yr in the LOS) west of about 31.2 ∘ E and east of about 33.5 ∘ E. Hussain et al. [2016] showed that creep estimates can be contaminated by vertical motions. To test this we use the estimated north-south component of motion from the interpolated GNSS velocities ( Figure S5) along with the creep estimates from both ascending and descending tracks to calculate the east-west and vertical components of motion using equation 3. We calculate the fault-parallel component of the creep rate assuming the fault strikes at N81 ∘ E. Figure 10c shows our estimated fault-parallel (in red) and vertical (in blue) components of motion for the fault creep rate. There appears to be little vertical motion along the creeping segment. The maximum fault creep rate is 14 ± 2 mm/yr along a portion of the fault located ∼30 km east of Ismetpasa. The average rate for the entire creeping section is 8 ± 2 mm/yr.

Iterative Unwrapping Benefits and Limitations
Our new iterative unwrapping procedure reduces the number of unwrapping errors in the overall small baseline network and thus improves the InSAR coverage as more correctly unwrapped pixels are added to the network instead of being discarded. However, it is clear that the process cannot fix all unwrapping errors ( Figure 5). We find that there is a sharp increase in the total number of error-free pixels within the first eight iterations after which the improvements are small. Therefore, to minimize unwrapping errors from the network, some interferograms with particularly poor unwrapping still need to be removed. An efficient procedure would be to run the unwrapping process for 8-10 iterations, remove any particularly bad interferograms (therefore modifying the input to the unwrapping algorithm), and repeat the iterations.
Traditionally, interferograms with unwrapping errors have either been discarded [e.g., Pinel et al., 2011;Hussain et al., 2016] or have been fixed manually [e.g., Hamlyn et al., 2014;Pagli et al., 2014]. Manual fixing requires drawing a polygon around the unwrapping errors in every interferogram and adding or subtracting an arbitrary integer multiple of 2 until the phase sum around an interferometric loop is equal to zero. This can be a very time-consuming and labor intensive process. The strength of our procedure is that the process is automated. However, as we show in Figure 4, our procedure cannot fix all unwrapping errors and so requires some manual intervention in discarding (or correcting) particularly bad interferograms.
An important limitation using our technique is that it requires a redundant small baseline network in order to compute the phase sum around closed interferometric loops. We cannot automatically detect unwrapping errors in individual isolated interferograms.
The aim of this method is to fix pixels that are unwrapped correctly. By adding a high cost to amending the unwrapped values for these pixels, the hope is that the next iteration of unwrapping will correctly unwrap the phase of nearby pixels. The method does not address the cause of the unwrapping error, however, which in some cases cannot be overcome simply by repeating the unwrapping process. Hence, some pixels remain badly unwrapped after any number of iterations.
Another limitation is that we inherently assume an "error-free" pixel, i.e., a pixel that undergoes loop closure, is unwrapped correctly. There may be special circumstances in which this may not be the case. Consider the simplest loop consisting of three acquisitions A, B, and C with interferograms AB and BC along the forward arc and CA on the return arc. If a particular set of pixels in either one of the forward arc interferograms (AB or BC) has an unwrapping error and these exact same pixels have the same magnitude error but with the opposite sign in interferogram CA, then those pixels will still undergo loop closure and be classed as "error-free" in our technique.
However, in reality most interferograms are a part of multiple interferometric loops. And so if this error occurs in one loop and not the other, our method can still detect it; i.e., interferogram BC is part of triangular loops ABC and BEC. Our unwrapping procedure becomes more robust with greater network redundancy. However, care should be taken not to introduce interferograms with large perpendicular and/or temporal baselines as they are likely to have unwrapping errors.

Interseismic Slip Rates
Our horizontal velocity field created by combining velocities from seven InSAR tracks, in both ascending and descending geometries in a GNSS-fixed Eurasia reference frame (Figure 7) confirms the right-lateral sense of motion expected from the North Anatolian Fault. Our simple elastic dislocation models fit the fault-parallel velocities within the 95% confidence range (Figure 8) with a statistically significant decrease in fault slip rate from 30 mm/yr (25-34, 95% CI) in the east to 28 mm/yr (23-33, 95% CI) to 21 mm/yr (14-27, 95% CI). Our estimated locking depths of 13 km (6-20, 95% CI), 13 km (5-22, 95% CI), and 17 km (10-25, 95% CI) show no such pattern. Our statistical test to discard the hypothesis of a constant slip rate assumes that the uncertainty attributed to the data is correct. If the uncertainty were underestimated due to the possibility that the apparent change in slip rates could result from other physical mechanisms such as other deformations or change in crust rheology, the level of confidence could be overestimated [e.g., Duputel et al., 2014].
The positive trade-off between the fault slip and locking depths means that a decreasing fault slip can be compensated by a decreasing locking depth near the fault. This would explain the large confidence intervals for these parameters and could explain the lateral variation in these parameters. However, if we assume the velocities in the far field to the north are zero, as we would expect with velocities in a Eurasia-fixed reference frame, then the far-field plate velocities (velocities to the far south on each profile) do appear to be decreasing eastward along the fault, from ∼30 mm/yr in profile A-A ′ to ∼20 mm/yr in profile C-C ′ (Figure 11), implying that the lateral change in these parameters are real variations along the fault. This pattern is also observed in the GNSS velocities (Figure 8).
There is a relatively sharp change in fault-parallel velocity south of the NAF (Figure 7) that coincides with the B-B ′ profile line. The feature does not correspond to a track boundary (Figure 1). Figure 12 shows the fault-parallel velocities projected onto profile D-D' that show this gradient between 100 km and 140 km. It is clear that the variation along the profile broadly matches the GNSS velocities, although the gradient at 120 km is steeper in the InSAR than the GNSS. This might be due to local atmospheric residuals in the InSAR velocities. The gradient does not correspond to any topographic changes along the profile. Ergintav et al. [2009] showed that the 1999 earthquakes resulted in postseismic deformation as far as Ankara, which is less than 100 km south of the NAF in this region. Therefore, the faster velocities to the west of the study region could be due to postseismic deformation from the 1999 earthquakes with the sharp gradient representing the eastern limit of postseismic deformation. Figure 11. Fault-parallel velocities for each profile shown in Figure 8 with the velocities in pale blue, pale red, and pale green corresponding to profiles A-A ′ , B-B ′ , and C-C ′ , respectively. Our best fit (MAP solution) model is shown by the bold line through the velocities. It is clear that there is a far-field decrease in velocity from profile A-A ′ to profile C-C ′ .
The largest recent earthquakes on the central portion of the NAF in recent times were the 1943 Tosya (M w 7.7), the 1944 Bolu-Gerede (M w 7.5), and the 1951 Kursunlu (M w 6.9) earthquakes ( Figure 13). Our fastest slip rate of 30 mm/yr corresponds to the peak coseismic slip region of the 1944 earthquake, while the central profile with 28 mm/yr corresponds to the 1951 earthquake slip, and the easternmost profile with the slowest slip rate of 21 mm/yr covers the 1943 earthquake rupture. In the case of the two largest earthquakes the coseismic surface slip decreases to the east. Previous studies have shown that overall coseismic slip decrease is indicative of off-fault strain dissipation [e.g., Manighetti et al., 2005]. If this pattern of off-fault strain dissipation also occurs during the interseismic period, then our model, which assumes all the slip occurs on the fault, would overestimate the slip rate on the fault. However, it remains unclear if distributed off-fault fault deformation occurs during the interseismic period. A dense network of long-term continuous GNSS measurements around the fault would help determine if this is an important mechanism of long-term strain dissipation.
Given the 95% confidence intervals, there is no significant statistical difference in the MAP slip rates for profiles A-A ′ and B-B ′ . These profiles also have the same MAP locking depth (13 km). Whereas the MAP slip rate and locking depth for profile C-C ′ , which crosses the 1944 earthquake rupture, are significantly different to those of the profiles over the 1943 earthquake. Similarly, the velocity change observed south of the fault (profile D-D' in Figure 12) roughly coincides with the limit between the two broken segments in the earthquakes. It is therefore possible that this difference arises due to large-scale fault segmentation coinciding with the boundary between the two large earthquakes [e.g., Manighetti et al., 2015;Perrin et al., 2016].
The change in slip rate along the fault could also arise from east-west extension within Anatolia. Earthquake moment tensors show significant number of earthquakes within Anatolia (Figure 7b), several with normal faulting mechanisms, implying that there is ongoing internal deformation within Anatolia. Aktug et al. [2013] also found significant ongoing deformation within Anatolia from detailed analysis of GNSS velocities in central Anatolia, which were more consistent with east-west elastic elongation rather than a rigid-body rotation [Reilinger et al., 1997;McClusky et al., 2000] or simple transport [Reilinger et al., 2006].
The average fault slip rate across the central NAF from our three profiles is 26 mm/yr, which is similar to the slip rate determined using GNSS alone for the region [e.g., Reilinger et al., 2006;Nocquet, 2012].

Fault Creep
Our estimates of fault creep rate by direct offset measurements of LOS velocity across the fault reveal that a ∼130 km portion of the central NAF is undergoing aseismic creep that reaches the ground surface.
Over the InSAR time interval, the fault creep rate has a maximum of 14 ± 2 mm/yr around 30 km east of Ismetpasa, which is slightly slower than the value determined by Cetin et al. [2014], who found the maximum creep to be 20 ± 2 mm/yr at the same location. This discrepancy can be explained by the fact that they used LOS velocities from a single look direction (descending). Using our descending velocities alone, which is the same data set used by Cetin et al. [2014], we estimate a similar maximum fault creep rate of 21 ± 2 mm/yr. This study is a confirmation that where available, both ascending and descending information can be used to estimate accurate and unbiased values of creep or other surface deformation that is not contaminated by vertical motions Our average creep rate for the entire portion of the creeping sections is 8 ± 2 mm/yr. This is similar to our MAP solution from our elastic model for profile B-B ′ (10 mm/yr). Our estimate for the average fault creep rate is similar to recent estimates by Karabacak et al. [2011], Ozener et al. [2013], Kaneko et al. [2013], and Cetin et al. [2014] who estimate average creep rates of 6-9 mm/yr, 7.6±1, 9 mm/yr, and 8±2 mm/yr, respectively. Our MAP solution for the depth extent of aseismic fault creep (9 km) is deeper than the 5 km estimated by Cetin et al. [2014] and the 4 km estimated by Rousset et al. [2016]. However, our 95% confidence bound on this parameter is large (1-20 km). It is possible that we are biased toward deeper depths because we resample our velocities to a 1 km by 1 km grid, which could be insensitive to very shallow creep depths. However, Hussain et al. [2016] showed that changing the creep depths over a large range (4 km to 12 km) only results in a small difference in the shape of the profile close to the fault, which is below the estimated uncertainty in the fault-parallel velocities. Therefore, it is more likely that the large confidence bound on the creep depth extent is due to the noise in the data. Bilham et al. [2016] used creepmeter measurements across the Ismetpasa section of the NAF to show that interannual surface slip is episodic and consists of periods of no slip (47% of the time in the past 2 years), interrupted by months of slow slip (44% of the time in the past 2 years) at rates of about 3 mm/yr or by abrupt slip events with transient velocities exceeding 3 mm/h with slip durations of many days, and, in the case of multiple events, with cumulative amplitudes of many millimeters. They determined near-fault average creep rate of 6.1 mm/yr with creep events extending down to depths of 3-7 km. The creep rate estimates are slightly lower than our estimate of 8±2 mm/yr, but this may be due to the creepmeter's incompletely sampling the full width of the surface shear zone. As discussed above, the locking depth determined by the creepmeter study is comparable to previous studies with our estimate of 9 km toward the upper bound of these estimates. Figure 13 shows that the majority of the creeping section is located on the eastern section of the 1944 M w 7.5 earthquake with creep mostly occurring where coseismic slip was lower. The first measurements of aseismic creep along this section of the fault were made by Ambraseys [1970], who estimated a creep rate of ∼20 mm/yr near the town of Ismetpasa. Although it is not known whether the fault was creeping before the 1944 earthquake, numerous studies have shown that the surface creep rate follows an exponential decay through time to a current steady state value of ∼8 mm/yr [e.g., Cakir et al., 2005;Kutoglu et al., 2010;Kaneko et al., 2013;Cetin et al., 2014], implying that aseismic creep was initiated as postseismic deformation following the large earthquake. Cetin et al. [2014] also showed that aseismic surface creep can, to some extent, be correlated with the geology along the North Anatolian Fault. The majority of the creeping segment is correlated with an Upper Jurassic-Lower Cretaceous limestone unit and could have been initiated due to pressure solution.
The average creep rate is about a third of the average fault slip rate (26 mm/yr) for this portion of the NAF, implying that strain is still accumulating along the fault. Shallow aseismic creep reduces the rate of interseismic strain accumulation by 30-40% compared to if the fault was fully locked. However, fault creep can increase the stresses at the edges of the creeping zone and thus bring the adjacent fault segments closer to failure. Assuming a uniform steady state creep rate of 8 ± 2 mm/yr down to 6 ± 3 km depth (average of Cetin et al. [2014], Rousset et al. [2016], and our MAP solution) along the entire 130 km creeping segment of the fault and 26 mm/yr (21-32, 95% CI) down to a locking depth of 14 km (7-22, 95% CI), in 200 years (approximate earthquake repeat time [Stein et al., 1997]) the creeping segment of the fault will have accumulated strain equivalent to an earthquake with moment magnitude between 7.4 and 8. This large range is mostly due to the large confidence range for our model parameters. Using the average MAP solution from the three profiles gives a strain deficit equivalent to a moment magnitude 7.7 earthquake in a 200 year period.

Conclusion
We have presented a new iterative unwrapping technique for small baseline InSAR processing that can be used to iteratively identify and mitigate unwrapping errors, therefore increasing the number of correctly unwrapped pixels in the small baseline network and improving the InSAR coverage compared to methods where unwrapping errors are rejected or masked. We have used this technique to process Envisat SAR data from seven tracks in both ascending and descending geometries spanning the time window between 2003 and 2010. The footprint of our tracks cover the entire central portion of the North Anatolian Fault in both viewing geometries. We combine the InSAR LOS velocities with published GNSS to create a horizontal velocity field for the region (assuming negligible vertical motions). Profiles through the fault-parallel velocities reveal an eastward decreasing fault slip rate (30 mm/yr, 28 mm/yr, and 21 mm/yr) with no such pattern in the locking depths (13 km, 13 km, and 17 km). Direct offset measurements of LOS velocity across the fault reveal that a ∼130 km portion of the central NAF is undergoing aseismic fault creep that reaches the ground surface at an average rate of 8 ± 2 mm/yr. The maximum creep rate of 14 ± 2 mm/yr is slower than previous estimates, which were biased by using data from only a single satellite look direction. We conclude that shallow aseismic creep on the central section of the NAF reduces the rate of interseismic strain accumulation by 30-40% compared to if it was fully locked. Nevertheless, the fault is still accumulating strain and remains capable of producing a large earthquake in the future.