Geodetic observations of postseismic creep in the decade after the 1999 Izmit earthquake, Turkey: Implications for a shallow slip deficit

The relationship between aseismic slip and tectonic loading is important for understanding both the pattern of strain accumulation along a fault and its ability to generate large earthquakes. We investigate the spatial distribution and temporal evolution of aseismic creep on the western North Anatolian Fault (NAF) using time series analysis of Envisat interferometric synthetic aperture radar (InSAR) data, covering the full extent of the 1999 Izmit and Düzce earthquake ruptures and spanning 2002–2010. Discontinuities in the line‐of‐sight velocity across the fault imply that fault creep reaches the Earth's surface at an average fault‐parallel rate of ∼5 mm/yr along an ∼80 km section of the NAF. By combining InSAR and published GPS velocities, we are able to extract the vertical and east‐west components of motion and show that the Adapazari basin is subsiding at a rate of ∼6 mm/yr. Vertical motions have biased previous estimates of creep in this region. The displacement time series close to the fault is consistent with an afterslip model based on rate‐and‐state friction, which predicts a rapid deceleration in fault creep rate after the Izmit earthquake to a near‐steady state ∼5 mm/yr after 5 years. Projecting our model 200 years into the future we find that the cumulative displacement of 1–1.3 m is insufficient to account for the shallow coseismic slip deficit observed in previous studies. Distributed off‐fault deformation in the shallow crust or transient episodes of faster slip are likely required to release some of the long‐term strain during the earthquake cycle.


