High‐Resolution Surface Velocities and Strain for Anatolia From Sentinel‐1 InSAR and GNSS Data

Measurements of present‐day surface deformation are essential for the assessment of long‐term seismic hazard. The European Space Agency's Sentinel‐1 satellites enable global, high‐resolution observation of crustal motion from Interferometric Synthetic Aperture Radar (InSAR). We have developed automated InSAR processing systems that exploit the first ~5 years of Sentinel‐1 data to measure surface motions for the ~800,000‐km2 Anatolian region. Our new 3‐D velocity and strain rate fields illuminate deformation patterns dominated by westward motion of Anatolia relative to Eurasia, localized strain accumulation along the North and East Anatolian Faults, and rapid vertical signals associated with anthropogenic activities and to a lesser extent extension across the grabens of western Anatolia. We show that automatically processed Sentinel‐1 InSAR data can characterize details of the velocity and strain rate fields with high resolution and accuracy over large regions. These results are important for assessing the relationship between strain accumulation and release in earthquakes.

well-instrumented regions such as California and Japan, the typical spacing between GNSS observation points of 10-50 km may still be insufficient to resolve strain localization at the scale necessary to distinguish between faults that are locked at the surface and those that are creeping aseismically (Elliott et al., 2016). The gaps in GNSS coverage are likely to persist and they have a major effect on the corresponding estimates of strain rate; regions of inferred high strain rate are controlled by the distribution of observations, potentially resulting in inaccuracies. Furthermore, temporal variations in strain accumulation around active faults may go undetected if velocities and strain rates are based on old or non-continuous observations (Bilham et al., 2016;Cetin et al., 2014;Rousset et al., 2016).
Interferometric Synthetic Aperture Radar (InSAR) provides spatially continuous measurements of surface motions, without instruments on the ground, with precision approaching that obtained from GNSS, and at a resolution that ranges from meters to hundreds of meters (e.g., Bürgmann et al., 2000;Hooper et al., 2012;Hussain et al., 2016;Walters et al., 2014;Wright et al., 2001). However, estimating interseismic strain remains challenging particularly in slowly deforming regions where ground displacements are small and error sources can dominate the differential radar phase (Elliott et al., 2016;Hooper et al., 2012;Shen et al., 2019). Recently, the number of InSAR-capable satellites and volume of associated data have increased and improvements in data quality and processing techniques now permit routine measurements of surface velocities over spatial scales appropriate for studying tectonic plate motions, regional fault systems, and the growth of mountains (e.g., Fattahi & Amelung, 2016;Grandin et al., 2012;Pagli et al., 2014;Tong et al., 2013;Wang et al., 2019;Wang & Wright, 2012). In particular, the European Commission's Sentinel-1 constellation, operated by the European Space Agency, with two near-polar orbiting synthetic aperture radar (SAR) instruments and a revisit period of 6-12 days for most active tectonic belts, has the potential to be a powerful hazard mapping and monitoring tool, which the geoscience community has begun to exploit (e.g., Elliott et al., 2015;González et al., 2015;Grandin et al., 2016;Shirzaei et al., 2017;Xu et al., 2020). By analyzing large volumes of short-revisit Sentinel-1 data, we can produce displacement time series with reduced impact from atmospheric noise.
In order to manage and process the large data volumes produced by Sentinel-1, we have developed opensource, automated workflows to efficiently produce interferograms and line-of-sight (LOS) time series and velocities (Lazecky et al., 2020;Morishita et al., 2020), which are valuable for a range of applications. Here we demonstrate our ability to measure large-scale interseismic deformation across Anatolia, an area encom-passing~800,000 km 2 and including the highly seismogenic North Anatolian Fault (NAF) Zone. We combine InSAR observations from the first~5 years of the Sentinel-1 mission with published GNSS data to create high-resolution surface velocity and strain rate fields for the region.

