Constraining crustal velocity fields with InSAR for Eastern Turkey: Limits to the block‐like behavior of Eastern Anatolia

The Sentinel‐1 satellite mission will enable global strain rate mapping from interferometric synthetic aperture radar (InSAR) and GPS data, and methods to combine these data in velocity fields will become increasingly important. Here we use InSAR to measure interseismic deformation in Eastern Turkey, across a ∼250,000 km2 area that spans the Arabia‐Eurasia plate boundary zone. From our InSAR data we first estimate slip rates and locking depths for the North and East Anatolian Faults (NAF and EAF) of 20 ± 3 mm/yr and ∼16 km and 11 ± 3 mm/yr and ∼16 km, respectively, but we also combine the InSAR data with existing GPS velocity measurements to construct high‐resolution velocity and strain rate fields across the region for the first time. We calculate 2‐D and 3‐D velocity fields and find that strain is mainly localized across the NAF and EAF and that there is negligible differential vertical motion across the Eastern Anatolian plateau. We also show that high‐resolution 2‐D strain rate fields can be calculated from InSAR alone, even in the absence of GPS data. We fit a block model to our velocity field and estimate slip rates of ∼21 mm/yr and ∼8 mm/yr for the NAF and EAF, showing that our previous estimates differ from these values because they neglected crustal rotation. Although this rotation is an important component of the velocity field in Eastern Turkey, systematic residuals between our velocity field and the best fitting block for Anatolia suggest that the region is not block‐like as proposed by previous authors.


Introduction
Interferometric synthetic aperture radar (InSAR) is a space geodetic tool that is well suited to mapping interseismic strain at high spatial resolution across large areas. While the earliest InSAR interseismic studies were based on small stacks of interferograms [e.g., Wright et al., 2001], current such studies commonly involve analysis of hundreds of interferograms spanning several overlapping tracks and across regions hundreds of kilometers wide [e.g., Jolivet et al., 2012;Kaneko et al., 2013;Tong et al., 2013;Cavalié and Jónsson, 2014]. New SAR missions such as Sentinel-1, with regular acquisitions and frequent revisit times, will enable InSAR to be used to map crustal deformation on the global scale, and to be incorporated into maps of global strain rate, which are currently based on GPS measurements alone [e.g., Kreemer et al., 2003]. The combination of GPS and InSAR measurements of crustal deformation will therefore become increasingly important in order to fully exploit these new data and produce global high-resolution velocity and strain rate fields.  developed a method to construct strain rate and velocity fields from InSAR and GPS data and tested it in Western Tibet to examine the role that the Altyn Tagh and Karakorum faults play in Tibetan crustal deformation. They found little strain localization on these major tectonic structures and found instead that strain was distributed more widely across Western Tibet. In this paper we employ similar velocity mapping methods in Eastern Turkey, where the North and East Anatolian Faults mirror the Karakorum and Altyn Tagh faults. Eastern Turkey is an ideal area for the application of these methods, as the fault slip rates are high and the spatial coverage of GPS measurements is relatively dense.

ANATOLIA E a s t A n a t o li a n F a u lt N o rt h A n a to li a n F a u lt
O v a c ik  Reilinger et al. [2006] (red), Ozener et al. [2010] (blue), Alchalbi et al. [2010] (purple), Yavaşoglu et al. [2011] (yellow), Reilinger and McClusky [2011] (black), Tatar et al. [2012] (green), and Aktug et al. [2012] (pink). All GPS velocities have all been rotated into the reference frame of Reilinger et al. [2006], using common stations, with the exception of data from Yavaşoglu et al. [2011]. All GPS data are shown in a fixed-Eurasia reference frame. Active strike-slip faults fromŞaroglu et al. [1992] are shown in black. The top inset figure shows the location of the study area within Turkey, and the bottom inset figure shows the interpolated horizontal GPS velocities that are used as an initial model to correct the InSAR data for orbital errors.
We use Envisat SAR data from 2003 to 2010 to measure interseismic strain accumulation across the whole width of the plate boundary zone in Eastern Turkey, from the Arabian plate in the south, across Anatolia, to the Eurasian plate in the north. We use data acquired on three descending and two ascending tracks to measure strain accumulation across both the NAF and the EAF and to determine whether or not there is significant strain accumulation away from the major faults. We first use a multi-interferogram network approach to create five overlapping maps of line of sight (LOS) velocity and then combine these rate maps with existing GPS measurements in order to derive a velocity field and a map of strain rate for Eastern Turkey. We also construct a 3-D velocity field across Eastern Turkey to investigate vertical deformation and show that InSAR data alone can be used to generate 2-D strain rate fields, even in the absence of GPS measurements. Finally, we use simple models to interpret our geodetic velocity field and investigate the nature of crustal deformation in the region.
WALTERS ET AL.