Introduction
While the upper crustal portions of most active faults are locked, storing elastic strain energy for decades to centuries and releasing it almost instantaneously in earthquakes, some faults instead slip steadily at slow rates over various time scales and spatial distributions. These steadily slipping (creeping) faults may store little to no strain energy and are unlikely to produce significant earthquakes if aseismic creep occurs throughout the seismogenic crust and the creep rate is equal to the tectonic loading [Bürgmann et al., 2000].
However, most faults do not undergo aseismic creep at all depths in the crust at the full tectonic loading rate. Examples from the Hayward Fault [Schmidt et al., 2005], the Longitudinal Valley Fault [Champenois et al., 2012;Thomas et al., 2014], the central section of the San Andreas Fault [Maurer and Johnson, 2014;Jolivet et al., 2015], and the Ismetpasa section of the North Anatolian Fault [Kaneko et al., 2013;Cetin et al., 2014] show that aseismic fault slip occurs in the upper crust on some segments of major strike-slip faults at rates that are significantly less than the full tectonic loading rate, implying that not all the accumulated strain energy is being released aseismically. In these cases, parts of the fault are fully or partially locked and can still generate moderate to large earthquakes [e.g., Avouac, 2015]. The spatial and temporal distributions of fault creep rate are therefore important for understanding the pattern of strain accumulation along a fault and its ability to generate large, damaging earthquakes.
The North Anatolian Fault (NAF) is a major continental right-lateral strike-slip fault located in northern Turkey. Together with the East Anatolian Fault, it facilitates the motion of Anatolia away from the Arabia-Asia collision zone toward the Hellenic subduction zone. Since the 1939 M w 7.9 Erzincan earthquake in eastern Turkey, the NAF has slipped in a sequence of large (M w > 6.7) earthquakes with a dominant westward progression in seismicity [Barka, 1996;Stein et al., 1997]. This sequence of earthquakes has been interpreted as a result of stress transfer along strike, where one earthquake brings the adjacent segment closer to failure [Stein et al., 1997;Hubert-Ferrari et al., 2000]. The most recent events were the M w 7.4 Izmit and M w 7.2 Düzce earthquakes in 1999. The Izmit earthquake ruptured ∼140 km of the western section of the North Anatolian Fault on 17 August 1999 [e.g., Wright et al., 2001b;Barka et al., 2002] and was followed by the Düzce earthquake on 11 November 1999, which ruptured another ∼45 km of the fault [e.g., Akyuz, 2002;Bürgmann et al., 2002a] ( Figure 1). Cakir et al. [2012] were the first to document postseismic fault creep along the Izmit rupture. Using interferometric synthetic aperture radar (InSAR) analysis of Envisat satellite images from ascending tracks 157 and 386 between 2003 and 2009, they observed a discontinuity in the InSAR velocities across the section of the rupture between Izmit and Akyazi. Their analysis of ERS satellite images showed no evidence of creep prior to the 1999 earthquakes. The authors concluded that the observed fault creep is postseismic deformation initiated by the 1999 Izmit earthquake.
Using elastic dislocation models, Cakir et al. [2012] estimate the aseismic creep rate to reach a maximum of 27 mm/yr and to extend from the surface to a depth of 12 km. This estimated creep rate is comparable to published geodetic slip rates for the NAF at this longitude, which lie between 11 and 26 mm/yr [Straub et al., 1997;Hubert-Ferrari et al., 2000;Ayhan et al., 2002;Meade et al., 2002;Le Pichon et al., 2003;Reilinger et al., 2006;Aktug et al., 2009;Ergintav et al., 2009Ergintav et al., , 2014. If accurate, their estimate suggests that little or no elastic strain is currently accumulating along this section of the NAF.
In this study we measure surface velocities between 2002 and 2010 across the region encompassing the Izmit and Düzce ruptures. We use Envisat advanced synthetic aperture radar (ASAR) images from three descending and two ascending tracks along with published GPS velocities. Each satellite track roughly covers a 100 km by 400 km area. We use this velocity field to investigate the spatial distribution and temporal evolution of aseismic creep along the Izmit rupture. We use an elastic half-space dislocation model to determine the depth and rate of the creep and examine its temporal behavior using an afterslip model based on rate-and-state friction [Ruina, 1983;Rice and Ruina, 1983;Tse and Rice, 1986;Rice, 1993;Segall, 2010].

InSAR Processing and Applied Corrections
Our data set consists of 96 Envisat images from five overlapping tracks acquired in descending (64, 293, and 21) and ascending (157 and 386) geometries ( Figure 1b). The images span the period between 2002 and 2010 and fully cover the Izmit and Düzce ruptures. Details of the data processed for each track are given in Table 1.
We focus the raw synthetic aperture radar (SAR) image products using the Jet Propulsion Laboratory/Caltech ROI_PAC software [Rosen et al., 2004] and constructed 229 interferograms using the DORIS software HUSSAIN ET AL.
IZMIT CREEP 2981  [Kampes et al., 2004]. The interferograms were chosen to minimize the time difference between acquisition dates (the temporal baseline) and the spatial separation of the satellite orbits (the perpendicular baseline). We correct topographic contributions to the radar phase using the 3 arc sec Shuttle Radar Topography Mission digital elevation model [Farr et al., 2007] and account for the known oscillator drift for Envisat according to Marinkovic and Larsen [2013].
We remove incoherent pixels and reduce the noise contribution to the deformation signal by applying the Stanford Method for Persistent Scatterers (StaMPS) persistent scatterer small baseline time series InSAR technique [Hooper, 2008;Hooper et al., 2012], which takes advantage of the spatial correlation between pixels and does not impose a temporal deformation model when identifying targets with stable phase characteristics through time. The StaMPS software selects only those pixels that have stable phase noise characteristics in time and uses this subset to compute velocities and time series.
The small baseline network allows for unwrapping error checks by summing the phase around closed interferometric loops [Biggs et al., 2007]. Interferograms showing obvious unwrapping errors were corrected manually; any others that could not be corrected were removed. In this way we ensure that we have a redundant network of interferograms with minimal unwrapping errors, which enables us to make a more robust estimate of the time-averaged line-of-sight (LOS) velocity. Figure 2 shows that we are left with a good redundant network spanning the time series for each track. Over the five tracks we use a total of 127 interferograms in the final redundant small baseline networks. The uncertainties on the final velocity for each pixel are calculated using bootstrap resampling and are presented at the 1 sigma level in the following work.
As the InSAR phase delay is a superposition of multiple signals, including tectonic deformation, atmosphere, and orbital errors, additional corrections are required. In section 2.1 we elaborate on our atmospheric corrections. As InSAR is a relative measurement, we simultaneously account for orbital errors and any remaining long-wavelength signals by combining the InSAR velocities with published GPS velocities in a Eurasia-fixed reference frame (section 2.2).

Atmospheric Delay Corrections
The spatial and temporal variations in tropospheric humidity, pressure, and temperature are 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].
We calculate the estimated phase delay due to the atmosphere for each of our interferograms using the ERA-Interim global atmospheric model reanalysis product [Dee et al., 2011] obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF). The ERA-Interim product provides atmospheric information at approximately 80 km spatial resolution on 60 vertical levels from the surface up to 0.1 hPa every 12 h. We then remove this atmospheric signal from each interferogram, implementing the correction using the TRAIN (Toolbox for Reducing Atmospheric InSAR Noise) software package [Bekaert et al., 2015c] after the method of Doin et al. [2009] andJolivet et al. [2011].
On average over the five InSAR tracks, the tropospheric correction reduces the standard deviation of each interferogram by about 6%. Note that as the standard deviations are mainly reflecting the capability to correct for the long-wavelength tropospheric signal, we removed a ramp from each interferogram prior to the standard deviation computation. The average reduction in standard deviation is small, implying that the ERA-Interim weather model is not capturing the full tropospheric variation, which means that some residual atmospheric signals still remain after the ERA-Interim correction ( Figure S1 in the supporting information). The average reductions in standard deviation per track are 13% for track 64, 3% for track 293, 8% for track 21, 7% for track 157, and −2% for track 386 ( Figure S2). Note that the weather model correction makes the interferograms for track 386 slightly noisier; this is also the track with the least number of interferograms in the final small baseline network. In total, the standard deviation is reduced in 62% of our interferograms after the ERA-I atmospheric correction. Figure 3 shows the calculated average line-of-sight (LOS) velocity through the InSAR time series with the reference for each track shown by the orange star. For all tracks, blue is motion toward the satellite and red is away. Our results are consistent with a right-lateral sense of motion across the fault. The difference between track 293 and track 21 in the overlap region is likely due to residual atmosphere. The higher uncertainties in the overlap region ( Figure S4) reflect this discrepancy.

InSAR Line-of-Sight Velocity Field in a GPS Reference Frame
To obtain a consistent velocity field across the region, we transform our InSAR velocities for each track into a Eurasia-fixed reference frame as defined by the Global Strain Rate Model project [Kreemer et al., 2014], from which we download the compilation of input GPS data. The GPS velocities immediately around the Izmit rupture are those published by Reilinger et al. [2006], which were derived from pre-1999 earthquake observations and therefore do not include postseismic or coseismic deformation [Reilinger et al., 2000;Ergintav et al., 2002].
We transform the InSAR into the GPS-Eurasia reference frame by first averaging the InSAR velocities that fall in a 1 km radius around every GPS station within the boundaries of the InSAR track. We project the GPS velocities into the local satellite line of sight and calculate the difference from the InSAR velocities. We then determine the best fit plane through the residual velocities using a weighted linear least squares adjustment. We remove this plane from the InSAR velocities to transform the LOS velocities into a Eurasia-fixed GPS reference frame.

