Coseismic Displacements from Moderate-Size Earthquakes Mapped by Sentinel-1 Differential Interferometry: The Case of February 2017 Gulpinar Earthquake Sequence (Biga Peninsula, Turkey)

We study the tectonic deformation from the February 2017 shallow earthquake sequence onshore Biga Peninsula (NW Turkey, NE Aegean region). We use InSAR interferograms (Sentinel-1 satellites) to identify the seismic fault (striking N110◦E) and seismological data (parametric data and Moment Tensor solutions from NOA and KOERI catalogues) so as to refine its geometry and kinematics using inversion techniques. Despite the moderate magnitudes of the main events of the sequence (5.0 ≤ Mw ≤ 5.2), the total surface deformation is 2.2 fringes (or maximum 6.2 cm along LOS) and it is well visible with InSAR because of the shallow depth of the four main events (6–8 km) and the good coherence of the signal phase. Our geodetic inversion showed that the fault has normal-slip kinematics, dimensions of 6 by 6 km (length, width) and dips at 45◦. The InSAR data are fitted by a uniform slip of 28 cm. In addition, 429 earthquakes were relocated with the HypoDD software and the use of a 1-D velocity model. The dip-direction of the fault is not retrievable from InSAR, but a south-dipping plane is clear from seismology and the aftershocks distribution. The spatial distribution of relocated events indicates the activation of one fault with a rupture zone length of about 10 km, a result of the occurrence of off-fault aftershocks along strike the main rupture. A stress inversion using 20 focal mechanisms (M ≥ 3.6; NOA solutions) indicates that faulting accommodates a N196◦E extension. It is confirmed that moderate (5.0 ≤ M ≤ 5.2) shallow events can be traced in InSAR studies and can produce surface displacements that provide useful data in fault inversion.


Introduction
In February 2017, a shallow seismic sequence of moderate events hit the region of Gulpinar-Babakale, onshore northwestern Anatolia, Turkey (Table 1, Figure 1). The four main events with magnitudes between 5.0 ≤ M w ≤ 5. Fourteen weaker events occurred between February 6 and 28 with Mw between 3.6 and 4.6, then the seismicity decreased quickly. Two events of Mw = 4.4 and 4.0, that occurred on January 14 and 15, 2017, might be seen retrospectively as foreshocks. The moment tensor inversion of 20 events (Table  1; NOA solutions) indicated primarily normal-slip kinematics. The main events occurred at depths 6-8 km and with dip angles ranging between 34° and 59°. Red triangles indicate the location of the broadband seismic stations that provided phase data for relocation of seismicity. Vectors indicate GPS velocities (red after [1]; green after [2]) represented with respect to stable Europe. NAT stands for North Aegean Trough. Thin black lines are traces of major faults from the NOAFaults database [3]. Red star indicates the NOA epicentre of the M6. 3 12 June 2017 earthquake south of Lesvos. Boxes with dashed boundaries indicate the frames of the Sentinel-1 satellite SAR images used (white-ascending, black-descending). Table 1. January to March 2017 NOA location and moment tensor parameters of earthquakes of Mw ≥ 3.6 in the Gulpinar-Babakale area (NW Turkey). The four main events, with Mw ≥ 5.0, are shown in bold (no 3, 4, 6 and 15 in column 1). Details on the MT solutions can be found online at http://bbnet.gein.noa.gr/HL/seismicity/mts/revised-moment-tensors.  Red triangles indicate the location of the broadband seismic stations that provided phase data for relocation of seismicity. Vectors indicate GPS velocities (red after [1]; green after [2]) represented with respect to stable Europe. NAT stands for North Aegean Trough. Thin black lines are traces of major faults from the NOAFaults database [3]. Red star indicates the NOA epicentre of the M6. 3 12 June 2017 earthquake south of Lesvos. Boxes with dashed boundaries indicate the frames of the Sentinel-1 satellite SAR images used (white-ascending, black-descending). The onshore location of the Gulpinar events points towards the activation of a seismic fault that runs about 10 km north of the Adramytion (Edremit) Bay branch of the North Anatolia Fault (NAF) near the island of Lesvos, where a branch of the NAF enters the Aegean Sea (also known as Skyros-Edremit fault; [2,[4][5][6]; Figure 1). The earthquake sequence occurred in the Biga Peninsula, which is mostly made of metamorphic rocks that are unconformably overlain by Eocene and Oligo-Miocene volcanic and volcano-sedimentary rocks [7,8]. The Adramytion (Edremit) branch of the NAF transcends the Biga Peninsula, before entering the Aegean Sea about 50 km to the east of the 2017 sequence ( Figure 1; [5]).
In the following, we analyze seismicity relocations and seismic moments from NOA (Table 1) together with Synthetic Aperture Radar interferograms (InSAR) made with images from the ESA Sentinel-1 satellites to constrain the parameters of the activated fault. Despite the moderate magnitudes, the surface deformation is well visible with InSAR which indicates that the fault is shallow. We then model the interferograms to retrieve the parameters of the normal fault and discuss them. The study of this seismic sequence offers important insights on the geometry and location of the seismic fault near Gulpinar in the context of its proximity to a major transcurrent structure.