SAR Data and Processing
We processed 145 Envisat SAR acquisitions to form 373 interferograms across five tracks (Table 1), using the Jet Propulsion Laboratory (JPL)/California Institute of Technology (Caltech) ROI_PAC software v3.01 [Rosen et al., 2004]. We used Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) by satellite orbits from the European Space Agency to correct interferograms for orbital effects and removed topographic effects using a filled 3 arc sec digital elevation model (DEM) from the Shuttle Radar Topography Mission [Farr et al., 2007], obtained from Consultative Group on International Agricultural Research-Consortium for Spatial Information (http://srtm.csi.cgiar.org). In order to extract the maximum coherence from each interferogram while not introducing unwrapping errors, we employed different levels of spatial averaging for different interferograms. The number of samples taken varied from 4 to 64, but we processed the majority of interferograms at 32 or 64 samples (∼640-1280 m pixels). We used a weighted power spectrum filter [Goldstein and Werner, 1998] to improve the signal-to-noise ratio and unwrapped phase using the branch-cut method [Goldstein et al., 1988]. We used phase closure around closed loops of interferograms to identify any unwrapping errors introduced during processing, and these were corrected manually [e.g., Biggs et al., 2007].

Atmospheric Correction of Interferograms
Phase delays from spatiotemporal changes in atmospheric refractivity are the largest source of error when measuring small ground deformations with InSAR. To estimate the magnitude of the atmospheric signals that we should expect to see in interferograms produced from our data set, we estimated phase delay maps using the ERA-Interim (ERA-I) global atmospheric model reanalysis product [Dee et al., 2011] obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF). To estimate the combined delay from differences in water vapor (wet delay) and atmospheric pressure (hydrostatic delay), we first created delay maps using the method of Doin et al. [2009] andJolivet et al. [2011]. The ERA-I atmospheric model provides estimates of 16 meteorological parameters including temperature, relative humidity, and geopotential at 6 h intervals from 1989 onward. Spatially, the outputs are given at grid nodes with ∼75 km horizontal spacing, and at 37 pressure levels between ∼0 and 50 km above sea level. We calculate atmospheric delay maps by producing vertical profiles of total path delay at ∼75 km resolution using the atmospheric delay formulations of Baby et al. [1988], interpolating these horizontally and vertically with bilinear and spline interpolation functions respectively and then truncating these vertical profiles of delay with a DEM to find the delay at the correct elevation for each point at the surface. Delay maps are linearly interpolated in time from the closest delay map before and the closest delay map after the SAR acquisition time.
We produced differential delay maps for each possible interferogram (i.e., each combination of acquisition dates), interpolated at 1.6 km resolution, and then calculated the maximum spatial variation in differential delay ( ij ) across each map [e.g., Walters et al., 2013]: where kl ij is the double difference in atmospheric delay between SAR acquisition dates i and j, and spatially between pixels k and l. For each interferogram, ij is an estimate of the expected maximum atmospheric signal, as the spatial variability of differential delay across the scene is the important factor, not absolute values. The mean wet ij across all possible interferograms on all tracks was ∼10 cm, and ∼3 cm for the hydrostatic ij ( Table 2). The mean total (wet + hydrostatic) ij was ∼10 cm for any one interferogram, significantly larger than the tectonic signal expected in any of our interferograms. We subtracted the total ERA-I atmospheric correction (wet + hydrostatic) from each of the 373 interferograms and then removed a best fitting linear orbital plane from each corrected interferogram. We compared the RMS of this corrected, flattened interferogram with the RMS of the interferogram with the same orbital plane removed, but with no ERA-I correction applied. In this case the RMS reductions are on average 73%, 71%, 65%, 66%, and 73% for tracks 78, 307, 35, 171, and 400, respectively. Overall, for all interferograms we calculate a 70% reduction in RMS. This is similar to the 73% quoted by Jolivet et al. [2011]. However, this estimate of RMS reduction could be considered as biased because the orbital trend that is removed from both corrected and uncorrected interferograms is derived from flattening the corrected interferograms. In the majority of InSAR studies a planar (or higher-order) empirical orbital correction to interferograms is applied irrespective of whether the atmosphere has been corrected for or not.
An alternative method is to compare the RMS of the ERA-I corrected and flattened interferograms with the uncorrected interferograms flattened independently with a best fitting linear orbital plane. In this case, the RMS reductions are 5%, 3%, 9%, and 1% for tracks 78, 307, 35, and 400, respectively, while track 171 saw an RMS increase of 9%. Overall for all interferograms, the RMS reduction is 4%. This result suggests that the ERA-I correction only makes a small improvement in RMS reduction over making an empirical linear orbital correction alone. This is more similar to the 9% reduction in RMS observed by Garthwaite et al. [2013] from applying ERA-I corrections to 51 interferograms over central Tibet. These results do not necessarily imply that the ERA-I delay maps are incorrect and instead could suggest that most of the ERA-I atmospheric correction is long wavelength which would otherwise be absorbed by the empirical orbital correction. We suggest that although the improvement from using ERA-I corrections seems small, for data sets comprising of several tens to hundreds of interferograms it may be worth adopting routinely as on average it removes some topographically correlated signals that cannot be removed by linear orbital correction and will also likely result in more accurate estimation of orbital errors. Therefore, in the following section we correct all interferograms with ERA-I-derived delay maps but also discuss briefly the effect on our results of not applying this correction.