Sentinel-1 Data and LiCSAR Processing
We process Sentinel-1 SAR data acquired on 14 overlapping tracks (7 ascending and 7 descending) over Anatolia, which were selected to cover the entire region from the intersection of the NAF and East Anatolian Fault (EAF) in the east to the Aegean Sea in the west (Figures 1 and S1 in the supporting information [SI]). Sentinel-1 data are acquired on every 12-day revisit from the beginning of the Sentinel-1A operational mission in October 2014 and every 6 days since Sentinel-1B became fully operational in September 2016.
Our InSAR data set includes 40 spatially and temporally consistent frames (~250 × 250 km) that we define as part of the Sentinel-1 processing system LiCSAR (Figures 1 and S1) (González et al., 2016;Lazecky et al., 2020). By default, we construct temporal baseline interferograms to the six closest acquisitions in time (three forwards and three backwards) and ad hoc additional longer-time span interferograms to help deal with low coherence due to vegetation in summer months and snow cover in winter months. For each frame, this results in a network of~600-800 interferograms derived from~200 acquisitions ( Figure S2). Interferograms are downsampled (i.e., multilooked) by a factor of 20 in range and 4 in azimuth producing ground pixels of~80 × 80 m (resampled to~100-m spacing during geocoding), and the interferometric phase is unwrapped using a statistical-cost, network-flow algorithm (i.e., SNAPHU;Chen & Zebker, 2000. We partially mitigate atmospheric contributions to apparent displacement signals by applying the iterative troposphere decomposition model implemented in the Generic Atmospheric Correction Online Service for InSAR (GACOS) (

Geophysical Research Letters
reduces the interferogram phase standard deviations by 20-30% ( Figure S3) , which should reduce the uncertainty in our LOS velocities by a similar amount compared to the uncorrected velocities. Additional LiCSAR data processing details can be found in the SI.

Interseismic LOS Velocity Field Estimation and Uncertainties
We use LiCSBAS, an open-source InSAR time series analysis package integrated with the LiCSAR processing system , to derive InSAR LOS displacement time series and velocities. Our LiCSBAS workflow for Anatolia consists of further downsampling the data by a factor of 10 to a pixel size of~1 km, which is sufficient for large-scale tectonic applications. We perform statistical quality checks  prior to the small-baseline inversion, which yields incremental and cumulative displacements and the mean displacement velocity. Despite the short spatial and temporal baselines that generally characterize Sentinel-1 data, gaps in the small-baseline network may still be present due to severe decorrelation (e.g., due to snowfall), extended periods of time with no acquisitions, and after unwrapping consistency checks ( Figure S2). LiCSBAS circumvents this problem by imposing the constraint that displacements are linear in time (i.e., constant velocity) across the gaps (e.g., Doin et al., 2011;López-Quiroz et al., 2009). Finally, we estimate the uncertainty in the velocity from its standard deviation using the percentile bootstrap method (Efron & Tibshirani, 1986) ( Figure S4), and we mask pixels based on several noise indices ( Figure S5). We also test for potential velocity biases associated with short temporal baseline interferograms in a Sentinel-1 network (e.g., Ansari et al., 2020) by removing 6-and 12-day pairs for one LiCSAR frame prior to LiCSBAS velocity inversion ( Figure S11); the standard deviation of the difference between these results is small (~2 mm/year).
After LiCSAR/LiCSBAS processing, each frame has its own independent reference point for velocity determination (e.g., Figure S6). We transform the LOS rate maps into a Eurasia-fixed reference frame using a regional GNSS velocity compilation (Figure 1a and SI) following the method outlined in Hussain et al. (2018); for each frame, we estimate and remove the best fitting second-order polynomial between an interpolated, smoothed GNSS-derived horizontal velocity field projected into the satellite LOS and the InSAR velocities ( Figure 1 and SI). This transformation yields a velocity field where the longest wavelength signals are tied to the GNSS data, but it does not affect features at the~100-km length scale and below.
Fault-perpendicular profiles from the overlap zones of adjacent tracks provide an indication of how well the rate maps agree after the reference frame transformation (Figure 2). We present one profile taken from ascending-track data crossing the NAF near Ismetpasa and extending southward through the Konya Basin (Figures 1 and 2) and another taken farther east from descending-track data crossing the NAF and EAF. Both profiles show good agreement between adjacent frames and clear changes in LOS velocity across major fault zones, consistent with the localization of interseismic strain (Cavalié & Jónsson, 2014;Walters et al., 2014).
The bootstrap-derived uncertainties are generally considered to be underestimates, particularly if the network is not fully connected . Therefore, we also assess LOS velocity uncertainties by calculating the difference between our LOS velocities and a velocity field created by interpolating horizontal GNSS data (see SI), and an associated semivariogram γ at separation distances h ranging from 0 to 150 km for two off-fault frames ( Figure S6). Our ffiffiffiffiffiffiffiffiffi γ h ð Þ p values serve as an estimate of velocity uncertainty that is robust up to length scales of~150 km (see SI) (Bagnardi & Hooper, 2018). We use this approach to examine the evolution of uncertainty in our residual LOS measurements by estimating ffiffiffiffiffiffiffiffiffi γ h ð Þ p for progressively longer time intervals, and we find general consistency with the theoretical model derived for error analysis of GNSS time series data (Zhang et al., 1997) for the first~3 years of our Sentinel-1 time series ( Figure S6). At longer time intervals the uncertainty estimates on our Eurasia-fixed velocities reach a minimum of 2-3 mm/year, likely because our interpolated GNSS velocities are only accurate to this level, whereas the bootstrap-derived estimates continue to decrease with increasing time series length. However, this exercise is useful for determining our ability to measure small amounts of displacement, the time necessary to achieve a certain level of accuracy across different length scales , and how detection limits on interseismic velocities evolve with time ( Figure S6).  Figure S8). (c, d) Same as above but for a profile that crosses the NAF and EAF (Figure 1c).

Geophysical Research Letters
As an additional estimate of uncertainty, we also calculate the velocity residuals in the overlap areas for all frames. We do this by assuming horizontal motion only and by correcting for variable LOS by dividing the LOS velocities by the sine of the local incidence angles before multiplying by the sine of the incidence angle at the center of each track (e.g., Hussain et al., 2018;Walters et al., 2014). Histograms of the overlap residuals are approximately Gaussian with means close to zero and standard deviations of 3.1-3.7 mm/year ( Figure S7). Because LOS velocities are not purely horizontal and due to uncertainties in the GNSS velocities used to transform the LOS information into a Eurasia-fixed reference frame, these values can be considered upper-bound estimates of ffiffi ffi 2 p times the velocity uncertainties for the frames, giving an average LOS velocity standard deviation of~2.4 mm/year.

East-West and Vertical Surface Velocities for Anatolia
The Eurasia-fixed ascending and descending LOS velocities ( Figure 1) provide a detailed picture of Anatolian surface motions. The most prominent feature is the pronounced gradient in velocity across the NAF, from negligible motion north of the NAF to rapid westward motion of Anatolia relative to Eurasia south of the fault (e.g., Figure S6). Additional features include localized regions where there is apparent motion away from the satellite in both ascending and descending geometries indicating subsidence (Figures 1b, 1c, and S6).
To remove some of the ambiguity associated with LOS measurements, we follow the approach of Wright, Parsons, and Lu (2004) and decompose the LOS velocities into east-west and vertical components for pixels with both ascending and descending information: where V LOS is the Eurasia-fixed LOS velocity, θ is the local radar incidence angle, α is the azimuth of the satellite heading vector, and [V E V N V U ] T is a vector with the east, north, and vertical components of motion, respectively. This equation has three unknowns, and we have two observational constraints in the form of ascending and descending LOS velocities. To calculate the full 3-D velocity field, we note that both viewing geometries are relatively insensitive to north-south motion and use the interpolated, smoothed north-south component of the GNSS velocity field (Figure 3a) to constrain V N before solving for V E and V U . This approach does not result in smoothed east-west or vertical velocities because of the LOS north-south insensitivity.
The resulting decomposed east-west velocity field (Figure 3) is easier to interpret than the LOS rate map mosaic and shows large-scale westward motion of Anatolia at a rate of 20-25 mm/year relative to Eurasia, with visible strain (a localized velocity gradient) across the entire NAF and portions of the EAF (Figure 3b). Along-strike variations in the width of the velocity transition are also evident and correspond to portions of the NAF near Izmit and Ismetpasa, where shallow aseismic slip (i.e., creep) has been previously documented (Figures 1, 3, and S9) (Ambraseys, 1970;Bilham et al., 2016;Cakir et al., 2014;Hussain et al., 2016;Jolivet & Frank, 2020;Kaneko et al., 2013;Rousset et al., 2016).
The decomposed velocity field reveals that portions of Anatolia are experiencing rapid vertical motions. The clearest example is the large zone of subsidence with rates >50 mm/year surrounding the Konya Basin in south central Turkey (Figures 2a, 3c, 3e, S6, and S11), which is attributed to rapid aquifer compaction due to groundwater extraction (Caló et al., 2017;Üstün et al., 2015).

Velocity and Strain Rate Fields From Sentinel-1 InSAR and GNSS Data
To estimate rates of tectonic strain accumulation, we can calculate velocity gradients directly from the decomposed velocity field (see SI; Figure S10), but our preferred method (see SI for a detailed justification) involves combining InSAR LOS velocity maps with GNSS data and inverting for a velocity and strain rate model using the VELMAP approach (Wang & Wright, 2012) (see SI). The technique consists of dividing the study area into a mesh of arbitrary spherical triangles ( Figure S13), assuming the velocity varies linearly (i.e., the strain rate is constant) within each triangle, and using shape functions (England & Molnar, 2005) to solve for the unknown velocities at the vertices of each triangle using the observed InSAR and GNSS measurements. The associated strain and rotation rates are calculated using the spherical approximation equations of Savage et al. (2001). The inversion is regularized using Laplacian smoothing, the strength of which has an impact on the resulting strain rate magnitudes (Figures S13 and S15) including slightly underestimating the strain rates associated with active faults. The approach also does not allow for steps in the velocity field. Additional VELMAP modeling information can be found in the SI.
Comparison of our preferred Sentinel-1-and GNSS-based model with one based on GNSS data alone (see SI; Figures 4 and S15) reveals that the inclusion of InSAR data improves the accuracy of the velocity field 10.1029/2020GL087376

Geophysical Research Letters
( Figure S14) and better captures velocity gradients (and therefore also estimates of strain accumulation) along the major faults ( Figure 2). In the GNSS-only model, the second invariant of the horizontal components of the strain rate tensor (a measure of the total magnitude of the strain rate) indicates the NAF is characterized by a patchy distribution of regions straining at rates >100 nanostrain/year with even higher strain rates (≥150 nanostrain/year) primarily near clusters of GNSS sites around the western and eastern strands of the fault. Furthermore, central Anatolia is inferred to be essentially undeforming, but in Western Anatolia where earthquake focal mechanisms and the GPS-derived velocity and strain rate fields of Aktug et al. (2009) show that normal faulting and extension is prevalent, portions of the major grabens are straining at rates >50 nanostrain/year (Figures 3 and 4). In contrast, the combined InSAR and GNSS strain rate model shows spatially coherent strain rate magnitudes ≥150 nanostrain/year localized along nearly the entire length of the NAF. The previously identified creeping sections of the fault (Figure 2b) are also associated with elevated strain rates compared to the GNSS-only map, which exhibits high strain rates in the Izmit region (Figure 2d) but much lower rates near Ismetpasa ( Figure S9). (d-f) Strain rate components from a joint inversion of the GNSS and Sentinel-1 LOS velocities. Black and magenta bars in (c) and (f) represent the contractional and extensional principal strain rates, respectively. The maximum shear strain rates imply focused deformation along the NAF and EAF, and wholesale positive dilatation across western Anatolia is indicative of extension, whereas short-wavelength features in the dilatation field likely reflect anthropogenic vertical signals that result in subsurface expansion and contraction and contribute to noisy patches in the second invariant estimates. See Figures S15 and S19 for additional components of the strain rate tensor and a comparison with seismicity rates, respectively.

Geophysical Research Letters
For comparison, we also derive VELMAP strain rates using the alternative Global Strain Rate Model (GSRM) GNSS data set (see SI; Kreemer et al., 2014), which are characterized by localized patches of high strain along the NAF and in central Anatolia, largely controlled by GNSS site density ( Figure S17).
Another characteristic of our combined Sentinel-1 InSAR and GNSS result is that the inferred strain rates along the NAF (Figure 4) are typically half of those stemming from an analysis of Envisat InSAR data by Hussain et al. (2018), who took a different approach to estimating strain rate by modeling fault-parallel velocities using 1-D elastic dislocation theory. A main conclusion of Hussain et al. (2018) is that strain rates are essentially uniform along the entire length of the fault, implying that the interseismic strain rate is constant in time except in the first decade or two after a major earthquake. We attribute most of the strain rate magnitude discrepancy to the factor of 2 difference between shear strain rates obtained by computing the full strain rate tensor (Savage et al., 2001;Savage & Burford, 1973) and those obtained by taking the spatial derivative of the smoothed, decomposed east-west surface velocity field (i.e., the velocity gradient; Figure S10; see SI for a detailed explanation) or as in Hussain et al. (2018), the gradient of fault-parallel velocity profiles associated with slip on a dislocation in an elastic half-space. Once the factor of 2 is taken into account, our strain rate magnitudes are still lower than those of Hussain et al. (2018) but exhibit a similar first-order pattern, suggesting that the nearly constant along-strike strain rate is a real and robust feature of the NAF ( Figure S19). This result has important implications as it suggests that geodetic strain rate can be used as a long-term estimate of future seismic hazard independent of time since the last earthquake. Second-order differences in strain rate magnitudes are due to the smoothing implemented in VELMAP (Figures S13, S15, and S19) and not explicitly accounting for fault creep. For example, if we examine the NAF-parallel velocities in a profile that crosses the creeping zone near Ismetpasa, we see that our preferred solution does not capture the sharp velocity gradient evident in the GSRM GNSS velocities (Figure 2b). Rougher VELMAP models (e.g., Figures S15 and S16) better reproduce this gradient and return strain rate magnitudes more consistent with the dislocation-based estimates of Hussain et al. (2018) but also introduce unacceptably high levels of apparent noise in the central Anatolian strain field (e.g., Figure S15). Future efforts will focus on developing an improved approach to model regularization that includes spatially variable smoothing and accounts for fault creep.
While a velocity gradient across portions of the EAF is visible in the decomposed east-west velocities (Figure 3), our combined strain rate model infers relatively low levels of strain along this fault zone compared to the NAF (Figure 4), consistent with previous InSAR-based studies (e.g., Cavalié & Jónsson, 2014;Walters et al., 2014). Furthermore, we find appreciable, localized strain accumulation only along the northeastern half of the EAF that is not apparent in the GNSS-only model. This is also where the east-west velocity contrast is most apparent (Figure 3b). While there is some seismicity associated with the EAF (Figure 3a), the recently compiled 1900-2012 earthquake catalog for Turkey (Kadirioğlu et al., 2018) indicates that the associated magnitudes and thus total moment release are much lower than those along the NAF, supporting the notion that less strain is accumulating along the EAF than along the NAF (Bletery et al., 2020). The 24 January 2020 M w 6.7 Elaziğ earthquake (Melgar et al., 2020) occurred on the portion of the EAF where we resolve both an east-west velocity gradient and elevated strain rates on the order of~70 nanostrain/year (Figures 3 and 4). We infer maximum shear strain and dilatation rates ≥100 nanostrain/year associated with active grabens and normal faulting within a broad zone of positive dilatation across the Western Anatolian Extensional Province but relatively low levels of strain along the Central Anatolian Fault Zone (Figures 3 and 4).

Conclusions
We have produced, to our knowledge, the largest regional interseismic measurement from InSAR to date, covering a~800,000-km 2 area and the majority of Anatolia. Our strain rate model displays high strains along the major tectonic features, which is consistent with the distribution of seismicity (Figures 2a and S20). While the availability of abundant GNSS and Sentinel-1 InSAR data for Anatolia combined with favorable fault orientations make it ideal for such a study, our results demonstrate the potential of Sentinel-1 data for enhancing the picture of surface deformation and hazard in other regions. A key factor is the equal geographical coverage of Sentinel-1 ascending and descending data, which permits the retrieval of 2-D and 3-D deformation fields for tectonic zones globally even without the benefit of a dense GNSS data set (see SI; Figure S19). In addition, the relatively low uncertainties on Sentinel-1-derived interseismic velocities ( Figure S7) are beneficial for estimating strain across slowly deforming regions and for resolving small temporal changes in deformation throughout the earthquake cycle. Although some challenges still remain for fault systems where the majority of motion is in the north-south direction, Sentinel-1 represents a major improvement over past SAR data sets. This improvement is crucial for monitoring vertical motions from anthropogenic activities and for constraining earthquake hazard, particularly across regions with millennial earthquake recurrence intervals, where seismic hazard assessments based on incomplete historical earthquake records can dangerously underestimate the true hazard (Stein et al., 2012;Stevens & Avouac, 2016).

Data Availability Statement
All interferograms are available for download from comet.nerc.ac.uk/comet-lics-portal, the time series analysis software LiCSBAS can be accessed via Morishita et al. (2020), and information regarding accessing GACOS corrections can be found in Yu, Li, Penna, and Crippa (2018).