InSAR
The ESA Sentinel-1A (S1A) and Sentinel-1B (S1B) satellites carry a C-band synthetic aperture radar of wavelength 56 mm and four operating modes. We used images of 5 × 20 m resolution acquired in the Interferometric Wide (IW) mode with a 250 km swath achieved using the Terrain Observation with Progressive Scans (TOPS) technique. In this technique, subsets of echoes, called bursts, are recorded sequentially along three sub-swaths (IW1, IW2 and IW3). The 6-day repeat time of S1 imagery leads to very coherent interferograms ( [33]).
We used six Single Look Complex (SLC) IW scenes acquired on the ascending tracks 29 and 131 and four on the descending track 36, covering the pre-, co-, and post-seismic periods ( Table 2). We calculated four ascending differential interferograms, two for track 131 and two for track 29, and five descending differential interferograms (Table 2, Figure 2). The ascending interferograms are shown in Figure 3 and the descending ones in Figure 4. The incidence angles for all scenes used together with the resulting line-of-sight (LOS) vectors are reported in Table 3.  Table 2. S1 interferograms produced to estimate the fault parameters. Acquisition times are in the form yyyy-mm-dd-hh-mm. The satellite tracks are shown in Figure 1.

Satellite (m-s) Track Master Image (m) Slave Image (s) Perpendicular Baseline (m)
Ascending S1A-S1B 131 201701311607 201702061606 72 S1A-S1B 131 201701311607 201702181606 −20 S1A-S1A 29 201702051615 201702171615 −8 S1A-S1B 29 201702051615 2017021151615 42 Descending S1B-S1A 36 201701310414 201702060414 6 S1A-S1B 36 201702060414 201702120414 96 S1A-S1A 36 201702060414 201702180414 −62 S1B-S1A 36 201702120414 201702180414 −159 S1B-S1A 36 201701310414 201702180414 −56 Table 3. Incidence angles for the Gulpinar area as computed from S1 orbital information. Heading is angle counterclockwise from North, subswath indicates the part of the S1 frame that was processed (frames are shown in Figure 1). LOS vector it is in unit of 1, east, north, Up components, respectively.   ). In the importing step, a radiometric calibration was applied, using the noise and calibration meta-data provided with the SLC data. The state vectors included in the meta-data are accurate enough to be  ). In the importing step, a radiometric calibration was applied, using the noise and calibration meta-data provided with the SLC data. The state vectors included in the meta-data are accurate enough to be  [34]). In the importing step, a radiometric calibration was applied, using the noise and calibration meta-data provided with the SLC data. The state vectors included in the meta-data are accurate enough to be used without further refinement. Overall, it is found that the quality of the state vectors provided in the S1 SLC data products is already high and efficient for co-registration, geocoding and interferometry. However, in the case of anomalous results, concatenated scenes or using advanced interferometric techniques (e.g., Persistent Scatter Interferometry), which is not the case in the present study, Precise Orbit Determination (POD) precision state vectors should be applied ( [34]). The co-registration step was carried out in order the slave scenes to obtain the same geometry as the master scene (Table 2). For TOPS interferometry, the accuracy of the co-registration is crucial. Specifically, in the azimuth direction an accuracy of a few thousandths of a pixel is required, otherwise phase jumps between subsequent bursts are observed. To ensure this very high co-registration accuracy, we used a method by exploiting the effects of the scene topography, with a matching procedure refined by a spectral diversity analysis that considers the interferometric phase of the burst overlap region ( [34]). The coherence and phase standard deviation thresholds used are 0.6 and 0.8, respectively. The procedure iterates the matching refinement until the azimuth correction is <0.01 pixel. After reaching this quality, it iterates the spectral diversity until the azimuth correction is <0.0005 pixel. The achieved accuracies of the co-registration are of the order of 1/100 azimuth pixel. To reduce the phase noise and obtain approximately square pixels, the co-registered scenes were multi-looked using 10:2 looks (range: azimuth) with pixel spacing of about 23 m × 27 m (range × azimuth). The final interferograms are formed by subtracting the topographic term estimated using the SRTM digital elevation model (DEM) at 30 m resolution and filtered ( [35]). Final products are geocoded with a final resolution of 30 m. Thanks to the short temporal separation, the coherence of the nine interferograms is very high. In some of them tropospheric effects can be seen but they do not mask the co-seismic deformation which is visible in all.