Rate Map Estimation for Eastern Turkey
We use a network-based approach to mitigate residual orbital errors in all interferograms and to combine all interferograms for each track to estimate interseismic LOS velocities and associated uncertainties. This method is based on the methods developed by Biggs et al. [2007], Elliott et al. [2008], and Wang et al. [2009] and to implement it we use the π-RATE (Poly-Interferogram Rate And Time-series Estimator) MATLAB-based software package . We correct orbital errors on an epoch-by-epoch basis by removing a six-parameter quadratic phase ramp from each interferogram, and in order to avoid the erroneous removal of long-wavelength deformation signals, we remove a GPS-derived a priori initial deformation model before orbital error estimation and then add back the subtracted model before rate map estimation [e.g., Wei et al., 2010;Garthwaite et al., 2013]. We assume that velocity is constant over the total time span of all our InSAR observations. More details of the method are presented in Appendix A.
The resulting rate maps are shown in Figure 2. In all five rate maps there is a clear change in LOS velocity across the NAF and the EAF that is qualitatively consistent with right-lateral and left-lateral strike-slip interseismic motion, respectively. We also calculate how well neighboring rate maps agree in overlapping areas and plot histograms of the velocity differences in Figure 2 (bottom row). We first make a correction for variable LOS based on the assumption of horizontal motion only, by dividing velocities by the sine of the local incidence angle, and then multiplying by the sine of 23 • , the incidence angle at the center of each track.  The three histograms have mean values of approximately zero and have standard deviations of the best fitting Gaussian distribution of 2.12, 1.70, and 2.61 mm/yr for T78-T307, T307-T35, and T171-T400, respectively. These values can be considered as estimates of √ 2× the velocity uncertainties for the individual rate maps, giving 1σ values of 1.50, 1.20, 1.85, and 1.85 mm/yr for T78, T35, T171, and T400, respectively. T307 overlaps with both T78 and T35, so we take the mean of the 1σ values for T78 and T35, which is 1.35 mm/yr.
We also investigated the effect of our choice of orbital and atmospheric correction on the results. We calculated sets of rate maps under three different additional combinations of processing choices; linear orbital correction and ERA-I atmospheric correction, linear orbital correction and no atmospheric correction, and quadratic orbital correction and no atmospheric correction. In general, these choices do not significantly affect our results, and the rate maps for tracks 78, 307, 35, and 171 are very similar for all combinations of processing choices. Track 400, however, which has the least number of interferograms, does not show a clear tectonic signal across the NAF or the EAF under any of these three additional processing choices and is strongly dissimilar to track 171 where the two tracks overlap. In general, using a linear orbital correction and/or no atmospheric correction results in a greater discrepancy between overlapping rate maps and larger standard deviations of the histograms of velocity differences. For the quadratic orbital correction, the ERA-I correction decreases the Gaussian standard deviations for T78d-T307d by 3%, T307d-T35d by 23%, and T171a-T400a by 24%. This is lower than the 73% improvement suggested by Jolivet et al. [2011], but higher than the estimate of only 4% made in the previous section. This may be because in the previous section we did not take into account the covariance of orbital corrections and instead estimated them . Rate map profiles across the NAF and EAF, assuming all motion is horizontal and fault parallel. (left, first and second panels) Rate map pixels from all descending data sets projected onto profile lines perpendicular to the NAF and EAF. Data points are colored by track, and each track has a correspondingly colored label in the rate map mosaic (right). Model profile lines are for a slip rate of 20 mm/yr and a locking depth of 16 km for the NAF, and for a slip rate of 11 mm/yr and a locking depth of 16 km for the EAF. (left, third and fourth panels) The same as for the top profiles but instead showing the two ascending data sets. Descending data points are shown in grey in the background. GPS velocities and 2σ error bars are overlain on the profiles. independently for each interferogram, whereas here we are using a network approach which will estimate orbits self-consistently.

Rate Map Inversion and Modeling
We first analyze profiles across each fault and fit these with a simple 1-D model. We project data from all five rate maps onto fault-perpendicular profiles for both the NAF and EAF. The profiles for all five data sets are consistent with each other ( Figure 3) and with the GPS velocities. We calculated a mean profile and 1σ bounds for each data set by inverting all the rate map data within 10 km along-profile bins, weighting the inversion using a spatial variance-covariance matrix (VCM) to account for spatial correlation between rate map pixels. The spatial VCM was estimated by fitting a 1-D autocovariance function of the form to signals in a 100 km × 100 km nondeforming region of each rate map, where C ij and d ij are the covariance and distance between pixels i and j, A 2 is maximum variance, and b is 1-D e-folding distance [Wright et al., 2003;Lohman and Simons, 2005;Parsons et al., 2006]. However, while the estimates of b are reasonable at around 10 km, the estimates of A from this method are 0.5-1.4 mm/yr and are all lower than the 1.20-1.85 mm/yr 1σ values estimated from the differences between rate maps. We therefore still use the b values of ∼10 km but use the variance estimates derived from overlapping regions of the rate maps.
We use the assumption that all deformation is fault parallel and horizontal and model the fault as a buried infinite screw dislocation in a homogeneous, isotropic elastic half-space, where during the interseismic interval, aseismic strike-parallel slip occurs at a rate s below a locking depth d. The rate of horizontal displacement y at a perpendicular distance x from the fault is y = (s∕π) × arctan(x∕d) [Savage and Burford, 1973]. For each profile separately, we varied the locking depth between 0 and 40 km at 1 km intervals, and for each locking depth, found the best fit slip rate and adjustable additive constant from a linear inversion that minimized the total RMS misfit between the model and data profiles.
The results are presented in columns 2 and 3 of Tables 3 and 4. The variability in locking depth is large, so we also repeat the inversions but use the same locking depth for all tracks; i.e., we assume that locking depth does not vary along the sections of the NAF and EAF in this region. To find the best locking depth for each fault, we project all the descending data onto one profile line and then calculate the mean profile line and find the best fit slip rate and locking depth using the method above. From this technique, the best fit locking depth for both the NAF and EAF is 18 km. The resulting slip rates from fixing the locking depths are presented in column 4 of Tables 3 and 4.
The mean and standard deviations of the five estimates of slip rate and locking depth for each fault can be taken as the average estimates and the corresponding uncertainties for the eastern NAF and the EAF. Fixing the locking depth for all tracks does not change the mean slip rate from the mean obtained with variable locking depths but does reduce the standard deviation by a small amount for both faults.
The T400a rate map is produced from the least amount of data and is likely to be the least reliable rate map, as shown in the previous section by its sensitivity to different choices of orbital and atmospheric correction. The mean slip rate and locking depth on the NAF are 20 ± 3 mm/yr and 18 ± 9 km for all five rate maps. The mean slip rate and locking depth excluding the T400a rate map are 20 ± 3 mm/yr and 16 ± 10 km, which perhaps are more representative of the true slip rate and locking depth.
The mean slip rate and locking depth on the EAF are 11 ± 3 mm/yr and 15 ± 5 km from all tracks, and 11 ± 3 mm/yr and 16 ± 5 km without track 400. Tables 3 and 4, there is no clear spatial pattern of slip rate variation along either the NAF or the EAF. Variation in slip rate between neighboring tracks is more likely  Table 6) are shown by red vectors, evaluated at nodes of the model mesh in a fixed-Eurasia reference frame. Error ellipses for both sets of velocity vectors are presented at the 2σ level. (middle) Model mesh, overlain with data, model, and residual (blue) velocity vectors. (right) From top to bottom, InSAR rate maps, model rate maps, and residuals for preferred velocity field model. The black cross in the top row is placed in the same reference location in each panel.

