Space Geodetic Views on the 2021 Central Greece Earthquake Sequence: 2D Deformation Maps Decomposed From Multi-Track and Multi-Temporal Sentinel-1 InSAR Data

Pioneering efforts well studied the deformation decomposition of single earthquake using a pair of ascending (ASC) and descending (DES) track interferometric synthetic aperture radar (InSAR) data. However, deformation decomposition of sequent events is rarely discussed and hard to implement. That's because it's hard to ensure deformations related to each earthquake can be recorded by a pair of ASC and DES track data. Three sequent earthquakes (Mw>5.5) just hit Central Greece in March 2021, and this earthquake sequence provides us with a perfect case to study 2-D (east-west and up-down) deformation decomposition when the mentioned premise cannot be satisfied. In this context, we proposed a Multi-track and Multi-temporal 2-D (MTMT2-D) method. Its novelty and behind rationale are to decompose 2-D deformations of each event through fusing multitrack and multitemporal interferograms. Based on the decomposed deformations, we invert the slip distribution of three earthquakes respectively. We found that the decomposed deformations can better constrain the fault geometry than the single InSAR interferogram. Furthermore, our geodetic inversion results also suggest a domino-like triggering rupture process for this earthquake sequence. It indicates that our MTMT2-D method can potentially reveal more details about earthquake sequence.