Along-Strike Variation in Fault Creep Rate
We observe a discontinuity in the line-of-sight velocities across the fault, seen as a sharp color contrast in Figure 4b. This is most clearly seen in tracks 64, 293, and 157 between 29.9 ∘ E and 30.7 ∘ E. This discontinuity is superimposed onto a longer wavelength smooth variation in velocity across the fault. Figure 4a shows HUSSAIN ET AL. IZMIT CREEP 2983  profiles of line-of-sight velocities through tracks 293 and 157. These two components of the deformation are most clearly seen in track 293 where a velocity discontinuity at the fault location is imposed onto a longer wavelength variation across the fault.
We interpret the long-wavelength signal to result from the relative motion of Anatolia with respect to Eurasia in the lower crust and upper mantle. The velocity discontinuity is due to shallow fault creep.
We calculate the LOS variance-covariance matrix of the noise for each track by computing the average radial covariance versus distance (autocorrelation) using the velocities in a 50 km by 50 km region ∼200 km to the south of the fault. This region is assumed to have no tectonic deformation and contain only atmospheric noise. We fit a covariance function, C(r), of the form estimating the variance, 2 , and the characteristic length , which gives the spatial correlation of noise as a function of distance between pixels, r. Our best fit values for each track, and the east-west velocities used in section 5, as well as the center of the region used to calculate the covariance function are shown in Table 2.
We estimate the rate of fault creep at various locations along the fault trace for both the Izmit and Düzce ruptures using the LOS velocities from each track. We first make short profiles of LOS velocity extending 5 km either side of the fault and then project velocities onto this line from within a 2.5 km window either side of HUSSAIN ET AL. IZMIT CREEP 2985 the profile. We fit two straight lines through the LOS velocities on either side of the fault and determine the offset at the fault trace (v creep ), i.e., the LOS fault creep rate, using the following linear set of equations: where v N (1∶i) and v S (j∶n) are the LOS velocities north and south of the fault, respectively; x N (1∶i) and x S (j∶n) are the perpendicular distance north and south of the fault, respectively; a N and b N are, respectively, the gradient and offset of the best fit line through the velocities north of the fault, a S the gradient of the best fit line through the southern velocities, and represents errors in the model.
We solve these equations and determine the error distribution of each parameter using the percentile bootstrap method [Efron and Tibshirani, 1986]. Unlike the regular bootstrap algorithm, where n random samples are taken from n observation, we select n 1 and n 2 random samples from n 1 and n 2 observations (LOS velocities) north and south of the fault, respectively [Bekaert et al., 2015b], where n 1 + n 2 is the total number of observations on both sides of the fault. Each bootstrap simulation provides an estimate for the unknown parameters (a N , a S , b N , and v creep ). We do this calculation at each location marked with a black circle in Figure 5b for every track.
In most cases the Best Linear Unbiased Estimator (BLUE) value for the fault creep approximates the mean from bootstrap resampling. We use the bootstrap results for our creep estimates because in some cases where the errors on each point are particularly large, the BLUE technique underestimates the uncertainty on the fault creep rate (see Figure S3).
Our results (Figure 5a) show that a section of the Izmit rupture extending about 80 km, from the Gulf of Izmit in the west to as far as 30.7 ∘ east, has undergone shallow creep during the period 2002-2010, at an average LOS rate of ∼2.3 mm/yr and ∼3.3 mm/yr in the descending and ascending tracks, respectively. The western extent of the creep is unknown due to the lack of near-fault geodetic data in the gulf.
The LOS fault creep rates east of 30.7 ∘ E (beyond 100 km in Figure 5a), covering the eastern end of the Izmit rupture and all of the Düzce rupture, appear to show little to no resolvable creep in the LOS descending geometry. The ascending fault creep rates (from track 386) in this region have large errors due to a sparse pixel coverage resulting from InSAR decorrelation, particularly in the mountains south of the fault.
The difference in sign between the creep rates on ascending and descending tracks ( Figure 5a) is a result of the two different satellite-viewing geometries. Vertical motion manifests as ascending and descending InSAR signals with approximately equal magnitude and the same sign in the LOS while east-west motions, by contrast, result in signals with opposite sign. Therefore, signals that show a positive correlation between ascending and descending tracks (e.g., between 20 km and 100 km distance in Figure 5a) are indicative of vertical motion.
The Izmit rupture is oriented east-west with little north-south coseismic deformation in the region adjacent to the fault; the fault plane solution for the earthquake (Figure 1) shows pure east-west right-lateral strike-slip HUSSAIN ET AL.
IZMIT CREEP 2986 motion, and GPS velocities adjacent to the fault are also approximately east-west ( Figure 1). To simplify our calculations, we therefore assume that deformation in the north-south direction is negligible. This allows us to decompose the displacement rates at points where we have both ascending and descending information into east-west (approximately fault-parallel) and vertical rates (Figure 5c).
West of 30.7 ∘ E, we find an average horizontal fault-parallel creep rate of ∼5 mm/yr with a maximum horizontal creep rate of 11 ± 2 mm/yr near the city of Izmit. This maximum rate is significantly slower than the estimate of 27 mm/yr between 2003 and 2009 from Cakir et al. [2012] but is still more than a third of the long-term slip rate on the fault. We find this to be due to contamination of the estimated horizontal velocities by vertical motions, which were assumed negligible by Cakir et al. [2012] on which we further elaborate in section 4.1.
In general, the fault creep rate decreases along strike toward the east with a small increase near the town of Akyazi. East of 30.7 ∘ E, including the Düzce rupture, the rate can be considered to be zero within uncertainty.