From the slip rates in
reflective of the level of precision to which these estimates can be made. Assuming that the slip rate does not vary along the sections of the NAF and EAF covered by these five InSAR tracks, the uncertainty estimates on locking depth and slip rate presented here can be considered realistic levels of precision for InSAR studies of this type involving tens of SAR acquisitions per track. With data sets of this size a slip rate estimate from a single track alone is likely to differ from the true slip rate by up to 3 mm/yr.
Our slip rate of 20 ± 3 mm/yr for the NAF is slightly lower than the 25-26 mm/yr estimated by Reilinger et al. [2006] for this section of the NAF, but in agreement with the more recent GPS results of Ozener et al. [2010] and Tatar et al. [2012], who find 16-24 mm/yr and 20 ± 2 mm/yr from denser GPS networks across the fault. Our slip rate is also compatible with the previous InSAR results of Wright et al. [2001] and Walters et al. [2011], who find ranges of 17-32 mm/yr and 20-26 mm/yr. For the EAF, our slip rate of 11 ± 3 mm/yr is in agreement with the 10 mm/yr estimated by Reilinger et al. [2006] and the 14 ± 2 mm/yr estimated by Westaway [1994]. Our slip rate is also in agreement with a recent InSAR estimate of 13 mm/yr further to the east along the EAF [Cavalié and Jónsson, 2014]. However, this study shows a much shallower 4.5 km locking depth for the EAF than we find here, possibly indicating significant variation in locking depth toward the junction with the NAF, or local earthquake cycle-related transient effects.

A Velocity Field for Eastern Turkey
We combine the GPS velocities shown in Figure 1 with our five InSAR rate maps ( Figure 2) in order to create a combined geodetic velocity field for Eastern Turkey. We use the methods of , who adapted the technique of Molnar [1997, 2005] to incorporate InSAR LOS velocities as input for the velocity field inversion. We divide the region of interest into a mesh of arbitrary spherical triangles (Figure 4, middle), making the assumption that strain rate is constant within each triangle.
We perform the velocity field inversion as described in Appendix B using the VELMAP package of MATLAB codes . Our GPS data are the combined set of vectors discussed in section A1 and shown in Figure 1, and our InSAR data are the five rate maps shown in Figure 2, resampled to 3.7 km pixels. Our methodology can be used to calculate both 2-D and 3-D velocity fields. Initially, we adopt the common assumption that there is no vertical motion across the region of interest, but we also investigate the validity of this assumption in section 5.2.
The inversion is regularized using Laplacian smoothing; the strength of which is controlled by a smoothing factor 2 (see Appendix B), and we jointly invert the GPS and InSAR data using a range of 2 values between 1 and 0.01. There is a trade-off between RMS misfit of the solution to the data and the roughness of the  [Ekström et al., 2012]. The focal mechanisms within the northern dashed rectangle have been colored blue, and their moment tensors have been summed to give the large blue focal mechanism labeled NAF. The focal mechanisms within the southern dashed rectangle have been colored red, and their moment tensors have been summed to give the large red focal mechanism labeled EAF. The two pairs of horizontal strain rate axes result from the Kostrov summation of the colored focal mechanisms within each rectangle. The black arrows represent compression, and the grey arrows represent extension. The large focal mechanisms and the strain axes for each region have been normalized to be the same size, so only represent style of seismic strain release, not rate. solution. For our preferred solution (model V1, see Table 6 for details of velocity fields and block models presented in this paper) we choose 2 = 10 −1.0 , which corresponds to a compromise between misfit and solution roughness ( Figures S1 and S2 in the supporting information). The preferred solution has GPS residuals that are small (1.6 mm/yr RMS). The largest residuals are concentrated around the structurally complex junction of the NAF and EAF, where several GPS data sets have been merged together (Figure 4, middle). These residuals could therefore either be due to mismatches and disagreement between different GPS data sets in this region or from short-wavelength variations in velocity that cannot be fit by a smoothed velocity field. The InSAR data are also fit well by the preferred solution, with a 1.2 mm/yr RMS misfit between data and model rate maps (Figure 4, right).
The velocity field for 2 = 10 −1.0 is presented in Figure 4 (left), and the strain rate field for this model is shown in Figure 5 (top left). Strain accumulation is strongly localized on the NAF and the EAF. The Anatolian Plateau and the Turkish-Syrian border region to the south of the EAF are straining at relatively low rates in comparison to the two major strike-slip faults.