Orbit
We do not expect a significant variability of the tropospheric conditions due to turbulence or layered tropospheric delay, because (a) the spatial extent of the deformation pattern is small (almost 10 km × 10 km) with good coherence, (b) most of the study area is flat (the altitude variability is below 200 m), (c) the fringes' spatial extent is large, thus the spatial (X-Y) error of the sampling procedure (used in fault inversion) is small and so is the (relative) deformation error at the sampling points, and (d) the earthquakes occurred in winter, with a tropospheric effect usually smaller than other seasons. Finally, the pattern of the deformation is regular. We can only spot small irregularities on some of the interferograms; in particular at the NNE of the outer fringe of the deformation lobe. Figure 3a shows the deformation produced by the first two main events, Figure 3c the first three, and  This first analysis of those nine interferograms, all co-seismic, shows that the two stronger events, in terms of ground deformation signature, are the first (EQ1) and the third (EQ3). It shows also that the location of the deformation remains the same all along the sequence, suggesting that all events are embedded almost at the same place, with lateral variations, if any, indistinguishable in our interferograms. Moreover, for the first and the last event, we could produce interferograms with only one event (Figure 2), but for the two intermediate events we do not have interferograms corresponding to a single event; they contain at least two events. For this reason, we develop a minimization strategy to separate each single event using a wide number of different interferograms (those of Table 4). For this purpose, we downloaded forty-nine (49) interferograms of the area produced by the COMET LiCSAR service (http://comet.nerc.ac.uk/COMET-LiCS-portal/), 10 from the ascending track 29, 12 from the ascending track 131, 13 from the descending track 36, 14 from the descending track 109 (Table 4). Among them, 40 are co-seismic (with 1, 2, 3 or 4 main events captured), 3 are pre-seismic, and 6 are post-seismic.
For the fault inversion we used the LiCSAR data because the LiCSAR archive was complete with no need of any other processing. The LiCSAR interferograms are of very good quality and comparable to those we produced ourselves (Figures 3 and 4). The studied area is small and with good coherence, so all interferograms were clear enough to identify fringes (or absence of fringe) with no ambiguity.
We measured the number of fringes in each interferograms, and estimated, by L1 norm minimization the best fitting contribution of each of the four events that minimizes the overall scatter model-observations. For each interferogram, the earthquakes that contribute to the deformation are indicated by a Y in columns 6 to 9 of Table 4, while values in columns 4 and 5 are the observation (visual inspection of LiCSAR data), and model (sum of the 4 contributions), respectively. The accuracy of the measurement made on each interferogram is primarily governed by the tropospheric noise, and we estimate that this accuracy is ±0.2 fringes. Our minimization also implicitly assumes that all events occur roughly on the same plane (we wrote earlier that this is an assumption we make) according to the fact that we could not see any lateral fluctuation of fringes within the population of interferograms. Finally, we also assume here that the peak fringe value is dominated by the vertical displacement of the co-seismic deformation, and roughly equivalent for ascending and descending views. In the case of normal faults, theoretical modeling shows that this assumption is valid at better than the 5% level and robust in the near field-in practice, everywhere above the projection of the fault at the surface. Table 4. List of forty-nine (49) interferograms of the Gulpinar-Babakale area produced by the COMET LiCSAR (http://comet.nerc.ac.uk/COMET-LiCS-portal/). First two columns: dates of the master and slave image. Third column: track. Fourth column: number of fringes measured on the interferograms. Fifth column: best fitting value that minimizes the scatter with the observations, obtained by summing the contributions of each earthquake captured by the interferogram (indicated by a Y in columns 6 to 9; the numbers in column 5 indicate the best fitting deformation, in number of fringes along the line of sight) in the minimization.  Our modelling approach itself does not require phase unwrapping. The edges of the fringes are enough for sampling for the inversion procedure. The inversion has as output the model which can indicate the zero deformation zones (and predict the modelled deformation at the sampling points, wrapped and/or unwrapped) and thus the reference point.
The separation of the four sources made with those 49 interferograms leads to the contributions listed in Table 5. The relative contributions of the four earthquakes (Table 5) become: 33 ± 7% (first event), 16 ± 7% (second event), 37 ± 7% (third event), 14 ± 7% (fourth event). The 7% uncertainty on the percentage of contribution is determined based on the 0.2 fringes uncertainty used in input for the observations (read above) and the average scatter of the whole minimization 0.26 fringes). We also expect the tropospheric contribution to be minimized in the inversion procedure, due to the fact that we used in total 31 different epochs in the 49 analyzed interferograms, therefore we expect the tropospheric variability to be averaged by using such a large number of epochs. Also due to the small extent of the deformed area and modest topography, we estimate the tropospheric effect to be not more than 0.3 fringe for each single interferograms (after having visually inspected the interferograms) and the overall tropospheric effect to be below 0.1 fringe globally after the minimization process. If all events occur at the same depth, then the contribution of the fourth event is less that what is should be, and thus its moment (and M w = 5.1 magnitude) are overestimated. Conversely, if its moment (and magnitude) is correctly estimated, and consistent with the others, then our result suggests that this event of February 12 is deeper that the three previous events, possibly one or two km deeper.
This second analysis (using the LICSAR service) leads to a best fitting estimate of 2.2 fringes for the whole sequence, thus slightly less than the 2.6 fringes estimated in our initial analysis of our nine interferograms, but consistent. For the modelling in the following step, we will use 2.2 fringes, thus 6.2 ± 0.6 cm in the LOS, equivalent to 6.8 ± 0.6 cm of total subsidence for the entire sequence.