InSAR Velocity Decomposition
To further investigate the spatial distribution of the apparent vertical motions, we decompose our full InSAR velocity field into east-west and vertical components (Figures 6a and 6b). We do this first by resampling our InSAR LOS velocities 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 2 km of the center of each grid point, and we reference each track to a Eurasia-fixed GPS reference frame, as described in section 2.2. For every pixel where information from both ascending and descending HUSSAIN ET AL. IZMIT CREEP 2987 geometries are available, we invert for the east-west and vertical components of the velocity using the method described by Wright et al. [2004], taking into account the local incidence angles and assuming that there is no north-south motion.
The east-west component of the InSAR velocity ( Figure 6a) clearly shows fault creep as a velocity discontinuity on the Izmit earthquake rupture continuing east along the fault from the Gulf of Izmit. The velocity discontinuity becomes unclear east of about 30.7 ∘ E, implying no shallow creep reaches the Earth surface east of this longitude. The vertical component (Figure 6b) shows no clear discontinuities across the fault other than in the Adapazari basin region.

Adapazari Basin Subsidence
The Adapazari basin is roughly 30 km wide and located around the town of Adapazari to the north of the NAF. It is located in a transtensional region between the Izmit-Adapazari and Mudurnu valley segments of the NAF [Emre et al., 1998;Ünay et al., 2001]. After detailed analysis of the regional geology, fault pattern, morphology, and stress and strain pattern in the region, Neugebauer [1995] suggested a pull-apart mechanism for the basin formation. The superficial geology consists of alluvial and fluvial sediments deposited by the Sakarya River, which flows though the center of the basin [Yigitbas et al., 2004;Ulutaş et al., 2011]. The thickness of the alluvium in the midsection of the basin is estimated to be more than 200 m according to groundwater wells drilled by the State Hydraulic Works [1983;Ulutaş et al., 2011] and possibly up to 2 km thick determined from a seismic refraction survey across the basin [Karahan et al., 2001].
HUSSAIN ET AL. IZMIT CREEP 2988 The displacement in the Adapazari basin (Figures 6c-6e) is seen only in the vertical component, with the basin subsiding at a rate up to 6 mm/yr. The subsidence is bound to the north and west by a region of high topography ( Figure 6e) and by the NAF to the south. The eastern limit of subsidence correlates with the extent of Holocene sedimentary deposits [Yigitbas et al., 2004].
The cause of this subsidence is unclear, but there are several possible explanations. The Adapazari basin is a region of high agricultural productivity [Erinç and Tunçdilek, 1952;Gedikli, 2004;Ikiel et al., 2012]. The two regions of rapid subsidence (the south-east corner and a smaller region in the north-west corner of the basin) are former wetland/swamps, which were drained in the last 30 years [Bilgin, 1984]. The subsidence could, therefore, be a surface response to this drainage. To first order, the shape of the subsiding region from our InSAR velocities matches the shape of the basin (Figures 6d and 6e); hence, the subsidence could also be due to sediment compaction. However, there may also be a tectonic influence on the subsidence in this region; Neugebauer et al. [1997] and Poyraz et al. [2015] showed that microseismicity fault plane solutions within the basin have significant normal faulting components. Solutions for the 1967 M w 7 earthquake (strike-slip) and its main aftershock (M w 5.5, normal) [Jackson and McKenzie, 1984] are also in good agreement with a pull-apart mechanism (Figure 6a) for this basin formation. Cakir et al.'s [2012] estimate of the fault creep rate was derived using a single ascending satellite orbit, in which eastward motion and subsidence both cause motion away from the satellite. Their estimate of the fault creep rate is therefore contaminated by subsidence from these vertical signals.