. Previous researches reported that the well-known and -studied faults in Central Greece did not generate this earthquake sequence [1]. These earthquakes are likely to occur on unmapped or named as blind faults [2], [3]. Determining the fault geometry of these blind faults is the significant premise for studying earthquake physics and understanding regional seismogenic mechanism [4], [5]. Densely spatiotemporal interferometric synthetic aperture radar (InSAR) measurements lit a lamp for determining fault parameters in those regions only with sparse Global Navigation Satellite System data, but simultaneously determining the striking and dipping direction for an unmapped or named as blind fault are still challenging, especially for cases only with moderate magnitude (<Mw7) [6]. Moreover, InSAR technique can only capture one-dimensional (1-D) surface displacement along the line-of-sight (LOS) direction. When earthquakes do not rupture the surface, it is hard to identify fault trace and even the dipping direction from a "visually deceptive" 1-D displacement because reported cases revealed that different mechanisms can cause the same 1-D LOS displacements [7], [8]. The issues mentioned above can be potentially solved using multidimensional deformations [9], [10], [11], [12], [13]. Integrating ascending and descending track InSAR measurements with different viewing geometries can generate 2-D (east-west and up-down) deformations to avoid visual fraud in parameters determination [14]. However, there are still left clues in 2-D decomposition of earthquake sequence events. A sequence of earthquakes occurring on adjacent blind faults in a short time period can make determination more complicated. This is because an interferogram of a 6-day temporal baseline (like Sentinel-1A/B) may include co-seismic deformations associated with multiple earthquakes and even their early post-seismic deformations. The ideal case is that both ascending and descending tracks data are perfectly acquired between each earthquake event. But even for a polar-orbiting satellite constellation like Sentinel-1,  its repeat cycle is at least 6 days [15]. Exactly corresponding ascending and descending data are hard to achieve. The Central Greece earthquake sequence exactly presents this challenge. We can see from Table I and Fig. 1 in which all possible Sentinel-1 multitrack and multitemporal interferograms are listed, only the Mw5.6 event has a pair of ascending and descending interferograms only including the related event's co-seismic deformation. For another two events, the simultaneously existing ascending and descending data, which only corresponds to a single event, are absent. Pioneering research on the 2021 central Greece (Thessaly) earthquake sequence inverted 2-D deformations of total three earthquakes using two interferograms of T80 and T175 (Mar. 2, 2021-Mar. 14, 2021) [1]. The inverted 2-D deformation maps contain the cumulative co-seismic deformations of three earthquakes and early post-seismic deformations [3]. Subsidence dominates the maximum displacement in UD, and a small cumulative uplift locates around the subsidence. Similarly, the asperities of inverted cumulative EW deformation also overlap with each other. It is hard to determine the deformations caused by each earthquake. Thus, to isolate the deformation of three respective earthquakes, previous efforts only use interferograms containing deformations related to single event. Those interferograms containing multiple seismic deformations were discarded. In previous studies, it has been pointed out that there may be a domino-like triggering mechanism between the three earthquakes [16]. In addition, the Mw6.3 earthquake was accompanied by significant post-seismic deformations [1], and those deformations are likely to disturb the geodetic inversion of subsequent earthquakes. In this context, deformation maps of each earthquake which isolate the effects of other events is of great significance for further revealing the triggering mechanism and accurate slip inversion.
To this end, we proposed a MTMT2-D method which inherits the idea of the ordinary small baseline subsets technique (SBAS) [17]. Its novelty and behind rationale are to decompose 2-D deformations of each event through fusing multitrack and multitemporal interferograms. In the decomposition, MTMT2D include two significant steps: (Compressed sensing-based phase unwrapping error correction (CSPUEC) method [18] and strainmodel based InSAR for geo-hazards' monitoring Approach (SIGMA) [19].
CSPUEC and SIGMA are used to improve the accuracy of InSAR measurements without changing any SBAS workflow. Based on these good quality measurements, MTMT2D method can decompose 2-D deformations using the SBAS concept.
This article is organized as follows. Section II describes the tectonic settings in details. Section III briefly introduces the methodology of MTMT2D. We test its performance in Section IV. Discussions are given in Section V. Finally, conclusion is drawn in Section VI.  [20]. This fault zone consists of many parallel synthetic and antithetic fault segments with a roughly ENE-WSW trending direction of shortening and a dominantly NE-SW extensional stress. In consequence, central Greece is dominated by conjugate system of normal faults with general directions NW-SE and NE-SW [21], [22]. Moreover, seismological data indicate that strong earthquakes are mostly associated with the conjugate system [23]. From the large tectonic background [see Fig. 2(b)], the northward movement of the Arabian plate and the compression of the Eurasian plate resulted in active tectonics in Aegean Sea zone. Therefore, the northward motion of the Arabian plate pushes the smaller Anatolian plate westwards along the North Anatolian fault (NAF), continuing along the North Aegean trough region, which is the boundary between the Eurasian and south Aegean plates [24], [25], absorbing most compression of the northward motion. Normal faults with right-lateral shear motion associated with the NAF appears to become more distributed in north Greece boundary but transferred into N-S or NW-SE direction. Furthermore, the Holocene alluvial deposits in Upper Pleistocene accelerated the formation of distributed NW-SE trend basin in central Greece [26], which managed most of the earthquake events in this region. The occurrence of these three earthquakes may be triggered by the movement of NW-SE of regional tectonics.
This earthquake sequence has attracted extensive attention and discussions. The fault parameter solutions given by scholars from various institutions are shown in Table.S1 Ganas et al. [16] relocated the aftershocks and conducted field survey in the seismic area, and observed several surface SE-NW striking ruptures of about 5-10 km. Given the location transition of epicenters of three earthquakes, it was believed that the earthquake sequence may have occurred in the structural belt of the SSE-WNW direction. Papadopoulos et al. [27] derived the fault parameters using teleseismic P-waves and InSAR co-seismic deformations, and further proposed that this earthquake sequence was generated by three unknown normal faults that had not been activated before, and believed that the second event has a SW dipping direction, which is completely contrary to the results of Ganas et al. Novellies et al. [28] also used InSAR technology to studied this earthquake sequence, and their results confirmed the NE dipping direction of the second event. Yang et al. [3] studied this earthquake sequence using co-seismic and early postseismic InSAR data. Their results support the NE dipping direction of the second event. Different views on the fault dipping direction of the second event have brought controversy on the focal mechanism of the second event. That is the motivation of MTMT2-D to decompose deformations from multiple earthquakes, and to further study the fault geometry of the second event.

A. MTMT2D Step I: CSPUEC
Triplet phase closure check has proven its capability in PU error detection [29], [30]. For a triangle loop of three unwrapped interferograms ψ 1,2 , ψ 1,3 and ψ 2,3 obtained from three time nodes 1, 2, and 3, the corresponding phase closure Δψ 1,2,3 is defined as When the interferograms are contaminated by unwrapping errors, Δψ 1,2,3 will not be equal to 0. We can list the triplet phase closure of all triangles in the SBAS graph as Δψ. Then we can construct a linear equation to build the relationship between unwrapping errors to be corrected and closure phases. Its mathematical formulation can be described as where B is the incidence matrix generated by SBAS graph, round( Δψ 2π ) means the integer cycles of unwrapping errors, X represents the ambiguity cycles for each interferogram to be corrected, N is the number of triangle loops and M is the number of interferograms. In ordinary SBAS graphs, N is usually less than M. To illustrate, three interferograms can only construct one triangle loop. Therefore, the solution search for (2) is a rank-deficiency problem. To search for an optimal solution to (2) and ensure the solution are all integers, L-0 norm is the best objective function when we assume that part of interferograms contain no unwrapping errors. In this article, the solution search for (2) turns out to be a sparse signal recovery problem. L-0 is hard to approach because it is not a non-convex optimization problem. In our previous effort [31], we convert (1) into an integer linear programming (ILP) solution search, in which L-1 norm is used to approach the superiority of L-0 norm where X + M ×1 and X − M ×1 are natural numbers including zero. It can be seen from (3) that this reformulation of (2) limits the solution search range to natural numbers including zero instead of all integers. Therefore, the complexity of solution search has a sharp drop [32]. In this article, we use the subroutine "intlinprog.m" of MATLAB as a toolbox to solve this L-1 norm optimization problem.
Noted that round( Δψ 2π )here is to approximate unwrapping errors. In the generation of interferogram, multiple-looking and filtering operations can cause fading signal in triplet closure. Then, the triplet closure becomes the sum of the unwrapping errors and the fading signal. round(·) here assumes that the absolute value of fading signal does not exceed π. If it surpasses π, a pre-correction of fading signal is required [33], [34], [35], [36].

B. MTMT2D Step II: SIGMA-Based 2D Deformation Decomposition
Through the comparison of recorded event time from Fig. 1 and Table I, one can see that only T102 recorded independent co-seismic deformation of three earthquakes. Except for one interferogram in T7 that does not record any surface deformation (we excluded it in the below experiments), all other interferograms contain the surface deformation related to at least one earthquake. For each earthquake event, one ascending interferogram and one descending interferogram only including the associated deformation are unachievable. To recover surface deformation related to each earthquake, we borrow the idea of SBAS technique. For a better illustration, we give a mathematical description about its realization. This process can be formed by a linear equation of L = B ⊗ A · X shown as (4) shown at the bottom of the next page. where L j i is the jth interferogram in i track, α i and θ i are, respectively, the azimuth and incidence angle for i track, X i e and X i u each means the surface deformation associated with ith earthquake event in EW and UD direction, ⊗ is the Hadamard product, A is the incidence matrix which converts surface deformation into LOS displacement L in combination with projection matrix B. A is defined as (5) shown at the bottom of the next page, where T means matrix transpose operation.
Before getting (4) to work, there is still a left issue. That is the solution process of (4) is pixel-wise but the selected pixels of each track are not in unique positions. Strain model can be introduced to solve this issue. In its realization, a transformation matrix representing correlation between neighboring pixels is used to convert L = B ⊗ A · X which only supports for one pixel to neighboring pixel ensembles. That is (6) shown at the bottom of the next page, where Ω represents the neighboring m points in a defined distance threshold, G geo is the transformation matrix based on strain model, i means the ith pixel in Ω, and Δx i n , Δx i e , and Δx i u , respectively, represent the distance of ith pixel between the center pixel c. In this article, we set the distance range to 200 m. We used an iterative reweighting least squares (IRLS) method to complete it. In the first iteration, the initial weight is set to the unit matrix. In the next iteration, the weight matrix is updated by the previously estimated residuals. Once the subsequent solution is less than 10 −4 , or the number of iterations reaches 100, the iteration process will stop

C. Interferometric Processing and MTMT2D Workflow Summary
Four tracks of Sentinel-1 TOPS data covering Thessaly, Central Greece region are used for inverting 2-D deformation of the earthquake sequence [see Fig. 2 By following the standard TOPS interferometric processing, we applied coregistration toolbox for Sentinel-1 (Ct-Sent) software to generate all interferograms by means of our proposed minimum-spanning-tree enhanced spectral diversity (MST-ESD) method [37]. In MST-ESD, we first performed a geometrical co-registratn method [38] using external digital elevation model (DEM) and precise orbit to resample all images to a common reference, then we used MST method to select high signal-to-noise ratio interferogram pairs. For those robust pairs, we perform ESD technique to correct the residual azimuth mis-registration. Noted that when image number in a same track is less than 4, MST-ESD is degenerated to single-reference method. After co-registration tasks, we simulated the flat-earth phase using orbit and DEM and removed them from processed 14 SLCs. In this step, we also generated latitude, longitude, and incidence pixel-wisely. Then we generated all possible interferograms, 17 in all. We multilooked all interferograms with a factor of 1 by 4, respectively in the azimuth and range direction, resulting a grid resolution of ∼13m by ∼13m. 17 interferograms were then filtered by an adaptive Goldstein filter [39] and we masked out low coherence areas with a mean coherence less than an empirical value of 0.75. After pixel selection, we unwrapped all interferograms using all pair shortest path minimum cost flow method [31]. Then, MTMT2-D can be applied. In summary, MTMT2-D shown in Fig. 3 consists of two detailed steps: phase unwrapping error correction and constructing strain model and pixel-wise solution search. In the first step, triplet phase closure check is used to detect phase unwrapping errors occurring in redundant interferograms (1). Then, ILP is applied to solve the unwrapping errors and add back to all interferograms (3). In the second step, all ascending and descending tracks data are used to construct the strain model (6). The radius for neighboring point search is set to 200 m. IRLS is used to solve (6) for all points. Finally, EW and UD deformation maps for three earthquakes can be obtained.

A. 2-D Deformations of Three Sequent Earthquakes
After performing MTMT2-D, 2-D deformation maps of the three events were decomposed out as shown in Fig. 4. Due to the short time interval between event March 3 and March 4, part of the early post-seismic deformations of event March 3 [masked area in Fig. 4(b) and (d)] were captured in the solution process, which was beyond the co-seismic focus of this article, thus we   masked them out and did not further discuss it. Visual interpretation of the two-dimensional deformation field shows that the 2021 Central Greece earthquake sequence caused obvious eastward and downward displacements, which is consistent with the rupture mechanism of normal faults. Fig. 4(a) and (d) are the 2-D deformation maps of event March 3 Mw6.3 earthquake, in which the pattern of EW and UD deformation are basically identical, mainly distributed like an eye, with an area of about 15 km by 15 km. The maximum EW displacement is almost 20 cm, and the maximum subsidence displacement is 40 cm. Fig. 4(b) and (e) presents the 2-D deformations caused by the second event of Mw5.8 earthquake. The spatial distribution of the deformation is like those of the first event. The deformation pattern is dominated by subsidence, accompanied by a certain degree of eastward movements. One can observe obvious bull-eye-like deformations with a maximum EW displacement of 5 cm and a maximum subsidence displacement of 14 cm. Deformation pattern covers about a 10 km by 10 km area. Slightly different from the first event, the second event did not rupture to the surface, and no obvious fault traces can be found in the field investigation. Therefore, it is believed that the second event was caused by an unmapped blind fault in this region. The fault geometry and focal mechanisms of this blind fault will be further discussed. The March 12 event also did not rupture the surface, and there are no mapped faults in the deformation region. Visual inspection on the EW deformation map, we can see an obvious westward movement and an eastward movement. These two parts present a symmetrical butterfly shape, with the maximum magnitude of 3 cm. Deformations of these two wings show a characteristic of squeezing each other. The UD deformation is relatively simple, only shapes as subsidence. Its maximum magnitude is 8cm. Although the EW deformation is characterized by self-compression, the deformation pattern is still dominated by subsidence, so it can be inferred that this event can be explained by a normal fault rupture mechanism accompanied by a relatively small strike-slip component.

B. RMSE of Deformation Map
Quantitively, we calculate a RMSE map after solving (6). MTMT2D combines the advantages of unwrapping error correction and strain model to achieve a good accuracy of deformation estimation. The RMSE map is shown as Fig. 5. The overall root mean square error is shown in the subgraph, with an average magnitude of about 5 mm. The larger magnitude area is where the first earthquake ruptures the surface. Those large RMSEs may be due to the fact that the deformation gradient in this region is so large that even CSPUEC cannot completely correct it. Overall, the RMSE map reflects the deformation estimation accuracy of MTMT2D. Regarding to the large RMSE area, we have a further discussion on it as described in Section V-A.

C. Validations Using Burst Overlap Interferometry
Although we lack external validation data such as GPS, along-track deformation is used as validation data in this article. The rationale behind is that we treat the decomposed 2-D deformation maps as the constraint of geodetic inversion, and we forward model the 3-D deformation using the inverted slip models.
We presented the derived slip model and along-track deformation results in Fig. 6. In Fig. 6(a), we overlay the extracted along-track deformation map on the inverted slip model, and in Fig. 6(b) we show the forwarded along-track deformation over the slip model. It can be seen from Fig. 6(c) that their difference is little. We can therefore conclude that the 2-D deformation decomposed using MTMT2-D is correct. MTMT2-D is effective in 2-D deformation decomposition.

A. Significance of CSPUEC on Deformation Derivation and Slip Distribution Uncertainty Reduction
The effect of the unwrapping errors on the interferogram is that a 2π cycle can introduce a deformation estimation error of ∼2.8 cm (for C band), which is especially frequent in the near field where the deformation gradient is large. In this seismic sequence, there are obvious unwrapping errors in the interferograms Fig. 1(g) and (h). If CSPUEC is not performed, the final deformation results and even the geophysical results will be affected. To illustrate, we presented three cross sections on these three interferograms to show the significance of CSPUEC. Interferograms after CSPUEC (pink dots in Fig. 7(a)-(c) shows more continuous profiles than those results without CSPUEC (gray squares in Fig. 7(a)-(c). Their difference values (purple lines in Fig. 7(a)-(c) are all around multiples of 2.8 cm. It means that CSPUEC can compensate for these unwrapping errors.
In order to further explore the effects of the unwrapping errors on the geodetic inversion results, we adopted a Jack-Knife (JK) statistical method to calculate the slip uncertainty of two models constrained by 2-D deformation maps without and with  correction [40]. The imperfection of Green's function is caused by simplified fault model and other factors. These regions in the fault model are very sensitive to potential errors in the observations. To further investigate which region of the derived fault model is sensitive to the unwrapping errors, JK tests are performed. Larger slip uncertainty indicates larger observation errors. To only evaluate the effect of the unwrapping errors on it, we did not include the strain model in 2-D deformation inversion. We adaptively down-sampled the 2-D deformation maps without and with correction into irregular triangles [41]. The difference of 2-D down-sampled deformation maps is presented in Fig.S4. We repeat the geodetic inversion process 100 times, adding a white noise with an average of zero and a standard deviation of 5 cm to down-sampled deformation maps in each iteration. We calculate the standard deviation of the results 100 times and take it as the uncertainty of the model (see Fig. 8). We calculated the difference [see Fig. 8(c)] between the uncertainty of uncorrected results [see Fig. 8(a)] and corrected results [see Fig. 8(b)]. It can be seen from Fig. 8(c) that there is an overall uncertainty reduction for results with correction. The largest uncertainty reduction is about 0.15 cm. It is verified that unwrapping errors can disturb the inverted slip model and unwrapping error correction is helpful to mitigate its disturbance.

B. Is March 4 Earthquake North Dipping or South Dipping?
In a recent report on Central Greece earthquake sequence [42], a question has been raised about the dipping direction of March 4 earthquake. It is assumed that the seismogenic fault of this earthquake has two kinds of potential faults: north dipping and south dipping. This question arises from the fact that both faults can interpret the deformation pattern in interferogram (T102,  The analysis of 2-D deformation maps in Fig. 4(b) and (e) show that the March 4 earthquake caused significant eastward displacements and subsidence with a maximum value of ∼10 cm, and slight subsidence is shown in the southwest and northeast of the uplift region. Part of the westward displacements were also displayed in the northeast of the deformation fields, with a maximum of ∼4 cm. The EW displacement maps are dominated by the eastward deformation which can be interpreted as the right-lateral movement of the north dipping blind fault or the left-lateral movement of the south dipping blind fault. The vertical displacements are dominated by subsidence, so it supported that the seismogenic fault is normal fault. Therefore, the dipping direction the blind fault cannot be determined only from the visual inspection on 2-D deformation maps. A geodetic inversion analysis is further needed.
In geodetic inversion, 1-D interferogram and 2-D deformation maps are used as constraints, respectively. It is worth noting that, in order to avoid the influence of the post-seismic deformations of March 3 earthquake, we extracted out the major co-seismic deformation field of the March 4 earthquake using an ellipse. Two kinds of rectangular dislocation models of north dipping and south dipping are constructed based on elastic half space. The faults are divided into 700 m by 700 m rectangles along the striking and down dipping direction. Fig. 9(a) and (b) presents the inversion results of 1-D interferogram. It can be seen from the correlation between observation and modeled deformation that the inversion results of two potential faults are good. The lowest data-model correlation also reached 0.92. Their related slip models are presented in Fig. 9(e) and (f).
We also inverted two slip models based on 2-D deformations. Their two models based on two different dipping directions are shown in Fig. 9(g) and (h), respectively. The model [see    discrepancy indicates that the geodetic inversion results more support the north dipping hypothesis. It also reflects that the 1-D interferogram is not only visually deceptive but can even cause completely different inversion results. Therefore, 2-D deformations are better model constraints. For this case which only a pair of ascending and descending interferograms associated with one earthquake cannot be obtained, our proposed Pseudo SBAS method is significant to determining the fault parameters.

C. Major Mw6.3 Event Potentially Triggers the Subsequent Two Events
Given that we obtained 2-D deformations related to three sequent earthquakes using Pseudo SBAS method, three slip models related to three sequent earthquakes (see Fig.S1-3) can be obtained, making it possible to calculate coulomb stress change for each earthquake. We calculate the static coulomb stress change at 5km depth based on these three models [43] with a Young's modulus of 80 Gpa and a friction coefficient of 0.8. Co-seismic slip model of Mw6.3 earthquake [see Fig. 6(a)] is utilized to calculate the quasi-static stress change. The first Mw6.3 earthquake increases the coulomb stress in the northwest and northeast segments of the fault, of which northwestern region is the potential locations for future two earthquakes. We deducted that the rupture was facilitated by the March 3 earthquake at the fault locations of the Mw5.8 earthquake that occurred on March 4. Coulomb stress change calculated by co-seismic slip model of the March 4 earthquake also indicates an increase of the coulomb stress in the northwest of the fault where are the potential locations for March 12 earthquake (green contour lines in Fig. 10). The coulomb stress change results potentially can explain the cascading trigger mechanism of the three sequent earthquakes from the southeast to the northwest. The triggering mechanism may be affected by the rupture of kinematics of the first Mw6.3 event, a large slip asperity (orange contour lines in Fig. 10) located close to fault 2 potentially leads to the rupture of the adjacent fault 2 (March 4), also further causing the rupture of fault 3 (March 12) northwest to the fault 2. During the rupture of fault 2, the stress accumulation of fault 3 is aggravated and accelerated the occurrence of March 12 earthquake. The edges of the three slip asperities show a complementary pattern. Also, aftershocks (gray dots in Fig. 10) that occurred between the three earthquakes were distributed around the large slip areas. It indicates there could be a potential domino cascade rupture kinematic process.

VI. CONCLUSION
In this article, a MTMT2D deformation decomposition method has been presented which incorporates CSPUEC and SIGMA methods into the multi-track and multi-temporal SBAS workflow. The proposed MTMT2D method can accurate decompose 2-D deformation related to single event from redundant interferograms which include the deformation related to one or more earthquakes. The effectiveness of the proposed MTMT2-D workflow has been validated through a real dataset of Sentinel-1 TOPS interferograms over central Greece. In geodetic inversion, we find that the decomposed deformations can better constrain the fault geometry than a single interferogram. In addition, our geodetic inversion results also show that this earthquake sequence has a domino-like trigger rupture process. This shows that it is possible to reveal more details of seismic sequence events using our MTMT2-D method.
It is important to note that MTMT2D method ignores northsouth deformations. This hypothesis is feasible in 2021 Central Greece earthquake sequence because normal faults of this sequence are all in a roughly east-west dipping direction, but it is not applicable for earthquake cases where north-south deformation may dominate. To handle with this problem, more satellite data with different viewing geometries are required. NISAR and Sentinel-1C are highly expected in the near future.