Seismology
The main shock of 6 February 2017, 03:51 UTC (M L = 5.0, M W = 5.2), as well as the aftershock sequence of 546 catalogue events, were relocated (with magnitude 1.0 < M L < 4.4; 88% between M L = 2.1 and 3.5; Figure 5a). The phase data were derived from manually picking the arrival time of the P-and S-waves from the recordings of 17 permanent seismological stations of the Hellenic Unified Seismic Network (HUSN) and from the Kandilli Observatory and Earthquake Research Institute (KOERI) Seismological Institute (see Figure 1 for station locations). Excluding one KOERI permanent station (GPRN), all other stations are located at least 40 km away from the Gulpinar sequence. The seismic datasets cover the period 1-28 February 2017.
Two location algorithms were used for the initial relocation procedure, namely the HYPO71 earthquake location program ( [36]) and the double-difference algorithm, HYPODD ( [37]). The last algorithm uses the catalog data from the seismic networks and the cross-correlation of the events. This method works by minimizing the double-difference between observed and calculated travel-times for pairs of neighboring events and relocates the hypocenters by assuming that if the inter-event distance is much smaller compared to the approximately common ray-path then travel-time differences can be attributed to the inter-event distance. The algorithm can incorporate both catalogue and cross-correlation differential times. The final relocation was performed by applying the double-difference algorithm, HYPODD, using the 1-D velocity model proposed for the study area by [38] as it was considered the more representative.
The first part of the work involved routine analysis of the recorded seismograms that were carried out in order to precisely pick P-and S-wave arrival times. For each phase, a quality value and a time uncertainty in seconds were assigned. A preliminary hypocenter location was performed using the HYPO71 code considering only events having more than three P-and three S-wave arrival-time readings. By removing time-biased data and employing a subset of events located with accuracy better than 1.5 km and having solution RMS < 0.2 s, the average Vp/Vs ratio and the velocity model of the area were investigated, using the mean travel-time residuals and location uncertainties (RMS, ERH, ERZ) minimization method. The use of the optimal, as previously described, dataset yielded satisfactory horizontal location uncertainties. More specifically, the majority of the epicentral solutions present ERX (98.6%) and ERY (97.4%) less than 1 km (Figure 5d,e). The solution RMS (root mean square error) of the events mainly ranges between 0.1 s and 0.4 s, with 65.8% of the values being smaller than or equal to 0.2 s (Figure 5c).
The above procedure allowed the relocation of 429 earthquakes, revealing, in some cases, significant differences in both the depths and the epicenters of these earthquakes with respect to the initial (HYPO71-HYPOINVERSE) catalogue (Figure 6a). Most of the hypocenters are clustered in the depth-range 5-12 km (Figure 5b). The dip-direction of the fault plane is clearly marked towards the south (Figure 6b) as it is inferred by the alignment of aftershock foci. In addition, the seismicity occurred at the same region where ground deformation was imaged by InSAR.