Modeling Profile Velocities
A fault-perpendicular profile of east-west velocities (Figure 7) confirms that there are two different deformation signals: a long-wavelength signal, related to relative interseismic motion of Eurasia and Anatolia at depth, and a velocity discontinuity at the fault, consistent with shallow fault creep. We model the long-wavelength signal as a screw dislocation at a depth d 1 in an elastic half-space, equivalent to interseismic slip on a fault plane below locking depth d 1 at rate S [Weertman and Weertman, 1964;Savage and Burford, 1973]. For the shallow creep, between the surface and depth d 2 at rate C (Figure 8), we use a back slip approach [Savage, 1983] in which the shallow creep is modeled as the sum of slip on the entire fault plane (Heaviside function, (x)) plus a screw dislocation in the opposite sense to the plate motion at depth d 2 . We also solve for a possible static offset between the profile and the model. The fault parallel (east-west) velocity, v EW (x), is as follows: where x is the perpendicular distance from the fault.
We find best fit values for each model parameter (S, d 1 , C, and d 2 ) and offset a, using a Bayesian approach, implementing the Goodman and Weare distribution. The most important benefits of the algorithm from Goodman and Weare [2010] compared to the typically used Metropolis-Hasting algorithm is reduced convergence time and the ability to efficiently explore highly irregular probability distributions. The algorithm consists of running multiple parallel Markov chains, or walkers, where the next iteration for each chain requires randomly selecting another chain from the ensemble and choosing a new position that is a random linear combination of the positions of both walkers. Each individual chain is initiated from a randomly selected point in the parameter space within defined prior constraints.
An important prior constraint in our inversion is that the maximum depth of fault creep, d 2 , must be less than or equal to the interseismic locking depth, d 1 . In addition, our MCMC sampler explores the parameter space constrained by −60 < S < 0, 0 < d 1 < 40, −30 < a < 30, 0 < C < 15, and 0 < d 2 < 40, assuming a uniform prior probability distribution over each range.
Our initial runs assuming a uniform prior probability distribution over the locking depth range (0-40 km) revealed that the data do not constrain the maximum locking depth d 1 . We therefore include an extra prior constraint on the locking depth noting that published values are generally in the range 11-22 km [e.g., Reilinger et al., 2000;Wright et al., 2001b;Michel and Avouac, 2002], we assume a Gaussian prior constraint for d 1 = 17 km with a 1 value of 5 km.
For our model runs, we use 600 walkers to explore the parameter space over 1,000,000 iterations producing 160,800 independent random samples from which we estimate both the maximum a posteriori (MAP) solution and corresponding parameter uncertainties.
We perform the inversion on three profiles for which velocities from within 15 km perpendicular distance are projected onto the profile lines shown in Figure 6a. The results of our MCMC analysis are shown in Figure 9. For each profile the blue line is the maximum a posteriori probability solution with the 68% and 95% confidence bounds on the model represented by the light and dark grey bands. The blue line (our preferred solution) is the MAP solution while imposing the Gaussian prior constraint on the locking depth discussed above.
HUSSAIN ET AL. IZMIT CREEP 2990 The black line represents the MAP solution assuming a uniform prior constraint on the locking depth. The MAP solutions are consistent across all three profiles with the fault slip rate, locking depth, and creep depth for each profile within ±1 of the average.
Profile C-C' shows a lower creep rate than profile B-B' (2 mm/yr compared to 5 mm/yr), which is consistent with the creep rates inferred directly from the InSAR velocity discontinuities ( Figure 5). Although the MAP solutions for maximum fault creep depth are consistent between profiles (6 km, 8 km, and 7 km), the 95% confidence bounds are very large. Figure 10 helps explain the cause behind the large confidence intervals. Changing maximum creep depth between 4 km, 8 km, and 12 km, with the other parameters fixed, makes only a small difference to the shape of the profile relative to the estimated uncertainty of the east-west velocities.
HUSSAIN ET AL. IZMIT CREEP 2991 Figure 10. Profile B-B' from Figure 9 with the bold black line representing the MAP solution with a maximum creep depth of 8 km. The dark and light grey shading around the best fit line represents 2 and 1 model uncertainties, respectively. We fix the slip rate, locking depth, offset, and creep rate to the MAP solution for this profile and plot two other model calculations for creep depths of 4 km and 12 km, shown in red and green, respectively. The difference between each model is within the noise level of the velocities and is the principle cause behind the large confidence intervals for the MAP creep depths in Figure 9. Figure 11 shows the marginal probability distribution for each model parameter for profile line B-B' , which shows the clearest velocity discontinuity at the fault. The distributions for profiles A-A' and C-C' can be found in the supporting information ( Figure S5), but the following description applies to all three cases.
The sampled distributions for the fault slip rate, the static offset, and the fault creep rate are approximately normally distributed. Our distributions reveal that several parameters trade-off with one another. The clearest example is between the fault slip rate and the locking depth.

Fault Creep Time Series
Theoretical studies and laboratory experiments, in the framework of rate-and-state friction, suggest that the mechanism of fault creep is linked to steady state velocity-strengthening behavior [e.g., Ruina, 1983;Rice and Ruina, 1983]. Rate-and-state models predict that these regions slip stably under tectonic loading, whereas velocity-weakening regions produce stick-slip motion [e.g., Tse and Rice, 1986;Rice, 1993;Kaneko et al., 2013].
The relative motion of two points, ∼3 km either side of the fault west of Lake Sapanca, is constrained by a long time series of displacements from GPS [Cakir et al., 2012], as well as by InSAR displacements from our tracks ( Figure 12). We correct these for a static coseismic offset by aligning the projection of the preseismic and postseismic displacements to zero at the time of the earthquake. We correct for the small interseismic component due to slip on a vertical dislocation using the model of Savage and Burford's [1973] with a locking depth of 9 km and slip rate of 17 mm/yr (average values from our profiles in Figure 9). Any remaining relative displacement between the two stations can be attributed to shallow fault creep.
After this correction, we find that the preseismic fault creep rate (v p ) from the GPS observations is 1 ± 1 mm/yr. The time series shows a period of rapid displacement after the earthquake, which we interpret as afterslip, followed by a slow decay spanning our InSAR data time window.
We model the time evolution of fault creep using a rate-and-state afterslip formulation [Marone et al., 1991;Scholz, 2002;Segall, 2010], which approximates the behavior of the system using a simple spring slider model. We note that the characteristic decay time t c of transient afterslip depends on both frictional properties as well as the system stiffness k: where is the normal stress, a and b are rate-and-state friction parameters, and v p is the preseismic fault creep rate.
Integrating equation (4) with a substitution (for details, see Segall [2010]) gives the displacement where v max is the depth-averaged maximum afterslip velocity.
We use an MCMC approach (as described in section 5) to find the model parameters that best fit the data (Figure 12c). Our MAP solution gives a maximum afterslip velocity (v max ) of 1050±150 mm/yr or 2.9±0.4 mm/d and (a − b)∕k = 51 ± 2 mm. The calculated marginal probability distributions for each parameter (Figure 13) show that the sample distributions are approximately normally distributed.