A Strain Rate Field From InSAR Alone
While Eastern Turkey has a very high density of GPS velocities, most continental tectonic regions worldwide do not have similar coverage. We use our data to test if geodetic strain maps can still be produced in such regions, with only sparse or even no input from GPS data. We rerun our velocity field inversion three additional times: once with input from GPS data alone (model V2), once with input from InSAR plus three GPS sites only (model V3), and once with input from InSAR and from three "dummy GPS" velocities obtained from the Mid-Ocean Ridge Velocity (MORVEL) plate model [DeMets et al., 2010] (model V4). For the dummy velocities, we fix the northern two velocities on the Eurasian plate to be zero and derive the southernmost velocity from the Arabia-Eurasia Euler pole from MORVEL. For the region where ascending and descending InSAR data overlap, the strain rate fields in Figure 5 show that the InSAR plus three GPS (model V3) and the InSAR plus MORVEL (model V4) solutions reproduce the same first-order features as the GPS-only solution (model V2) that strain is strongly focused on the NAF and EAF. The joint InSAR and GPS strain rate field (model V1) combines features from the "InSAR-only" (models V3 and V4) and GPS-only (model V2) solutions. We also compare the velocity fields from the InSAR plus three GPS (model V3) and GPS-only (model V2) solutions (Figure 9, top left) and find that the RMS of the differences in east velocity is only 1.4 mm/yr for the region where ascending and descending InSAR data overlap. The mean 2σ uncertainty for the difference in east velocity between these two velocity fields (models V2 and V3) is 1.9 mm/yr for this region. We acknowledge that the InSAR plus MORVEL solution is not entirely independent of the GPS data, as the GPS data were used to constrain orbital errors on the InSAR data. However, orbital uncertainties will be smaller for the new Sentinel-1 SAR mission than for Envisat due to on board GPS positioning [Torres et al., 2012] Figure 7 (right), the background color shows differential vertical velocities from the inversion for 3-D velocities, where north-south velocities are constrained by GPS alone. In both panels, the solid black outline shows the extents of all the InSAR rate maps, and the highlighted and dashed trapezoidal region shows the area covered by both ascending and descending InSAR data. The average uncertainty in vertical velocity across the trapezoidal region is 2.6 mm/yr. constraint may not be required in future. In this case we have shown that it will be possible to create such 2-D strain rate maps from InSAR alone across all continental regions, even in the absence of GPS data.

Vertical Motions
The results presented above have been based on the assumption that there is no differential vertical movement across the NAF and EAF. For the NAF, this assumption of horizontal movement is supported by Walters et al. [2011], who found that there was minimal vertical or fault-perpendicular motion across the fault. We use the Kostrov relation [Kostrov, 1974], to show that there is no appreciable non-strike-slip component to the seismicity over the last 30 years. Figure 6 shows the results of a Kostrov moment summation across Eastern Turkey [Kostrov, 1974]. In the Kostrov summation, the ijth component of the strain rate tensoṙi j from N cumulative earthquakes within a volume V over a time interval T is given by the following: where is the shear modulus, taken as 3.23 × 10 10 Pa, and M n ij are the components of the earthquake moment tensor. The seismogenic thickness is taken to be 15 km, and the time interval covered by the seismic record is only 36 years. However, the magnitude of the strain rate axes shown in Figure 6 have both been normalized to the same total strain rate, as the observation time interval is too short to enable robust spatial comparisons of strain rate magnitude to be made. Therefore, the choice of V, T, and do not affect the style of strain release shown in Figure 6.
We also use both ascending and descending InSAR data, along with existing GPS data, to measure differential vertical motion across the Arabia-Eurasia plate boundary zone in Eastern Turkey for the first time. We carry out a separate inversion where north-south velocities in the solution are constrained by the GPS alone, east-west velocities are constrained both by the InSAR and GPS data, and vertical velocities are constrained by the InSAR data alone (model V5, Figure 7). Since InSAR can only measure relative ground motions, we set the mean vertical velocity to zero. For the region covered by both ascending and descending InSAR data, which is the region in which this inversion is well constrained, the standard deviation of vertical velocities is 0.7 mm/yr, indicating that there is no significant relative vertical motion. Also, east-west velocities within this region are almost identical to those from our preferred solution, so we therefore conclude that our earlier assumption of no relative vertical motion across this region is valid.

Discussion
The strain fields in Figure 5 show that deformation across Eastern Turkey is concentrated on the NAF and EAF. Profiles through the preferred velocity field confirm this with little deformation occurring away from these narrow linear structures. In the profile in Figure 8, the linear gradients away from the faults (e.g., between the NAF and EAF in profile A-A ′ ) do not necessarily indicate strain concentration and instead are the result of crustal rotation. However, the sharp changes in gradient in this profile across the NAF and EAF cannot be explained by the rotation of a single region and therefore indicate strain accumulation. The Ovacik fault does not appear to be resolved as accumulating strain in Figure 5, at least above the uncertainty of the velocity field (∼1 mm/yr). [2012] combine InSAR and GPS data to create a high-resolution strain map of Western Tibet and find that the majority of focused strain does not coincide with major structures such as the Altyn Tagh and Karakoram faults. They therefore suggest that the block modeling approach that is frequently used to predict seismic hazard [e.g., Meade, 2007;Thatcher, 2007] is likely to overestimate slip rates on the Altyn Tagh and Karakoram faults and underestimate slip rates on other structures, if used in Western Tibet.