Stress Field
We analyzed the data set of Table 1 with the open source STRESSINVERSE package. This package runs a joint inversion for stress and fault orientations from focal mechanisms ( [39]), using the Michael's method ( [40,41]). The objective is to calculate the stress tensor that best explains the seismological observations. The inputs are the orientation and slip direction of a set of fault planes as determined from moment tensor inversion (Table 1)

Stress Field
We analyzed the data set of Table 1 with the open source STRESSINVERSE package. This package runs a joint inversion for stress and fault orientations from focal mechanisms ( [39]), using the Michael's method ( [40,41]). The objective is to calculate the stress tensor that best explains the seismological observations. The inputs are the orientation and slip direction of a set of fault planes as determined from moment tensor inversion (Table 1) Table 1 locations. Relocated events are from this study. A south-dipping plane can be seen on the vertical cross section C1-C2 that shows the foci projections orthogonal to fault's strike.  Table 1 locations. Relocated events are from this study. A south-dipping plane can be seen on the vertical cross section C1-C2 that shows the foci projections orthogonal to fault's strike.

Figure 7.
Map of northeast Aegean with published focal mechanisms by MT inversion (NOA solutions using the ISOLA code by [43]; beachball numbers correspond to events in Table 1). Thin black lines are active faults from the NOAFaults database [3].

Figure 7.
Map of northeast Aegean with published focal mechanisms by MT inversion (NOA solutions using the ISOLA code by [43]; beachball numbers correspond to events in Table 1). Thin black lines are active faults from the NOAFaults database [3].