Fault Creep and Elastic Modeling
Our InSAR results show that fault creep reaches the surface along ∼80 km of the Izmit rupture, from the Gulf of Izmit to the eastern edge of the Adapazari basin. High creep rates on the fault where it enters the Gulf ( Figure 5) suggest that aseismic creep may continue into the Gulf, extending the length of the creeping segment of the Izmit rupture. Over the InSAR time interval, our time-averaged fault creep rate has a maximum value of 11 ± 2 mm/yr (between 2002 and 2010) near the city of Izmit.
Our simple elastic model fits the data on each profile (Figure 9) within the 95% confidence bound, with a maximum a posteriori probability (MAP) estimate of the fault creep rate of 5 mm/yr occurring from the surface to a depth of 8 km. Note that 5 mm/yr is the average fault creep rate along the width of the profile (30 km), while 11 mm/yr is the maximum local creep rate near the city of Izmit. This is less than half that estimated by Cakir et al. [2012] (27 mm/yr) using velocities from ascending tracks alone. Much of this difference, particularly in the Adapazari basin, is likely from contamination of the horizontal displacement rates by vertical motion, which was neglected by Cakir et al. [2012]. Using only our ascending velocities and ignoring vertical motion, we obtain creep rates of ∼25 mm/yr in the Adapazari basin.
For the first 10 km east of Izmit, our fault creep estimates agree well with that of Cakir et al. [2012]. Farther east of this and into the Adapazari basin our estimates of fault-parallel creep decreases to ∼6 mm/yr, while HUSSAIN ET AL. IZMIT CREEP 2994 Figure 14. Projected east-west velocities along the profiles B-B' and C-C' indicated in Figure 6a in the paper. The dark and light grey areas represent the 95% and 68% confidence bounds for our preferred MAP solution shown in Figure 9. The green line shows the best fit result when we fix the fault slip rate and locking depth to the values used in Cakir et al.
[2012] (27 mm/yr for the slip rate and 12 km for the locking depth) and invert for the creep rate and creep depth. In blue is a forward calculation using the best fit parameter values for this location from Cakir et al. [2012]. Cakir et al.'s [2012] begin to increase to about 27 mm/yr. In their model this coincides with a deepening of the depth to the top (0 km to 1 km) and bottom (5 km to 12 km) of the creeping segment; i.e., fault creep is extending down to deeper depths but not reaching the ground surface. In our models the maximum depth of the creeping segment (equivalent to Cakir et al. [2012], bottom depth) remains fairly constant at ∼7 km. We do not solve for a top depth because offsets in the velocity gradients across the fault (section 3) imply that fault creep reaches the ground surface along the entire 80 km of the creeping segment; therefore, we implicitly assume a top depth of 0 km. Figure 14 compares the confidence range (in grey) for our MAP solution with the prediction from the best fit parameters estimated by Cakir et al. [2012] (in blue). It is clear that they overestimate the fault creep rate at both locations. The green line is the MAP result when the far-field slip rate and locking depth are fixed to the values used by Cakir et al. [2012]. The fit, although within our 95% confidence range, is poor at far-field locations compared to our solution, and the fit is also poor at these locations to existing GPS data. However, the creep rate and depth estimates are similar to our MAP results implying that these parameters are fairly insensitive to the far-field fault slip rate and locking depth.

Cakir et al.
[2012] processed their InSAR data using a persistent scatterer single master approach while we use the persistent scatterer small baseline technique. In the single master case, unwrapping errors can only be assessed visually while we use a more robust loop closure technique to systematically identify and remove (or fix) interferograms with unwrapping errors increasing the chances that the final LOS velocity is estimated using good quality data (see section 2). An interferogram network containing unwrapping errors would translate into an inaccurate LOS velocity map. In areas of poor coherence we recommend using the small baseline processing technique in order to minimize the influence of these errors in the final LOS velocity estimate.
HUSSAIN ET AL.

IZMIT CREEP 2995
Furthermore, our creep estimates from direct offset measurements at the fault ( Figure 5) agree well with the estimates from the profile modeling ( Figure 9). For profile A-A' , the direct offset method gives ∼8 mm/yr compared to 5 (2-8, 95% confidence interval (CI)) mm/yr from the model; for profile B-B' , the direct offset method gives ∼6 mm/yr compared to 5 (3-8, 95% CI) mm/yr from the model, and for profile C-C' , the direct offset give ∼3 mm/yr compared to 2 (0-13, 95% CI) mm/yr from the model. The corrected GPS displacement time series (Figure 12b) allows us to directly estimate the fault creep rate for the Envisat period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). From the GPS alone we estimate a creep rate of 9 ± 4 mm/yr, which compares well with our estimate of ∼6 mm/yr (from the direct offsets at location of profile B-B').
Our estimates for interseismic locking depth and fault slip rate show a clear trade-off ( Figure 11, top left). This is an expected feature of Savage and Burford [1973]-type elastic models [e.g., Wright et al., 2001a;Walters et al., 2011], where a larger slip rate can be compensated by a deeper locking depth. This trade-off is the main cause of the large confidence interval on these parameters. Therefore, we cannot discount a deeper locking depth and faster slip rate than our MAP estimates as these still fall within the confidence limit of our solutions. We have also tried a model for profile B-B' where we fix just the locking depth to 15 km. Our best fit results for the slip rate, creep rate, and creep depth are 20 mm/yr, 6 mm/yr, and 10 km, respectively. The slip rate and creep rate are similar to our MAP solution of 20 mm/yr and 5 mm/yr, respectively. The creep depth is a little deeper than the MAP estimate of 8 km.