Wang and Wright
Our results here show that the converse is true for Eastern Turkey; unlike in Western Tibet, strain is mostly focused on the major faults. Therefore, we next consider the hypothesis that block-like behavior may be an appropriate description of the kinematics of this region and attempt to fit our velocity field (model V1) using a simple block model (model B1). We use the DEFNODE software [McCaffrey, 1995[McCaffrey, , 2002 to invert the velocities at the mesh nodes for interseismic coupling and tectonic block rotation. For comparison with the results of Reilinger et al. [2006], we use the same five blocks defined in that study for this region (Figure 8), but only solve for Euler poles for the three major blocks which cover our velocity field (Arabia, the Black Sea, and Anatolia). We fix the Euler poles for the other two blocks (Sinai and Iran) to those reported by Reilinger et al. [2006]. We fix the locking depth for the NAF and EAF to both be 15 km and the fault geometry and locking depths for all other block boundaries to the parameters in Reilinger et al. [2006].
The locations of our Euler poles for Anatolia and The Black Sea are in good agreement with the locations calculated by Reilinger et al. [2006] (to within 1σ and 2σ, respectively, Table 5), although our rotation rate for Anatolia is slower than the estimate presented in that study. Our Euler pole for northern Arabia is similar to the estimate from Reilinger et al. [2006] but does not agree within uncertainty. However, it might not be expected that the exact locations and rotation rates of the Euler poles should agree, since our study uses significantly more geodetic data (both additional GPS and InSAR) and also defines Euler poles for these regions based only on velocities from smaller geographic regions. The estimated slip rates across the NAF and EAF are shown in Figure 8 (top right). At a longitude of 38 • E, the NAF has a slip rate of 21.1 ± 0.1 mm/yr and the EAF has a slip rate of 7.5 ± 0.1 mm/yr. Both faults have estimated fault-perpendicular rates of less than 1 mm/yr. The formal uncertainties on these slip rates are underestimated because they only reflect the uncertainties in the relative Euler vectors and not those from assumed model parameters such as fault locking depth, location, and dip. The true uncertainty in slip rate is likely to be closer to the slip rate uncertainties estimated in section 4 (∼3 mm/yr). Although our two estimates for the slip rate on the NAF still agree to within 1σ uncertainties and the two estimates for the EAF agree to within 2σ uncertainties, the block model estimates for the NAF and EAF are slightly higher and lower, respectively, than the estimates from section 4. This is likely because the modeling in that section neglected the effects of rotation of Anatolia and Arabia. This can be seen in the profile in Figure 8, where a forward model of strain accumulation across the NAF and EAF which neglects crustal rotation shows significant misfit to the data. We suggest that the neglection of this rotation is also likely to bias the slip rate estimates for the NAF and EAF from InSAR data presented by Cavalié and Jónsson [2014]. In particular, the estimate of 13 mm/yr strike-slip motion across the EAF from that study is likely to be an overestimate of the true slip rate. The RMS misfit of the block model to the velocity field is 3.0 mm/yr for Anatolia, which is larger than the formal uncertainty for the velocity field in this region. In particular, there are large residuals in Anatolia (∼10 mm/yr) around the eastern end of the NAF  Savage and Burford [1973], with slip rates and locking depths of 20 mm/yr below 16 km and 11 mm/yr below 16 km, respectively. This model clearly neglects the velocity gradients across Anatolia and Arabia which correspond to rotation of these regions.
that are orientated perpendicular to the fault (Figure 8). These residuals fall outside the area covered by our InSAR data, but they show that although the block model predicts significant convergence across the fault at its eastern end, this is not supported by our velocity field and by the GPS data. We therefore suggest that although rotation is necessary to describe the overall kinematics of Eastern Anatolia, it cannot be well described as a single rigid block as proposed by previous authors [Reilinger et al., 2006;Cavalié and Jónsson, 2014], due to large systematic misfit to the best fitting block model.
We also investigate whether the velocity field obtained from InSAR and only three GPS sites in section 5.1 (model V3) can be used to obtain the same block model results. We rerun our block model inversion two additional times with different velocity fields as input: once using the velocity field derived from InSAR plus only three GPS sites (model V3) and once using the velocity field derived from GPS data alone (model V2). We only use the velocity field data in the region of InSAR coverage for both inversions (Figure 9). For both resulting block models (models B3 and B2, see Table 6) the estimates of slip rates across the NAF and EAF (Figure 9, bottom) and Euler poles for Anatolia and Arabia (Table 5) are similar to the results from our joint InSAR and GPS block model (model B1, Figure 8). At a longitude of 38 • E, the strike-slip rates estimated across the NAF and EAF are ∼20 mm/yr and ∼6 mm/yr for the block model derived from InSAR plus three GPS only (model B3), and ∼22 mm/yr and ∼7.5 mm/yr for the block model derived from GPS alone (model B2).
In section 5.2 we produced a high-resolution vertical velocity field across Eastern Turkey and found that there is no significant relative vertical motion across the region. This may appear surprising considering the topographic variation across the Anatolian Plateau (Figure 1, see top inset) but is also in agreement with recent seismicity, which records no significant vertical motion and also no thrust faulting in this region. The lack of thrust faulting is also be consistent with recent studies that suggest the Turkish-Iranian Plateau experiences dynamic support [e.g., Zor et al., 2003;Şengör et al., 2003], and the suggestion that the Plateau's elevation is due to dynamic uplift from a region of hot, low-density underlying mantle rather than crustal thickening [Şengör et al., 2003;Copley and Jackson, 2006]. The results from this study suggest that either this episode of uplift has finished or that the uplift is sufficiently slow or long wavelength that it is undetectable in ∼500 km long InSAR swaths. The lack of relative vertical motion across Eastern Turkey and the lack of convergence across the EAF confirms that there is no "Arabian push" that contributes to Anatolian motion [e.g., Reilinger et al., 2006] and that Anatolian crustal motion is driven only by forces relating to subduction at the Hellenic trench, either trench "suction" [e.g., Conrad and Lithgow-Bertelloni, 2004] or differences in gravitational potential energy between the Anatolian Plateau and the trench [e.g., McKenzie, 1972].
The localization of strain in Eastern Turkey has previously been used to support the concept that continental crust deforms by the differential motion of rigid, nondeforming blocks, with Eastern Anatolia defined as one such block [e.g., Reilinger et al., 2006;Cavalié and Jónsson, 2014]. However, our results show that there are large residuals to the best fitting block for Eastern Anatolia, implying that this region deforms. Molnar and Dayem [2010] suggest that major strike-slip faults, defined as those with slip rates >10 mm/yr, form adjacent to rigid regions because they are shallow, brittle manifestations of continuous strain in the lower crust and upper mantle, which necessarily

10.1002/2013JB010909
concentrates near discontinuities in strength. They also suggest that these major strike-slip faults are likely to have narrow shear zones at depth because of this concentration of strain. If we regard the Arabian shield and stable Eurasia as rigid, then Eastern Anatolia is bordered by two regions of strong crust, and strain concentration along these boundaries could cause Anatolia to have a low strain rate, without requiring that it is rigid. In addition, the fact that strain is concentrated around the NAF and EAF decades after the most recent large earthquakes supports the requirement of multiple relaxation timescales for the lower lithosphere, suggesting either a Burgers rheology [e.g., Hetland and Hager, 2005] or that these major faults are underlain by weak shear zones [Takeuchi and Fialko, 2013;Yamasaki et al., 2014].

Summary
We have generated five overlapping maps of LOS velocity in Eastern Turkey using InSAR, covering both the NAF and EAF. We have used these data to calculate slip rates for the NAF and EAF of 20 ± 3 mm/yr and 11 ± 3 mm/yr respectively, with associated locking depths of 16 ± 10 km and 16 ± 5 km. We have also used these InSAR rate maps together with a compilation of GPS data in order to calculate velocity and strain rate fields for Eastern Turkey. These show that strain is mainly localized across the NAF and EAF, with little strain accumulating in the Eastern Anatolian Plateau or in the northern Arabian plate. We have also calculated slip rates for the two major faults by fitting a block model to our velocity field, which gives slip rates of 21 ± 0.1 mm/yr and 8 ± 0.1 mm/yr on the NAF and EAF, respectively. The EAF slip rate is lower than that estimated from our velocity profiles as the block model accounts for crustal rotation. Although crustal rotation is an important component of the velocity field in Eastern Turkey, systematic residuals between our velocity field and the best fitting block for Anatolia suggest that the region is not block-like as proposed by previous authors. We have shown that there is negligible differential vertical motion across the NAF and EAF, which implies no active shortening in Anatolia, and the lack of an Arabian push driving the westward motion of Anatolia. In addition, we have demonstrated that while GPS data are spatially dense in this region and therefore can be combined with InSAR to recover 3-D deformation, high-resolution 2-D strain rate fields can be calculated from InSAR data alone, even in the absence of GPS data.

Appendix A: Rate Map Estimation With -RATE
In this study we use a network-based approach to mitigate residual orbital errors in all interferograms and to estimate interseismic LOS velocities and associated uncertainties from these InSAR data for each track. This method is based on the methods developed by Biggs et al. [2007], Elliott et al. [2008], and Wang et al. [2009], and to implement it we use the π-RATE (Poly-Interferogram Rate And Time-series Estimator) MATLAB-based software package .

A1. Estimation of Orbital Errors
We estimate residual orbital errors on an epoch-by-epoch basis using a network approach and a quadratic orbital error correction, whereby a six-parameter quadratic phase ramp is used to correct each interferogram [e.g., Wei et al., 2010;Garthwaite et al., 2013]. In order to avoid the erroneous removal of long-wavelength deformation signals, we remove a GPS-derived a priori initial deformation model before orbital error estimation and then add back the subtracted model before rate map estimation [e.g., Wei et al., 2010;Garthwaite et al., 2013]. This model contains predicted LOS velocities at each pixel, converted from horizontal GPS velocities in a fixed-Eurasia reference frame. The GPS data set is a combination of seven different GPS studies published between 2006 and 2012 ( Figure 1). In order to ensure consistency between the different data sets, we rotated all of the data sets into the fixed-Eurasia reference frame of Reilinger et al. [2006], using common GPS sites. For each of these common sites, we estimated the east and north velocity differences between the rotated velocities and those from Reilinger et al. [2006]. For all common sites, the mean of the absolute differences in north and east velocity are 0.53 mm/yr and 0.68 mm/yr respectively, which are of comparable magnitude to the stated uncertainties for most of the GPS measurements. One data set, that of Yavaşoglu et al. [2011], does not have any sites in common with any of the other data sets but was presented in a fixed-Eurasia reference frame, and so we have not applied a rotation to the velocities from that data set. In the combined data set, multiple measurements made at common stations were combined by calculating the mean velocity, weighted by the stated uncertainties of each measurement. The horizontal velocities were interpolated using a Delaunay triangulation of GPS sites followed by bilinear interpolation within each triangle from the velocities at the vertices.

10.1002/2013JB010909
These gridded velocities were then filtered with a 50 km diameter Gaussian filter and converted to LOS velocities using the local satellite incidence vector.

A2. Inversion for LOS Velocity
After correction for orbital errors, inversion of interferograms for LOS velocity is performed on a pixel-by-pixel basis, using graph theory to select the optimum combination of interferograms to use [Feigl and Thurber, 2009;. The stacking method uses Kruskal's algorithm [Kruskal, 1956] to find a minimum spanning tree for each pixel. Each branch corresponds to an interferogram and is assigned a length equal to the fraction of coherent pixels in the interferogram. The minimum spanning tree is therefore the network that connects all epochs using the most coherent interferograms, without including any closed loops of interferograms. Kruskal's algorithm is run on a pixel-by-pixel basis, so each pixel in the final rate map depends on different combinations of interferograms, and not all interferograms are required to be coherent for a given pixel. The inversion is weighted using a temporal variance-covariance matrix (VCM), which accounts for both the variance of individual interferograms and the covariance between interferograms that share a common acquisition date [Emardson et al., 2003]. The elements of the temporal VCM are the following: where 2 i is the variance of epoch i and 2 ij−i ′ j ′ is the covariance between the interferogram formed from epochs i and j and the interferogram formed from epochs i ′ and j ′ . The variances for each epoch are estimated from the variances of all interferograms in the network.
During rate map inversion π-RATE uses three thresholds, here termed Ψ, Ξ, and Υ, to exclude outliers from the inversion itself and from the final rate map: 1. If any pixel has less than Ψ coherent interferograms, then the pixel is excluded from the rate map. 2. The inversion for each pixel's velocity is performed on an iterative basis. First, a constant velocity for the given pixel is estimated from all interferograms, assuming a time-linear fit to all interferometric displacements. Then, if the largest residual to this linear fit for any interferogram is greater than Ξ times the standard deviation of that observation (from the temporal VCM), then that interferogram is excluded and the velocity is estimated again. This continues until all interferograms have residuals of less than Ξ times their a priori standard deviations, or there are less than Ψ interferograms remaining for that pixel. 3. If the standard deviation of a pixel's estimated velocity is greater than Υ, then the pixel is excluded from the rate map.
These three thresholds are designed to extract a coherent velocity for each pixel in the rate map, while excluding observations and pixels that are contaminated by residual unwrapping errors or atmospheric artifacts. In our analysis, Ψ is set as 10 interferograms, Ξ is set as 2, and Υ is set as 2 mm/yr. These values were chosen in order to maximize spatial coverage of the rate maps while also minimizing errors.
The time interval and area covered by our InSAR data includes no earthquake larger than M w 5.7, and no coseismic signals were observed in any of our interferograms. We assume that any postseismic deformation can be neglected due to the small size of these earthquakes. All the GPS data sets used in this study are either not affected by earthquakes larger than M w 5.5 or have been screened and corrected for any coseismic offsets.

10.1002/2013JB010909
is constant within each triangle. Therefore, the GPS velocity U or InSAR LOS velocity X LOS of any point within a triangle may be represented in terms of the velocities of its vertices: where and are longitude and latitude of the observation point, u m is the velocity of vertex m, Z LOS is the LOS vector joining the satellite to the point ( , ) on the ground, and N m is the interpolation or shape function. This shape function describes the contribution of u m to U or X LOS as a function of latitude and longitude: For m = 1, the shape function coefficients a m , b m , and c m are as follows: The equivalent shape functions for vertices 2 and 3 are found by cyclic permutation of the subscripts in the order 1, 2, and 3.
Combining equations (B1) and (B2), we get our combined observation equations: where d gps and d sar are observation vectors made up of GPS (U) and InSAR LOS (X LOS ) velocities respectively, A gps is the design matrix consisting of the shape functions for the GPS points, A sar is the design matrix composed of shape functions and local LOS unit vectors, A orb is the orbital design matrix to solve for any remaining orbital errors in the InSAR data, m vel is the model vector containing velocities for the vertices in the triangular mesh (u m ), m orb are the orbital parameters for each set of InSAR data, 2 is the smoothing factor that determines the strength of Laplacian smoothing, and ∇ 2 is a scale-dependent umbrella operator, a discretized approximation of the Laplacian smoothing operator given as [Desbrun et al., 1999] where u i and u j are velocities at vertices i and j, |e ij | is the length of the edge e ij joining the two vertices, and N 1 (i) is the set of vertices joined by only one edge to vertex i.
The InSAR data sets used in this study have already been corrected for orbital errors as they have effectively been tied to the Eurasian reference frame of the GPS data during rate map construction, so m orb only consists of a static offset for each rate map of relative LOS velocities. In this study the magnitude of the offset for each rate map is close to zero and is <0.4 mm/yr in all cases.
The observation equations can be solved using a least squares inversion, weighted by the full variance-covariance matrix for the observations, which is formed from the formal uncertainties of the GPS data and from the full covariance matrix of InSAR data modeled using an exponential function (equation (2)) and the maps of formal uncertainty that accompany each rate map. However, these formal uncertainties are unrealistically small at ∼1 mm/yr, while we previously estimated uncertainties of 1.20-1.85 mm/yr from the overlap of neighboring rate maps. We therefore scale these maps of formal uncertainty by dividing by their respective median values and multiplying by the uncertainties we estimated from the overlapping rate maps in section 4. This work was supported by the Natural Environment Research Council (NERC) through the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET) and Earthquakes without Frontiers (EwF) programs. All Envisat data were provided and are copyrighted by the European Space Agency. ECMWF ERA-Interim data were provided by the British Atmospheric Data Centre. We are grateful to JPL/Caltech for use of the ROI_PAC software. All figures were prepared using the GMT package [Wessel and Smith, 1998]. We thank John Elliott and Hua Wang for valuable discussions and Rob Reilinger and Kurt Feigl for constructive and thorough reviews. We also thank the Editor, Paul Tregoning, for his helpful comments and suggestions.