Modelling the Fault Dislocation
The modelling is made for the entire sequence assuming pure normal faulting, and an elastic and homogenous medium. As we did for the separation of the four main sources in Section 2.1, we assume here that only the four main events contribute to the observed ground deformation. InSAR does not permit to discriminate between the south and north hypothesis for the fault plane, but the relocated seismicity (Figure 6b) clearly indicates a south-dipping plane. In input of the inversion, we assume a south-dipping plane, with a priori azimuth 110 • , dip-angle 45 • , rake −90 • , and 7 km centroid depth. We assume a uniform slip model because the Gulpinar earthquakes are relatively small (5.0 ≤ M ≤ 5.2) and the seismic fault did not reach the surface. This type of faulting leads to interferograms with a smooth and circular shape. With such shape we know (from our experience) that we can capture only the low spatial frequency component of the deformation. If variable slip exists, it is either small or not visible in the geodetic data. As a consequence, we consider that the safest and more robust approach was to invert for very simple source.
The inversion is made with the ascending interferogram 5-17/2/2017 (Figure 3d; track 29) and the descending interferogram 6-18/2/2017 (Figure 4d; track 109). The ascending interferogram covers the whole period of major seismicity and it is also, among the various produced interferograms, one with less tropospheric noise. The descending interferogram misses the first event that we determined earlier to account for 33% of the total (i.e., 0.8 fringes; Table 5). The best fitting fault has parameters listed in Table 6, and it is plotted in Figure 9 over ASTER topography (downloaded from https://asterweb.jpl. nasa.gov/gdem.asp). A uniform slip amount of 28 cm fits the observed surface deformation. We also constructed all four synthetic interferograms using our fault model (shown in Figure 10; one per viewing geometry) predicting 5.3-6.1 cm of ground motion away of the satellite. In Figure 10 we report the maximum and minimum motions in the line of sight (negative corresponds to increase of the LOS distance). Please note that the predicted positive values (i.e., motion towards the satellite) are very small and below our capacity of detection. The data are fitted at ±2 mm which is at the level of the noise. According to this model, the centroid is at 5.4 km depth, thus 0.5-1.5 km higher than the one estimated from seismicity. The total moment, assuming a crustal rigidity of 30 GPa, is 30.6 × 10 16 N m consistent with the 29.6 × 10 16 N m inferred from seismology (using the events listed in Table 1). This discrepancy between seismic and geodetic depths is common for shallow earthquakes. This is probably because the elasticity hypothesis is a coarse approximation in the uppermost crust (0-3 km); the very shallow (non-rigid) layers just "follow" the motion of the deeper layers, thus giving (in the modelling) the impression that the source is shallower ( [44,45]). In consequence, the effective centroid depth could well be around 6.5-7 km, and the fault tip at 4.5 km depth. (negative corresponds to increase of the LOS distance). Please note that the predicted positive values (i.e., motion towards the satellite) are very small and below our capacity of detection. The data are fitted at ±2 mm which is at the level of the noise. According to this model, the centroid is at 5.4 km depth, thus 0.5-1.5 km higher than the one estimated from seismicity. The total moment, assuming a crustal rigidity of 30 GPa, is 30.6 × 10 16 N m consistent with the 29.6 × 10 16 N m inferred from seismology (using the events listed in Table 1). This discrepancy between seismic and geodetic depths is common for shallow earthquakes. This is probably because the elasticity hypothesis is a coarse approximation in the uppermost crust (0-3 km); the very shallow (non-rigid) layers just "follow" the motion of the deeper layers, thus giving (in the modelling) the impression that the source is shallower ( [44,45]). In consequence, the effective centroid depth could well be around 6.5-7 km, and the fault tip at 4.5 km depth.    Figure 9; tracks are shown in Figure 1). The amount of LOS motion is provided as a negative number for motion away from the sensor or "subsidence" and positive for "uplift" (for example for descending track 36, the "subsidence" vs. "uplift" is 54 mm vs. 3 mm, respectively).

Tectonic Deformation and Moderate-Size Earthquakes
The 2017 Gulpinar-Babakale sequence is a very interesting case of shallow earthquakes with moderate magnitudes in the Aegean area. The tectonic subsidence mapped with InSAR reached about 2.2 fringes motion away from the satellite (LOS; Figures 9 and 10), moving from outside towards inside the 2nd fringe lobe. We mapped permanent deformation greater than 1 cm on the hanging wall block, spanning a N-S distance of 7 km (Figure 11). The maximum hanging wall subsidence is estimated between 5.4-6.1 cm (depending on viewing geometry; see red curve in Figure 11 for track 36). The strike-orthogonal cross section ( Figure 11) also showed that the deformation resulting from co-seismic motion along the south-dipping normal fault does not have a strong relation with present-day topography, because the highest co-seismic subsidence does not correlate with topographic lows. The depth to fault-tip is at 3-4 km ( Table 6), indicating that this E-W fault is a young structure accommodating crustal extension onshore Biga Peninsula.
The 2017 events demonstrate the importance of a systematic monitoring of the ground deformations over continents because moderate (4 < M < 6) events are capable to generate surface displacements, provided they occur at relatively shallow depths such as Gulpinar (5-8 km), Pawnee, Oklahoma (4-9 km rupture depth; USA; Mw = 5.8; [46]), La Habra, California (2 km centroid depth; USA; M = 5.1; [47]), Ain Temouchent (5 km centroid depth, Algeria, Mw = 5.7, [48]) or Little Skull Mountain (9-10 km; USA; Mw = 5.6; [49]). We note that a successful mapping of tectonic deformation with ENVISAT data was during the 2011 seismic swarm at southwest Peloponnese, where three shallow normal-slip events with magnitudes Mw = 4.7-4.8 occurred within three months near each other ( [50]). Moreover, the 2017 Ischia earthquake (M3.9) caused a subsidence of 4 cm in line of sight that was detected by InSAR because of its shallow depth, (800 m; [51]). At hypocentral depths below 10 km, the surface deformation due to moderate-size earthquakes is perhaps marginally detectable or undetectable by C-band InSAR (e.g., the 2013 Mw = 5.4 normal-slip event in central Greece; [52]). For moderate earthquakes deeper than 10-km depth, the focal mechanism of the earthquake is the  Figure 9; tracks are shown in Figure 1). The amount of LOS motion is provided as a negative number for motion away from the sensor or "subsidence" and positive for "uplift" (for example for descending track 36, the "subsidence" vs. "uplift" is 54 mm vs. 3 mm, respectively).