Time Series
Our displacement time series model based on rate-and-state-dependent friction is consistent with both the GPS and InSAR displacements ( Figure 12). Note that this model assumes that the maximum coseismic slip is not collocated with the postseismic creep. This is required because in regions that have undergone coseismic slip the change in stress due to the earthquake, Δ , is negative and this rate-and-state model predicts relocking. However, if the maximum coseismic slip zone is not collocated with the creeping zone, then afterslip can occur. The assumptions hold in our case because the maximum coseismic slip occurred at depths of around 11 km [Feigl et al., 2002] while our observed creep occurs in the upper ∼8 km or less (Figure 9).
The model predicts that following the 1999 Izmit earthquake, the fault experienced a period of rapid afterslip with an initial maximum rate of 2.9 ± 0.4 mm/d or 1.1 ± 0.2 m/yr. This estimate is slightly lower than that of Bürgmann et al. [2002b] who found a maximum afterslip rate of ∼1.4 m/yr in this region using GPS alone. This discrepancy could be due to the simplicity of our model and the fact that our maximum afterslip velocity is a depth-averaged value and neglects any off-fault deformation. The marginalized probability distributions for the rate-and-state model parameters show a clear trade-off between the two model parameters (Figure 13), where a faster maximum afterslip, v max , rate can be compensated for by a smaller value of (a − b)∕k.
Over the 5 years after the earthquake, the afterslip rate decays to a rate of ∼5 mm/yr. Continued afterslip 5 years after the earthquake is not completely unexpected. Rubin and Ampuero [2005] define the end of the afterslip period to be the time at which the deceleration vanishes. Segall [2010] showed that an estimate of the duration of postseismic afterslip can be determined by where d c is the critical slip distance from rate-and-state-dependent friction laws, v plate is the long-term fault parallel loading rate, k c the critical stiffness, and k the stiffness. It is reasonable to associate k c ∕k with the ratio of the area of the rupture to that of the nucleation zone. Using laboratory estimates of d c [e.g., Scholz, 1988;Dieterich and Kilgore, 1994], the values for d c ∕v plate are of the order 10 4 to 10 5 s, so t post could range from months to years. Hearn et al. [2002] calculated a 1-3 MPa change in coseismic shear stress, Δ , for the uppermost 5 km of the crust at the fault location where we calculate the maximum fault creep. Assuming Δ is positive at steady state, we can use equation (7) [Segall, 2010] HUSSAIN ET AL. IZMIT CREEP 2996 From our displacement time series, we have (a − b)∕k = 51 mm, giving stiffness, k, equal to 3-8 MPa/m. The normal stress between 1 and 5 km is 30-130 MPa (assuming lithostatic pressures), which gives us an estimate for the rate-and-state parameter (a − b) of 0.001-0.014. Using a steady state fault creep rate of 5 mm/yr, we arrive at a frictional shear stressing rate of 0.015-0.04 MPa/yr. Our estimate for the rate-and-state parameter (a − b) compares well with the results of Kaneko et al. [2013] who determined a value between 0.002 and 0.008 for the Ismetpasa section of the NAF. Our estimates are larger than the range estimated for the Parkfield segment of the San Andreas Fault (0.0006-0.0018) by Johnson et al. [2006]. However, the authors of that study acknowledge that their estimates for (a − b) are low and explain it by assuming that the afterslip occurs in a transition zone between velocity weakening (negative a − b) and velocity strengthening (positive a − b). Laboratory experiments give values for a − b, in the range −0.0019-0.0070 for upper crustal rocks such as sandstones, slates, marbles, and granites [Ikari et al., 2011]. However, as noted by Marone [1998], laboratory friction experiments are highly idealized relative to natural faults, and it is not clear whether a comparison of laboratory-derived friction parameters are relevant to earthquake faulting.
The simple afterslip model used in this study has several important limitations. First, the model is 1-D and therefore assumes that the fault has homogeneous properties. There have been several successful attempts in describing the general aseismic behavior of faults and its spatial variability using 2-D spring slider fault models [e.g., Dieterich, 1992;Perfettini et al., 2003a] as well as new 3-D modeling techniques [e.g., Jolivet et al., 2015].
As noted by Perfettini and Avouac [2004], 1-D afterslip models like the one used in this study do not take into account the viscous relaxation of the deep crust or upper mantle. For large events like the M w 7.4 Izmit earthquake we can expect a transfer of stress to regions below the seismogenic zone.
Another limitation is that we do not consider the role of aftershocks after the Izmit earthquake, which may also accumulate slip.

Moment Release
The total moment released by the aseismic creep at the Izmit rupture between 1999 and 2012 (assuming fault creep extends down to 8 km along the 80 km length of the creeping segment) is of the order of 10 18 Nm. The total moment released in the Izmit earthquake was on the order of 10 20 Nm [Pinar et al., 2001;Delouis et al., 2002]. Bürgmann et al. [2002b] modeled the complete distribution of afterslip on and below the coseismic slip patch in the 87 day period after the Izmit earthquake and found that the total moment released due to afterslip on the entire fault was of the order 10 19 Nm. This is still only 10% of the coseismic moment release, but it would have altered the stresses in the surrounding crust and could have contributed to the triggering of the Düzce earthquake a few months later.
Currently, the near-steady state creep rate of 5 mm/yr is only releasing about 5-40% of the accummulating moment from plate loading (assuming creep extends down to 8 km depth and fault loading rate of 20 mm/yr (9-27, 95% CI) to a locking depth of 17 km (8-28, 95% CI)). Therefore, in terms of the moment budget, and thus stress transfer and seismic hazard, the effect of the shallow, aseismic slip in the past decade is small compared to that from plate loading.