Tectonic Deformation and Moderate-Size Earthquakes
The 2017 Gulpinar-Babakale sequence is a very interesting case of shallow earthquakes with moderate magnitudes in the Aegean area. The tectonic subsidence mapped with InSAR reached about 2.2 fringes motion away from the satellite (LOS; Figures 9 and 10), moving from outside towards inside the 2nd fringe lobe. We mapped permanent deformation greater than 1 cm on the hanging wall block, spanning a N-S distance of 7 km (Figure 11). The maximum hanging wall subsidence is estimated between 5.4-6.1 cm (depending on viewing geometry; see red curve in Figure 11 for track 36). The strike-orthogonal cross section ( Figure 11) also showed that the deformation resulting from co-seismic motion along the south-dipping normal fault does not have a strong relation with present-day topography, because the highest co-seismic subsidence does not correlate with topographic lows. The depth to fault-tip is at 3-4 km ( Table 6), indicating that this E-W fault is a young structure accommodating crustal extension onshore Biga Peninsula.
The 2017 events demonstrate the importance of a systematic monitoring of the ground deformations over continents because moderate (4 < M < 6) events are capable to generate surface displacements, provided they occur at relatively shallow depths such as Gulpinar (5-8 km), Pawnee, Oklahoma (4-9 km rupture depth; USA; M w = 5.8; [46]), La Habra, California (2 km centroid depth; USA; M = 5.1; [47]), Ain Temouchent (5 km centroid depth, Algeria, M w = 5.7, [48]) or Little Skull Mountain (9-10 km; USA; M w = 5.6; [49]). We note that a successful mapping of tectonic deformation with ENVISAT data was during the 2011 seismic swarm at southwest Peloponnese, where three shallow normal-slip events with magnitudes M w = 4.7-4.8 occurred within three months near each other ( [50]). Moreover, the 2017 Ischia earthquake (M3.9) caused a subsidence of 4 cm in line of sight that was detected by InSAR because of its shallow depth, (800 m; [51]). At hypocentral depths below 10 km, the surface deformation due to moderate-size earthquakes is perhaps marginally detectable or undetectable by C-band InSAR (e.g., the 2013 M w = 5.4 normal-slip event in central Greece; [52]). For moderate earthquakes deeper than 10-km depth, the focal mechanism of the earthquake is the most important factor as dip-slip events produce more surface deformation than strike-slip events in the LOS direction of the satellite (considering simple, elastic Okada type fault models; e.g., [53]). Because of the cm-size surface deformation, its detectability is also influenced by environmental conditions as they may introduce phase decorrelation effects. It is well known that the output of SAR image processing also depends on local atmospheric conditions during image acquisition (dry vs humid) and land cover type (arid vs vegetated).
These moderate, shallow events provide crucial information to the tectonics/InSAR communities as to understand active tectonics, fault kinematics and fault interactions in Europe and Mediterranean as we cannot work only on large and relatively infrequent M > 6 events. The availability of frequent Sentinel interferograms enhanced our imaging (mainly geometry, slip distribution) of shallow seismogenic structures with no surface expression, but still capable of generating strong ground motions.
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 20 most important factor as dip-slip events produce more surface deformation than strike-slip events in the LOS direction of the satellite (considering simple, elastic Okada type fault models; e.g., [53]). Because of the cm-size surface deformation, its detectability is also influenced by environmental conditions as they may introduce phase decorrelation effects. It is well known that the output of SAR image processing also depends on local atmospheric conditions during image acquisition (dry vs humid) and land cover type (arid vs vegetated). These moderate, shallow events provide crucial information to the tectonics/InSAR communities as to understand active tectonics, fault kinematics and fault interactions in Europe and Mediterranean as we cannot work only on large and relatively infrequent M > 6 events. The availability of frequent Sentinel interferograms enhanced our imaging (mainly geometry, slip distribution) of shallow seismogenic structures with no surface expression, but still capable of generating strong ground motions.

Sentinels Observation Strategy
The 2017 Gulpinar study also validates the strategy of ESA and partners of producing and broadcasting systematic products processed at 100-m resolution. This resolution of 100 m is enough, because it is almost two orders of magnitude less than the extent of the deformed area.
Furthermore, this sequence of moderate events might be used as a good test case for tuning the software provided on the ESA GEP/TEP platform (https://tep.eo.esa.int/) or externally and compare them in a rigorous manner, and for addressing many of the basic and fundamental issues that are required for clean inter-comparisons and a better standardization of the InSAR solutions from various providers, such as the conventions for coordinates, data formats, pixel value for masking the sea, and DEM used for topographic correction.
A last point is that the Sentinel missions removed the major difficulty for InSAR analysis that was the poor temporal coverage of InSAR data, which made impossible in the past to distinguish deformation due to different earthquakes at the same location (e.g., [53]). Another interesting implication is the capability to distinguish deformation due to early M5+ aftershocks following a M6 Figure 11. Cross-section E1-E2 showing LOS motion (in mm; red line) orthogonal to fault's strike along with topography (in m) from ASTER sensor. Blue line indicates the displacement according to the average incidence angle. ASTER GDEM is a product of METI and NASA.

Sentinels Observation Strategy
The 2017 Gulpinar study also validates the strategy of ESA and partners of producing and broadcasting systematic products processed at 100-m resolution. This resolution of 100 m is enough, because it is almost two orders of magnitude less than the extent of the deformed area.
Furthermore, this sequence of moderate events might be used as a good test case for tuning the software provided on the ESA GEP/TEP platform (https://tep.eo.esa.int/) or externally and compare them in a rigorous manner, and for addressing many of the basic and fundamental issues that are required for clean inter-comparisons and a better standardization of the InSAR solutions from various providers, such as the conventions for coordinates, data formats, pixel value for masking the sea, and DEM used for topographic correction.
A last point is that the Sentinel missions removed the major difficulty for InSAR analysis that was the poor temporal coverage of InSAR data, which made impossible in the past to distinguish deformation due to different earthquakes at the same location (e.g., [53]). Another interesting implication is the capability to distinguish deformation due to early M5+ aftershocks following a M6 event (e.g., L'Aquila 2009 [54]; 2012 Emilia, Italy, seismic sequence, [55]), so that mainshock slip models are better constrained.

Fault Kinematics and Rupture Zone Dimensions
The mostly onshore location of the foci of the Gulpinar-Babakale earthquake sequence (Figure 6a,b), its normal-slip kinematics as it is indicated by the focal mechanism data (NOA solutions, see Figure 7) and the surface deformation (subsidence) mapped by InSAR (Figures 9 and 11) suggest that this sequence ruptured one south-dipping normal fault segment, although from seismic data alone it is not excluded a second slip plane. Assuming along-dip continuity in the N-S direction, it seems reasonable to identify the 2017 seismic fault as the down-dip plane of the Tuzla normal fault, a 10 km+ south-dipping normal fault outcropping near the town of Tuzla ( [56]). This could be confirmed in the future with local geological and geophysical data (incl. drawing geological cross-sections) and using our model parameters.
The spatial distribution of relocated events indicates an E-W orientation of the rupture zone of this sequence, in agreement with the focal mechanisms, with a total length of about 10 km (Figure 6b; along strike section D1-D2). The size (length) of the rupture zone is expected to be larger than the mainshock 6-km rupture extent (Table 6) because of off-fault aftershocks that immediately follow the mainshock, due to static stress transfer effects on elastic Earth ( [57][58][59][60][61]). In normal-slip events a large amount of static stress is transferred to optimally-oriented (receiver) faults near the mainshock, that are located along strike the source fault ( [60,62]) and host the off-fault aftershocks.
The 10 km length (L) of the 2017 rupture zone is in agreement with empirical estimates for earthquake magnitudes 5.0 ≤ M w ≤ 5.2 in the Mediterranean region (9 km; using the formula LogL = 0.47M w − 1.49 from [63]). A larger size than the empirical estimates is expected because of the occurrence of four (4) moderate-size events.
The 2017 sequence raised concerns over a probable triggering of a large earthquake on the neighboring NAF-branch (Figure 1), during the following weeks, which did not occur eventually due to the moderate size of the main events of the sequence.

1.
We used S1 interferograms to map tectonic deformation due to the February 2017 shallow earthquake sequence onshore Biga Peninsula (NW Turkey). The total surface deformation is 2.2 fringes (or maximum 6.2 cm along LOS) and it is well visible with InSAR because of the shallow depth of the four main events (6-8 km) and the good coherence of the signal phase. 2.
429 earthquakes were relocated with the HYPODD software. The dip-direction of the fault is not retrievable from InSAR, but a south-dipping plane is clear from seismology and the aftershocks distribution, as seen on a vertical cross-section across strike. The spatial distribution of relocated events indicates the activation of one normal fault with a rupture zone length of about 10 km, an effect of the occurrence of off-fault aftershocks. 3.
Using dislocation modeling, we identified that the seismic fault is striking N110 • E, has a dip angle of 45 • towards the south, it has dimensions 6 by 6 km (length − width), it slipped by an amount of 28 cm and it is located between 4-10 km depth.

4.
A stress inversion using 20 focal mechanisms (M w ≥ 3.6; NOA solutions) indicates that faulting accommodates a N196 • E extension.

5.
It is confirmed that moderate (5.0 ≤ M ≤ 5.2) shallow events can be traced in InSAR studies and can produce surface displacements that provide useful data in fault inversion.
Author Contributions: P.K. and I.P. did the InSAR analysis, P.B. and P.E. did the dislocation modeling and wrote part of the discussion, A.M. did the seismicity relocation and A.G. did the stress analysis, wrote the introduction and the discussion-conclusion sections.

Conflicts of Interest:
The authors declare no conflict of interest.