Implications for the Shallow Slip Deficit
Inversions of coseismic surface deformation data from large (M w ∼7) strike-slip earthquakes indicate that the coseismic slip in the uppermost section of the seismogenic crust is often systematically less than that at seismogenic depths (4-10 km) [e.g., Simons et al., 2002;Fialko et al., 2005;Bilham, 2010;Dolan and Haravitch, 2014]. Despite variations between different inversion results, coseismic slip studies of the Izmit earthquake also show this so-called "shallow slip deficit." Delouis et al. [2002] and Feigl et al. [2002] slip inversions reveal that the average coseismic slip in the midcrust (between 6 and 15 km) beneath the western margin of Lake Sapanca was 5-7 m while in the shallow crust (<6 km) the average coseismic slip was 3-4 m, producing a coseismic slip deficit of 2-3 m in the shallow crust.
There are two possible explanations for the origin of the shallow slip deficit: the deficit arises due to distributed off-fault deformation in the uppermost crust, in which case it is not strictly speaking a deficit as such, or the shallow slip deficit is accounted for by aseismic slip on the shallow fault.
The mean recurrence time for earthquakes along the North Anatolian Fault is 200 years [Stein et al., 1997]. Projecting our afterslip model 200 years into the future, we find the total average displacement due to afterslip is 1-1.3 m, which is 30-65% of the slip deficit. Therefore, the total steady state afterslip (integrated over the HUSSAIN ET AL. IZMIT CREEP 2997 There is a clear discontinuity between the dominantly green and blue (−10 to −5 mm) pixels to the north and the red pixels to the south (7 to 13 mm) of the fault, consistent with fault creep with a right-lateral sense of motion across the fault. The distance between pixels either side of the fault that show a clear displacement discontinuity is around 100 m. earthquake cycle) is insufficient to account for all of the shallow coseismic slip deficit. If the deficit is to be explained by slip on the fault during the period between earthquakes (but not by afterslip), it must therefore occur during phases of transient creep. No such events have been observed at this location during the period of good geodetic observations (since ∼1995). Alternatively, the slip deficit could be due to distributed inelastic coseismic deformation in the uppermost few kilometers of the crust.
Whether this deformation occurs in the coseismic period or in the interseismic period is difficult to tell. Wright et al. [2001b] noted that the Izmit earthquake triggered shallow slip on a smaller fault south of the main rupture. Their coseismic interferogram also shows small off-fault movements on other local faults. The integrated slip from all the small triggered slip events since 1999 has not been clearly documented, but it is unlikely to be of the order of a few meters. Nevertheless, it could help explain a portion of the shallow slip deficit [Dolan and Haravitch, 2014].
Slow long-term distributed deformation [Fialko et al., 2005], possibly due to viscous creep mechanisms, is also difficult to resolve due to the noise in the InSAR data. Figure 15 shows that we can constrain the surface deformation to be occurring within ∼50 m either side of the fault, if not on the fault plane itself. Our fault creep time series uses GPS stations located ∼3 km from the fault. Thus, any off-fault deformation within this distance would have been subsumed and represented as fault creep in our time series model. Distributed deformation mechanisms active during the interseismic period in any case would be expected to produce deformation at longer wavelengths. A dense network of long-term continuous GPS measurements near the fault would help to constrain the role of viscoelastic relaxation.

Conclusion
We have investigated the spatial distribution of aseismic creep on the Izmit and Düzce sections of the North Anatolian Fault using Envisat ASAR images in both ascending and descending geometries between 2002 and 2010. Our results show a discontinuity in the LOS velocities at the surface trace of the 1999 Izmit earthquake, consistent with aseismic creep on the shallow fault. Fault creep on the Izmit rupture at the end of our InSAR observation period (2010) occurs at an average rate of ∼5 mm/yr and extends about 80 km from the Gulf of Izmit in the west to about 30.7 ∘ E in the east. However, it is likely that this extends farther west into the HUSSAIN ET AL.

IZMIT CREEP 2998
Gulf making our estimate of the length of the creeping segment a lower bound. We observe a time-averaged maximum fault creep rate of 11±2 mm/yr near the city of Izmit for the period 2002-2010. We also find that the Adapazari basin region is subsiding at a rate of about 6 mm/yr. The causes of this subsidence remain unclear but could be related to water pumping or tectonic subsidence related to movement on a nonplanar fault. The time series of relative displacement change between two points either side of the fault is consistent with an afterslip model that predicts a period of rapid deceleration from a maximum velocity of 2.9 ± 0.4 mm/d immediately after the earthquake to a near-steady state value of ∼5 mm/yr after about 5 years, implying that aseismic creep could continue for many years. The moment released by the shallow aseismic slip between the period 1999 and 2012 is small compared to the moment released by the Izmit earthquake. The current rate of moment release due to aseismic slip is about 5-40% of the rate of moment accumulation from plate loading. Therefore, we conclude that the NAF in this region is mostly locked and accumulating strain, and the long-term impact of shallow aseismic creep is negligible in terms of stress transfer and seismic hazard. Projecting our afterslip model 200 years into the future, we predict the total displacement due to shallow afterslip to be 1-1.3 m. This is insufficient to account for the 2-3 m slip deficit in the shallow crust (<6 km) observed in coseismic slip studies of the Izmit earthquake. Distributed inelastic deformation in the uppermost few kilometers of the crust or slip transients during the interseismic period are likely to be important mechanisms for generating the shallow slip deficit. project grant NE/I028017/1, which supports the lead author's research studentship via the FaultLab project at the University of Leeds. The Envisat satellite data are freely available and were obtained from the European Space Agency's Geohazard Supersites project. Our MCMC algorithm was developed using Aslak Grinsted's open source GWMCMC Hammer algorithm available on the mathworks website. The GPS data were obtained from the Global Strain Rate Model project website (http://gsrm.unavco.org). Most of the figures in this paper were made using the public domain Generic Mapping Tools (GMT) software [Wessel and Smith, 2001]. COMET is the Centre for the Observation and Monitoring of Earthquakes, Volcanoes and Tectonics. We would like to thank Romain Jolivet and Eric Lindsey for providing useful suggestions that helped improve the quality of the manuscript. Results can be obtained by contacting the lead author (eeehu@leeds.ac.uk).