Seismic and Geodetic Imaging (DInSAR) Investigation of the March 2021 Strong Earthquake Sequence in Thessaly, Central Greece

: Three strong earthquakes ruptured the northwest Thessaly area, Central Greece, on the 3, 4 and 12 March 2021. Since the area did not rupture by strong earthquakes in the instrumental period of seismicity, it is of great interest to understand the seismotectonics and source properties of these earthquakes. We combined relocated hypocenters, inversions of teleseismic P-waveforms and of InSAR data, and moment tensor solutions to produce three fault models. The ﬁrst shock (M w = 6.3) occurred in a fault segment of strike 314 ◦ and dip NE41 ◦ . It caused surface subsidence − 40 cm and seismic slip 1.2–1.5 m at depth ~10 km. The second earthquake (M w = 6.2) occurred to the NW on an antithetic subparallel fault segment (strike 123 ◦ , dip SW44 ◦ ). Seismic slip of 1.2 m occurred at depth of ~7 km, while surface subsidence − 10 cm was determined. Possibly the same fault was ruptured further to the NW on 12 March (M w = 5.7, strike 112 ◦ , dip SSW42 ◦ ) that caused ground subsidence − 5 cm and seismic slip of 1.0 m at depth ~10 km. We concluded that three blind, unknown and unmapped so far normal fault segments were activated, the entire system of which forms a graben-like structure in the area of northwest Thessaly.


Introduction
During March 2021 the area of northwest Thessaly, Central Greece, was ruptured by a sequence of three strong earthquakes [1] (Figure 1) that occurred on the 3 (10:16:08 UTC), 4 (18:38:19 UTC) and 12 (12:57:50 UTC) of the month. Source parameters of these earthquakes, determined by the National Observatory of Athens (NOA) as well as from other national and international seismological centers, are listed in Table S1. It is noteworthy that for the first and third strong earthquakes consistent moment magnitudes, M w~6 .3 and M w~5 .7, respectively, have been determined by the various centers. However, magnitudes determined for the second strong earthquake vary from M w = 5.9 to M w = 6.3, an issue examined further later. On the other hand, the various fault-plane solutions produced by seismological centers are consistent in that the three strong earthquakes shared focal mechanisms of normal faulting striking roughly NW-SE. This is also consistent with the active tectonics of the area, which is dominated by a NE-SW to N-S active extensional field [2][3][4][5][6][7][8]. The main faults that dominate the area are the Tyrnavos and Larissa Faults, which dip to N-NE, and the Rodia Fault, which dips to S-SW Figure 1. The Thessaly area, Central Greece, that ruptured by three strong earthquakes on 3, 4 and 12 March 2021 with magnitudes, Mw, of 6.3, 6.2 and 5.7, respectively, as determined in this study; stars show strong earthquake epicenters relocated in the present study. Blue and purple colored circles illustrate foreshocks and aftershocks occurring from 28 February to 3 March and from 3 to 15 March, respectively (epicentral determinations by NOA, http://www. gein.noa.gr/en/seismicity/earthquake-catalogs, last access 25 May 2021). Black lines and black triangles show the main normal faults in the area [2][3][4][5][6][7][8] and settlements reported in the main text, respectively; fault name codes: LF = Larissa Fault, TF = Tyrnavos Fault, RF = Rodia Fault. In the inset map, arrows illustrate directions of main lithospheric motions; thick and thin rectangles show the study area and the east Corinth Gulf area, respectively, the last one discussed in the main text.
The March 2021 earthquake sequence is quite challenging for a number of reasons. First, no strong earthquake has occurred in the entire Thessaly area since 1980, therefore, the last sequence is the first recorded by modern seismograph and other networks, e.g., CGPS. Besides, the March 2021 earthquakes ruptured an area which has been considered as a long-term seismic gap identified on the basis of the historical seismicity of the Thessaly area. The discussion on this issue opened after the observation that the north side of Thessaly, including the area of the March 2021 activity, has not ruptured by strong earthquakes since AD 1781, while the south side ruptured by several strong earthquakes from 1954 up to 1980, after an apparent quiescence since AD 1773 [9,10] (Figure 2). The main conclusion of that apparent space-time earthquake clustering was that should the seismicity pattern in south Thessaly repeat in the north side, then a migration of the activity would be expected northwards in the near future [9,10]. This important issue was Figure 1. The Thessaly area, Central Greece, that ruptured by three strong earthquakes on 3, 4 and 12 March 2021 with magnitudes, M w , of 6.3, 6.2 and 5.7, respectively, as determined in this study; stars show strong earthquake epicenters relocated in the present study. Blue and purple colored circles illustrate foreshocks and aftershocks occurring from 28 February to 3 March and from 3 to 15 March, respectively (epicentral determinations by NOA, http://www.gein.noa.gr/en/ seismicity/earthquake-catalogs, last access 25 May 2021). Black lines and black triangles show the main normal faults in the area [2][3][4][5][6][7][8] and settlements reported in the main text, respectively; fault name codes: LF = Larissa Fault, TF = Tyrnavos Fault, RF = Rodia Fault. In the inset map, arrows illustrate directions of main lithospheric motions; thick and thin rectangles show the study area and the east Corinth Gulf area, respectively, the last one discussed in the main text.
The March 2021 earthquake sequence is quite challenging for a number of reasons. First, no strong earthquake has occurred in the entire Thessaly area since 1980, therefore, the last sequence is the first recorded by modern seismograph and other networks, e.g., CGPS. Besides, the March 2021 earthquakes ruptured an area which has been considered as a long-term seismic gap identified on the basis of the historical seismicity of the Thessaly area. The discussion on this issue opened after the observation that the north side of Thessaly, including the area of the March 2021 activity, has not ruptured by strong earthquakes since AD 1781, while the south side ruptured by several strong earthquakes from 1954 up to 1980, after an apparent quiescence since AD 1773 [9,10] (Figure 2). The main conclusion of that apparent space-time earthquake clustering was that should the seismicity pattern in south Thessaly repeat in the north side, then a migration of the activity would be expected northwards in the near future [9,10]. This important issue was further examined and debated in several subsequent studies based on geological data and on space-time seismicity patterns [11][12][13][14][15]. further examined and debated in several subsequent studies based on geological data and on space-time seismicity patterns [11][12][13][14][15]. On the other hand, a major challenge is the investigation of the seismic fault(s) associated with the three strong earthquakes. The area extended to the south and southeast from the source of the first earthquake is dominated by mapped and well-studied faults, such as the Tyrnavos (TF) and Larissa (LF) ones, e.g., [1,2,4,5] (Figure 1). Α first attempt to investigate the seismotectonics of the three earthquakes was made on the basis of InSAR analysis [16]. The authors of that study suggested that all the three earthquakes have been associated with a normal fault system striking about SE-NW and dipping NE. As we will see later, we propose a different seismotectonic interpretation. The only surface-fault trace observed in association with the March 2021 activity was found in the area of Mesochori (Figure 1) (Prof. I. Koukouvelas, personal communication) but a relevant study is still in progress. Abundant but rather minor ground fissures were reported in the area [17]. Surface manifestations of soil liquefaction were also observed in several spots ( Figure 3) lying within the limiting distance, R, from the epicenter of the first earthquake, predicted by empirical relationships between R and earthquake magnitude, M [18]. However, such ground fissures do not provide evidence of surface-fault trace since liquefaction is controlled by other factors, such as M, R and the soil susceptibility. On the other hand, a major challenge is the investigation of the seismic fault(s) associated with the three strong earthquakes. The area extended to the south and southeast from the source of the first earthquake is dominated by mapped and well-studied faults, such as the Tyrnavos (TF) and Larissa (LF) ones, e.g., [1,2,4,5] (Figure 1). A first attempt to investigate the seismotectonics of the three earthquakes was made on the basis of InSAR analysis [16]. The authors of that study suggested that all the three earthquakes have been associated with a normal fault system striking about SE-NW and dipping NE. As we will see later, we propose a different seismotectonic interpretation. The only surface-fault trace observed in association with the March 2021 activity was found in the area of Mesochori ( Figure 1) (Prof. I. Koukouvelas, personal communication) but a relevant study is still in progress. Abundant but rather minor ground fissures were reported in the area [17]. Surface manifestations of soil liquefaction were also observed in several spots ( Figure 3) lying within the limiting distance, R, from the epicenter of the first earthquake, predicted by empirical relationships between R and earthquake magnitude, M [18]. However, such ground fissures do not provide evidence of surface-fault trace since liquefaction is controlled by other factors, such as M, R and the soil susceptibility. The present paper focuses on the investigation of the seismotectonics of the strong earthquake activity of March 2021 and the investigation of fault models for the three strongest earthquakes on the basis of seismological and satellite geodesy (DInSAR) observations. We are also interested to understand whether or not the strong earthquakes ruptured segments of the same fault or of different faults. Of equal interest is the investigation of possible links between the March 2021 seismogenic faults with the already known and mapped faults having surface expression in the area ( Figure 1). Results on such issues are of importance for the better assessment of the level of seismic hazard in the Thessaly area.

Materials and Methods
To investigate the seismotectonics of the three strong earthquakes that ruptured the Thessaly area during March 2021, we relocated their hypocenters, inverted teleseismic P-wave records, produced Interferometric Synthetic Aperture Radar (InSAR) images and inverted InSAR data, and took also critically into account moment tensor solutions produced by several seismological centers. Eventually we redetermined moment magnitudes and concluded with fault models for the three strong earthquakes and proposed a seismotectonic interpretation for the entire sequence.

Εarthquake Magnitudes
Consistent moment magnitudes, Mw, have been determined by various seismological centers for each one of the first and third strong earthquakes (Table S1). However, Mw determined for the second earthquake ranges from 5.9 to 6.3. It is remarkable that the worldwide recognized Global Centroid Moment Tensor (GCMT) Project and the GEO-SCOPE Observatory-French Global Network of broad band seismic stations did not produce moment tensor solutions for this earthquake. In addition, the United States Geological Survey (USGS), another renowned institute, has been able to determine only body-wave magnitude mb = 5.8 for the same earthquake. On the other hand, the solution produced by the GEOFON system of the Deutsche GeoForschungsZentrum (GFZ), Potsdam, provided Mw = 6.3 based on regional stations of the GEOFON system without the use of teleseismic records. The reason is that teleseismic records of the second The present paper focuses on the investigation of the seismotectonics of the strong earthquake activity of March 2021 and the investigation of fault models for the three strongest earthquakes on the basis of seismological and satellite geodesy (DInSAR) observations. We are also interested to understand whether or not the strong earthquakes ruptured segments of the same fault or of different faults. Of equal interest is the investigation of possible links between the March 2021 seismogenic faults with the already known and mapped faults having surface expression in the area ( Figure 1). Results on such issues are of importance for the better assessment of the level of seismic hazard in the Thessaly area.

Materials and Methods
To investigate the seismotectonics of the three strong earthquakes that ruptured the Thessaly area during March 2021, we relocated their hypocenters, inverted teleseismic P-wave records, produced Interferometric Synthetic Aperture Radar (InSAR) images and inverted InSAR data, and took also critically into account moment tensor solutions produced by several seismological centers. Eventually we redetermined moment magnitudes and concluded with fault models for the three strong earthquakes and proposed a seismotectonic interpretation for the entire sequence.

Earthquake Magnitudes
Consistent moment magnitudes, M w , have been determined by various seismological centers for each one of the first and third strong earthquakes (Table S1). However, M w determined for the second earthquake ranges from 5.9 to 6.3. It is remarkable that the worldwide recognized Global Centroid Moment Tensor (GCMT) Project and the GEOSCOPE Observatory-French Global Network of broad band seismic stations did not produce moment tensor solutions for this earthquake. In addition, the United States Geological Survey (USGS), another renowned institute, has been able to determine only body-wave magnitude m b = 5.8 for the same earthquake. On the other hand, the solution produced by the GEOFON system of the Deutsche GeoForschungsZentrum (GFZ), Potsdam, provided M w = 6.3 based on regional stations of the GEOFON system without the use of teleseismic records. The reason is that teleseismic records of the second earthquake have been covered by the records of an earlier large earthquake (M w = 7.4) occurring in the Kermadec-New Zealand region on 4 March 2021 at 17:41:24.2 UTC (Joachim Saul, GFZ; personal communi-cation). We redetermined moment magnitudes from the inversion of teleseismic records for the first earthquake and from the inversion of InSAR data for the three strong earthquakes.

Hypocenter Relocation
To improve the hypocentral determinations of NOA for the three largest earthquakes of the seismic sequence (Table S1) we relocated the hypocenters by implementing a 1D velocity model (Table 1) and by employing the Non-Linear Location (NLLoc) algorithm [19,20], which follows a non-linear location approach [21] and provides more reliable solutions and hypocenter error estimates as compared to linearized algorithms. The velocity model used is a revision of the one found in an earlier study for a nearby area [22]. NLLoc provides a complete, probabilistic solution of the earthquake location problem expressed in terms of the posterior density function (PDF) in the space and time domains. Only manually repicked P and S phases recorded by permanent stations of the Hellenic Unified Seismological Network (https://bbnet.gein.noa.gr/HL/real-time-plotting/husn/husnmap, last access 13 May 2021) situated at epicentral distances of up to~120 km from the first shock were utilized. The reason is that the crustal thinning from the plate boundary towards the backarc area creates significant errors in accurately locating the earthquake, especially when distant seismic phases are included in the analysis [23]. The difference in the onset times picked relative to NOA's solution is on average ± 0.15 s. Following an iteration procedure, the relocation was repeated, by including phase residuals at seismic stations obtained from the previous run, until no further significant decrease of the RMS value between two successive runs was achieved, i.e., until minimizing the location errors. To correct for irregular station distribution a relevant station weighting tool in NLLoc was utilized.

Earthquake Focal Mechanisms
Earthquake focal mechanisms for the three strong earthquakes have been obtained through moment tensor solutions produced by several seismological centers (Table S2). To understand similarities and differences between the various fault-plane solutions we calculated the average strike, dip and rake (slip vector) for each one of the two nodal planes (NPs) involved in the solutions published for each one of the three strong earthquakes. The results received are presented in Section 3.2.

Rupture Process from the Inversion of Teleseismic P-Waveforms
The temporal and spatial evolution of the seismic slip upon the fault plane of the first shock of 3 March 2021 has been modeled applying a non-negative, least squares finite-fault inversion scheme based on a kinematic parameterization of the fault [24][25][26]. The data set used consists of P waveforms from teleseismic records at distances ranging from 30 • to 90 • with good azimuthal coverage ( Figure 4) and a high signal to noise ratio. Waveform data were collected from GEOFON and other seismographic networks and downloaded through the IRIS (Incorporated Research Institutions for Seismology) Data Management Center (https://ds.iris.edu/ds/nodes/dmc/, last access 20 May 2021). The waveforms were cut nearly 2 s before the first P arrival and their duration used for inversion is 25 s. All waveforms were pre-processed to remove the mean offset and instrument response before the inversion. They were also band-pass filtered between 0.04 and 0.6 Hz using a Butterworth filter, re-sampled to 0.2 samples/s and finally integrated in time to obtain displacements. The calculated elementary synthetics were convolved with an attenuation operation under the assumption that t* = 1 s, where t* is the attenuation parameter of teleseismic body waves that represents the total body wave travel time divided by Q along the ray path for P waves. More details on the application of method for earthquakes in the Mediterranean region and elsewhere can be found in several papers e.g., [24,25,27].
downloaded through the IRIS (Incorporated Research Institutions for Seismology) Data Management Center (https://ds.iris.edu/ds/nodes/dmc/, last access 20 May 2021). The waveforms were cut nearly 2 s before the first P arrival and their duration used for inversion is 25 s. All waveforms were pre-processed to remove the mean offset and instrument response before the inversion. They were also band-pass filtered between 0.04 and 0.6 Hz using a Butterworth filter, re-sampled to 0.2 samples/s and finally integrated in time to obtain displacements. The calculated elementary synthetics were convolved with an attenuation operation under the assumption that t* = 1 s, where t* is the attenuation parameter of teleseismic body waves that represents the total body wave travel time divided by Q along the ray path for P waves. More details on the application of method for earthquakes in the Mediterranean region and elsewhere can be found in several papers e.g., [24,25,27]. This method could not be applied to the second strong earthquake of the 4 March since no adequate number of teleseismic P-wave records are available. As explained earlier, the reason is that in many stations the records have been covered by the records of a large earthquake that occurred in the New Zealand-Kermadec seismic zone. On the other hand, the earthquake of 12 March was not large enough to apply the method. Therefore, the investigation of the space and time seismic slip distribution from the inversion of teleseismic P-waveforms was restricted to the first shock of 3 March 2021.
The seismic fault dimensions were chosen large enough to permit the slip upon the fault in all possible directions if it is imposed by the waveform data. Namely, the rectangular fault created was of 32.5 km in length and of 16 km in depth measured from the surface. Taking a fault dip of ~45° the along dip width of the fault is about 25 km. The fault was discretized to 180 sub-faults (cells), 18 of them along strike and 10 along dip. Each point source response was computed with a code based on the generalized ray theory (e.g., [28]) and the use of the crustal velocity model with the specifications shown in Table 1. We adopted the relocated hypocenter obtained after our analysis, which was set at horizontal distance of ~19 km from the southeast edge of the fault. The exact way the synthetics were constructed followed the discussion by [29]. This method could not be applied to the second strong earthquake of the 4 March since no adequate number of teleseismic P-wave records are available. As explained earlier, the reason is that in many stations the records have been covered by the records of a large earthquake that occurred in the New Zealand-Kermadec seismic zone. On the other hand, the earthquake of 12 March was not large enough to apply the method. Therefore, the investigation of the space and time seismic slip distribution from the inversion of teleseismic P-waveforms was restricted to the first shock of 3 March 2021.
The seismic fault dimensions were chosen large enough to permit the slip upon the fault in all possible directions if it is imposed by the waveform data. Namely, the rectangular fault created was of 32.5 km in length and of 16 km in depth measured from the surface. Taking a fault dip of~45 • the along dip width of the fault is about 25 km. The fault was discretized to 180 sub-faults (cells), 18 of them along strike and 10 along dip. Each point source response was computed with a code based on the generalized ray theory (e.g., [28]) and the use of the crustal velocity model with the specifications shown in Table 1. We adopted the relocated hypocenter obtained after our analysis, which was set at horizontal distance of~19 km from the southeast edge of the fault. The exact way the synthetics were constructed followed the discussion by [29].
According to the fault-plane solutions available for the 3 March 2021 strong earthquake (Table S2) the fault plane dips towards either NE or SW. However, a fault dipping NE is consistent with the active tectonics of the area ( Figure 1). After several inversion iterations the strike of 312 • and dip of 45 • were found to better fit the data for the NE-dipping fault. As we will see later, this is consistent with the source solution found after InSAR analysis too. Although fault-plane solutions indicate a rake of about −90 • , which implies nearly pure normal fault, the rake vector was allowed to vary within the range from −60 • to −150 • during the inversion procedure. The procedure was performed repeatedly setting different values of rupture velocity, varying from 2.5 to 3.0 km/s. As the rupture velocity did not change during each one of the inversion iterations, six time windows were inserted with 0.3 s time lag. This option allowed a rise time of up to 1.8 s to be taken on each sub-fault, if required by the observations. Table 2 summarizes the most important parameters involved in the inversion process.

Rupture Process from InSAR
Copernicus is a European Commission initiative, which, in cooperation with the European Space Agency (ESA), developed the new family of Sentinel satellites. Sentinel-1 carries a sophisticated Synthetic Aperture Radar (SAR) instrument for capturing images of the Earth's surface. The radar instrument uses C-band of microwave radiation and can operate in four modes, but the standard product is the Interferometric Wide Swath (IWS). Concerning the IWS polarimetry the wave has a single mode polarization (Vertical-Vertical or Horizontal-Horizontal). The default IWS mode over land has a swath width of 250 km and a ground resolution of 5 × 20 m. This mode images in three sub-swaths using the Terrain Observation with Progressive Scans SAR (TOPSAR), while the actual repeat period is 6 days. One of the main applications of the IWS mode scenes concerns the monitoring of land deformation.
Sentinel-1 SAR scenes assure a series of advantages: (a) continuous, all-weather day and night imagery, (b) rapid revisit period in the same imaging mode (6 days), (c) constant and regular acquisition to build up a large global archive, (d) wide area coverage, thanks to the 250 km image swath width, and (e) narrow orbital tube.
In order to detect and measure surface deformation caused by the earthquake activity during March 2021, Sentinel-1 IWS Single Look Complex (SLC) products (Table S3), covering the study area were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/#/home, accessed on 14 April 2021). The sensing dates range from 25 February to 15 March 2021 for four ascending passes and from 24 February to 14 March 2021 for three descending passes. The digital elevation model (DEM) used is based on the NASA's Shuttle Radar Topography Mission (SRTM) 3 arcseconds DEM (https://srtm.csi.cgiar.org/10, accessed on 14 April 2021), which is of spatial resolution of approximately 90 m/pixel. Spaceborne SAR interferometry produces 3D topographic data of the Earth's surface directly from two SAR images [30]. An extension of the basic technique, called Differential SAR Interferometry (DInSAR), allows measurements of land deformation [31][32][33]. DInSAR exploits the phase difference between two or more coherent complex-valued images, i.e., SLC products, in order to derive path-length differences in the scale of the carrier wavelength and below. The capability of SAR interferometry to remotely monitor areas much wider than traditional surveying techniques, makes this technique particularly suitable for both regional and local scales [34].
In the case of Sentinel-1, DInSAR processing presents some peculiarities because of the special TOPSAR mode used by the SAR sensor to acquire the Interferometric Wide-Swath (IWS) data. The S1 IWS SLC product, used for interferometric applications, consists of three sub-swaths and each sub-swath image consists of a series of bursts. What we needed for our analysis was (i) the minimum of two SAR images taken before (master scene) and after (slave scene) the earthquake event and forming an interferometric pair, and (ii) a DEM of the area under study. Table S4 lists features of the master and slave scenes used in the present study.
The main processing steps followed are the following: co-registration of the two complex images, topography elimination using the DEM of the area, baseline estimation, generation of the interferogram, DEM/ellipsoid interferogram flattening, adaptive filtering/estimation and generation of the interferometric coherence, phase unwrapping, phase to displacement. These steps were performed with the use of the ENVI SARscape ® software (L3Harris Geospatial, Boulder, CO, USA).
After applying the above main steps, DInSAR was only able to measure the path length difference in its Line of Sight (LoS) direction. Therefore, DInSAR-derived displacements represent the 1D deformation along the LoS direction. LoS represents the direction/distance between the SAR sensor and the target at hand. This is an essential and general limitation of SAR systems. By using ascending and descending LoS deformation maps we decomposed LoS displacements [35,36] into vertical (up-down) and horizontal (east-west) deformation maps. After obtaining the shapefile of sampled points from the raster DInSAR map, the coseismic signal was modeled through Non-Linear and Linear inversions [37], thus allowing to infer the geometry, kinematics and slip distribution in the seismic fault.
The displacement model and the determination of a single fault with distributed slip are produced by following several steps. First an image sampling is required, using the LoS displacement, i.e., the observed displacement, to set a number of points to model. The image sampling was carried out using the equally spaced points approach. To produce the surface displacement, which is induced by a rectangular fault plane in a homogeneous and elastic half-space, a well-known elastic dislocation approach [38] was followed through Non-Linear Inversion. The best-fit source is determined by a set of fault parameters (dip, rake, strike, length), which better predict the observed displacement data. After finding the best-fit solution, inversion is carried out to determine the slip distribution over the fault plane constrained via Non-Linear inversion.
The above procedure was applied sequentially for the three strong earthquakes under investigation.

Hypocentral Relocation
The relocated hypocenters of the three strongest earthquakes of the seismic sequence are listed in Table 3 and displayed in Figure 2. As regards the first and second earthquakes the new epicenters are shifted towards SW about 2.5 and 1.5 km away from the respective NOA's epicenters. For the third earthquake, however, the relocated and the NOA's epicenters are very close, which is due to the fact that NOA's routine determinations are improved thanks to the installation and operation of six portable seismograph stations by the Geophysical Laboratory, Aristotelian University of Thessaloniki, before the occurrence of the third earthquake. The RMS found for the three relocated hypocenters are less than 0.22, while in NOA's solutions it is equal to 0.67, 0.47 and 0.41 for the 1st, 2nd and 3rd strong earthquakes, respectively. It is noteworthy that according to the calculated focal depths, h, the second and third earthquakes were found quite shallow, with h = 4.0 km and h = 4.5 km, respectively, as compared to the 3 March first shock (h = 10.5 km). The respective focal depths determined by NOA are 8.5 km, for the first shock, and 4.8 and 7.0 km for the next two earthquakes. However, the centroid depths found from the inversion of teleseisic records (Table 2) and of InSAR data (Section 3.3) as well as from the moment tensor solutions produced by the various centers (Table S1) are in general larger than the focal depths particularly of the second and third earthquakes.

Earthquake Focal Mechanisms
From the results summarized in Table 4, it is revealed that the solutions produced for each one of the three earthquakes are consistent each other since the average strikes, dip and rake have small standard deviations. Comparing the three earthquakes, we found that they are characterized by a strike which is about NW-SE or WNW-ESE. On the other hand, both nodal planes (NPs) of the first two earthquakes have nearly the same average dip, which ranges between 44 • and 49 • but it is either 39 • or 54 • for the third earthquake. The values of the average rake imply that the faulting style of the three earthquakes is nearly pure normal with a small strike-slip component. However, the strike-slip component is relatively larger in the third earthquake. Table 4. Average (A v ) geometric features of the nodal planes (NPs) determined in fault-plane solutions produced by several national and international seismological centers (see Table S2) for the three strong earthquakes.

Date
Time

Seismic Slip Distribution
The synthetics obtained for the 3 March 2021 first shock fit well-enough the recorded waveforms ( Figure 5). The heterogeneous spatial slip distribution, which is illustrated in a SE-NW cross-section of the fault area ( Figure 6), shows that seismic slip occurred in a main patch of the fault plane and in a secondary patch situated to the southeast but at shallower depth with respect to the main one. The main slip area is concentrated at depths from 6 to 15 km with the maximum slip of~1.20 m concentrated at depths of 10-12 km. The minimum slip of~20 cm occurred near the surface. The rupture length, L, changes with depth, being L~15 km at depths from 8 to 15 km and L~22 km from below the surface down to 8 km. The increase of L to the upper part of the fault plane is due to the fact that the secondary slip patch was developed at shallower depths with respect to the main patch. From the inversion procedure we found that the rake values upon the fault plane varies between −80 • and −115 • but the mean rake value in the main patch of slip is about −95 • , which is consistent with the average rake of −93 • (Table 4) calculated from the various fault-plane solutions (Table S2). The total seismic moment released was estimated at Mo = 3.5 × 10 18 N*m, which corresponds to Mw = 6.3. The moment rate function shows that the total duration of the seismic rupture process was ~10.2 s but the main moment release occurred within the first 4 s of rupture (Figure 7). From the time evolution of seismic slip, which is illustrated by six sequential snapshots with time step of nearly 1.7 s (Figure 6), it is evident that after the initiation of the process the rupture propagated mainly upwards for about 7 s. After reaching close to the surface the rupture propagated bilaterally but mainly towards SE where the secondary patch of slip developed.
The horizontal projection of the rupture area of the first shock ( Figure 8) shows that the most affected villages from the first earthquake, like Mesochori and Damasi [1], are situated within the rupture area. On the contrary, in the town of Tyrnavos and in the city of Larissa, which are situated outside the rupture area, less damage was reported. It is also evident that most of the foreshocks that preceded the first shock from the 28 February to the 3 March 2021 occurred within the main asperity that ruptured. The total seismic moment released was estimated at M o = 3.5 × 10 18 N*m, which corresponds to M w = 6.3. The moment rate function shows that the total duration of the seismic rupture process was~10.2 s but the main moment release occurred within the first 4 s of rupture (Figure 7). From the time evolution of seismic slip, which is illustrated by six sequential snapshots with time step of nearly 1.7 s (Figure 6), it is evident that after the initiation of the process the rupture propagated mainly upwards for about 7 s. After reaching close to the surface the rupture propagated bilaterally but mainly towards SE where the secondary patch of slip developed.
The horizontal projection of the rupture area of the first shock ( Figure 8) shows that the most affected villages from the first earthquake, like Mesochori and Damasi [1], are situated within the rupture area. On the contrary, in the town of Tyrnavos and in the city of Larissa, which are situated outside the rupture area, less damage was reported. It is also evident that most of the foreshocks that preceded the first shock from the 28 February to the 3 March 2021 occurred within the main asperity that ruptured.

Ground Deformation from INSAR
The results obtained from InSAR analysis are characterized by coherent values ranging from 0.0 to 1.0, the value 1 corresponding to very high coherence. The DInSAR process produced images in radian unit in the range −π to π, which caused ambiguity problems. Although the pattern of deformation could be detected, the main information regarding the deformation value could not be read properly. In order to get results of deformation containing metric value, an unwrapping process was performed, and the phase unit was transformed into metric units in LoS for every interferometric pair. This unwrapping step was undertaken with the use of the Minimum Cost Flow method [39].

The First Shock of 3 March 2021
The wrapped interferogram of the 3 March earthquake and the corresponding displacement map (Figure 9) based on the interferometric pair are composed by Sentinel 1 SLC images taken in ascending mode on 25 February 2021 and 3 March 2021 for the master and slave images, respectively. A clear pattern of 14 fringes is evident, which forms a lobe showing subsidence. From the displacement map the subsidence amplitude is estimated up to about −40 cm. The subsidence domain is nearly identical with the main patch of rupture obtained from the inversion of P waveforms (Figure 8). In the southern part of the lobe of fringes an incomplete fringe is also recognized corresponding to uplift, the amplitude of which from the displacement map is estimated up to ~3 cm. The comparison of LoS displacements obtained from the observed and modeled data are nearly identical given that only minimal residuals were received ( Figure 10).

Ground Deformation from INSAR
The results obtained from InSAR analysis are characterized by coherent values ranging from 0.0 to 1.0, the value 1 corresponding to very high coherence. The DInSAR process produced images in radian unit in the range −π to π, which caused ambiguity problems. Although the pattern of deformation could be detected, the main information regarding the deformation value could not be read properly. In order to get results of deformation containing metric value, an unwrapping process was performed, and the phase unit was transformed into metric units in LoS for every interferometric pair. This unwrapping step was undertaken with the use of the Minimum Cost Flow method [39].

The First Shock of 3 March 2021
The wrapped interferogram of the 3 March earthquake and the corresponding displacement map (Figure 9) based on the interferometric pair are composed by Sentinel 1 SLC images taken in ascending mode on 25 February 2021 and 3 March 2021 for the master and slave images, respectively. A clear pattern of 14 fringes is evident, which forms a lobe showing subsidence. From the displacement map the subsidence amplitude is estimated up to about −40 cm. The subsidence domain is nearly identical with the main patch of rupture obtained from the inversion of P waveforms (Figure 8). In the southern part of the lobe of fringes an incomplete fringe is also recognized corresponding to uplift, the amplitude of which from the displacement map is estimated up to~3 cm. The comparison of LoS displacements obtained from the observed and modeled data are nearly identical given that only minimal residuals were received ( Figure 10).  The best-fit solution found for the source (Table 5, Figure 11) is in general consistent with the solution determined from the inversion of teleseismic P waveforms. The geodetic source solution provides a seismic fault dipping to NE with strike SE-NW (317°), dip angle of 30° and rake −110°. The geodetic seismic moment (Table 5) calculated through the Okada formalism analyzed in [38] results in earthquake magnitude Mw = 6.3, which coincides with the seismically determined magnitude. The maximum slip found is 1.5 m at depth 5-7 km, which is shallower than the depth of 10-12 km determined from the seismic inversion method. This is probably due to that the geodetic fault dip angle is smaller than the one (45°) found from both the seismic inversion method and the average angle (48°) of the various fault plane solutions published (Table 4).   The best-fit solution found for the source (Table 5, Figure 11) is in general consistent with the solution determined from the inversion of teleseismic P waveforms. The geodetic source solution provides a seismic fault dipping to NE with strike SE-NW (317°), dip angle of 30° and rake −110°. The geodetic seismic moment (Table 5) calculated through the Okada formalism analyzed in [38] results in earthquake magnitude Mw = 6.3, which coincides with the seismically determined magnitude. The maximum slip found is 1.5 m at depth 5-7 km, which is shallower than the depth of 10-12 km determined from the seismic inversion method. This is probably due to that the geodetic fault dip angle is smaller than the one (45°) found from both the seismic inversion method and the average angle (48°) of the various fault plane solutions published (Table 4).  The best-fit solution found for the source (Table 5, Figure 11) is in general consistent with the solution determined from the inversion of teleseismic P waveforms. The geodetic source solution provides a seismic fault dipping to NE with strike SE-NW (317 • ), dip angle of 30 • and rake −110 • . The geodetic seismic moment (Table 5) calculated through the Okada formalism analyzed in [38] results in earthquake magnitude M w = 6.3, which coincides with the seismically determined magnitude. The maximum slip found is 1.5 m at depth 5-7 km, which is shallower than the depth of 10-12 km determined from the seismic inversion method. This is probably due to that the geodetic fault dip angle is smaller than the one (45 • ) found from both the seismic inversion method and the average angle (48 • ) of the various fault plane solutions published (Table 4).

The Earthquake of 4 March 2021
For the second strong seismic event of 4 March, the wrapped interferogram and the displacement map ( Figure 12) are based on the interferometric pair composed by Sentinel 1 SLC images taken in ascending mode on 3 March 2021 and 9 March 2021 for the master and slave images, respectively. The deformation pattern is illustrated by four fringes, which form a lobe showing ground subsidence of estimated amplitude up to about −10 cm. The comparison of LoS displacements obtained from the observed and modeled data are identical given that zero residuals were received ( Figure 13). The earthquake epicenter falls at the north-eastern margin of the subsidence domain, which is consistent with the small focal depth of the earthquake. However, the average centroid depth of ~13 km is projected well inside the subsidence domain.

The Earthquake of 4 March 2021
For the second strong seismic event of 4 March, the wrapped interferogram and the displacement map ( Figure 12) are based on the interferometric pair composed by Sentinel 1 SLC images taken in ascending mode on 3 March 2021 and 9 March 2021 for the master and slave images, respectively. The deformation pattern is illustrated by four fringes, which form a lobe showing ground subsidence of estimated amplitude up to about −10 cm. The comparison of LoS displacements obtained from the observed and modeled data are identical given that zero residuals were received ( Figure 13). The earthquake epicenter falls at the north-eastern margin of the subsidence domain, which is consistent with the small focal depth of the earthquake. However, the average centroid depth of~13 km is projected well inside the subsidence domain.

The Earthquake of 4 March 2021
For the second strong seismic event of 4 March, the wrapped interferogram and the displacement map ( Figure 12) are based on the interferometric pair composed by Sentinel 1 SLC images taken in ascending mode on 3 March 2021 and 9 March 2021 for the master and slave images, respectively. The deformation pattern is illustrated by four fringes, which form a lobe showing ground subsidence of estimated amplitude up to about −10 cm. The comparison of LoS displacements obtained from the observed and modeled data are identical given that zero residuals were received ( Figure 13). The earthquake epicenter falls at the north-eastern margin of the subsidence domain, which is consistent with the small focal depth of the earthquake. However, the average centroid depth of ~13 km is projected well inside the subsidence domain.   The best-fit geodetic source solution (Table 5, Figure 14) provides a seismic fault dipping to about SW with strike ESE-WNW (116°), dip angle of 40° and rake −101°. The geodetic fault solution is consistent with one of the two average NPs found from various seismic fault-plane solutions ( Table 4). The geodetic seismic moment (Table 5) results in earthquake magnitude Mw = 6.2. Seismic slip on the fault occurred at depths ranging from 5 to 20 km but the maximum slip of 1.2 m was found at depth ~7 km.

The Earthquakes of 3 and 4 March 2021 Combined
Displacement decomposition including both the seismic events of 3 and 4 March 2021 was also produced from interferometric pairs in ascending mode, taken on 25 February 2021 and 9 March 2021 for master and slave images, and in descending mode taken on 24 February 2021 and 8 March 2021 for master and slave images. The wrapped interferograms clearly show the migration of the ground deformation towards NW ( Figure  15). The opposite subsidence/uplift pattern in the two earthquakes is evident in the vertical displacement map (Figure 16a), which verifies that the first normal fault dips toward NE while the second dips towards SW. Τhe two earthquakes are also characterized by opposite polarity in the horizontal displacement (Figure 16b).  Table 4). The geodetic seismic moment (Table 5) results in earthquake magnitude M w = 6.2. Seismic slip on the fault occurred at depths ranging from 5 to 20 km but the maximum slip of 1.2 m was found at depth~7 km. The best-fit geodetic source solution (Table 5, Figure 14) provides a seismic fault dipping to about SW with strike ESE-WNW (116°), dip angle of 40° and rake −101°. The geodetic fault solution is consistent with one of the two average NPs found from various seismic fault-plane solutions ( Table 4). The geodetic seismic moment (Table 5) results in earthquake magnitude Mw = 6.2. Seismic slip on the fault occurred at depths ranging from 5 to 20 km but the maximum slip of 1.2 m was found at depth ~7 km.

The Earthquakes of 3 and 4 March 2021 Combined
Displacement decomposition including both the seismic events of 3 and 4 March 2021 was also produced from interferometric pairs in ascending mode, taken on 25 February 2021 and 9 March 2021 for master and slave images, and in descending mode taken on 24 February 2021 and 8 March 2021 for master and slave images. The wrapped interferograms clearly show the migration of the ground deformation towards NW ( Figure  15). The opposite subsidence/uplift pattern in the two earthquakes is evident in the vertical displacement map (Figure 16a), which verifies that the first normal fault dips toward NE while the second dips towards SW. Τhe two earthquakes are also characterized by opposite polarity in the horizontal displacement (Figure 16b).

The Earthquakes of 3 and 4 March 2021 Combined
Displacement decomposition including both the seismic events of 3 and 4 March 2021 was also produced from interferometric pairs in ascending mode, taken on 25 February 2021 and 9 March 2021 for master and slave images, and in descending mode taken on 24 February 2021 and 8 March 2021 for master and slave images. The wrapped interferograms clearly show the migration of the ground deformation towards NW (Figure 15). The opposite subsidence/uplift pattern in the two earthquakes is evident in the vertical displacement map (Figure 16a), which verifies that the first normal fault dips toward NE while the second dips towards SW. The two earthquakes are also characterized by opposite polarity in the horizontal displacement (Figure 16b). The solutions obtained imply that with the second earthquake the activated fault was subparallel and antithetic to the one activated with the first earthquake. This is a reminder of similar seismotectonics characterizing the strong earthquake sequence that ruptured the eastern Gulf of Corinth on 24 and 25 February and again on 4 March 1981 in Central Greece [40,41] to the south of Thessaly (Figure 1). This case is discussed further later.

The Earthquake of 12 March 2021
The ground deformation caused by the third seismic event of 12 March 2021 is illustrated in wrapped interferograms, in LoS displacement maps as well as in vertical and east-west decomposed displacement maps (Figures 17-19). A clear pattern of three fringes forming a relative small lobe showing subsidence is evident (Figure 17). From the displacement maps the ground subsidence is estimated up to about −10 cm, while from the decomposition displacement the subsidence amplitude is estimated up to −9 cm. The comparison of LoS displacements found from the observed and modeled data shows The solutions obtained imply that with the second earthquake the activated fault was subparallel and antithetic to the one activated with the first earthquake. This is a reminder of similar seismotectonics characterizing the strong earthquake sequence that ruptured the eastern Gulf of Corinth on 24 and 25 February and again on 4 March 1981 in Central Greece [40,41] to the south of Thessaly (Figure 1). This case is discussed further later.

The Earthquake of 12 March 2021
The ground deformation caused by the third seismic event of 12 March 2021 is illustrated in wrapped interferograms, in LoS displacement maps as well as in vertical and east-west decomposed displacement maps (Figures 17-19). A clear pattern of three fringes forming a relative small lobe showing subsidence is evident (Figure 17). From the displacement maps the ground subsidence is estimated up to about −10 cm, while from the decomposition displacement the subsidence amplitude is estimated up to −9 cm. The comparison of LoS displacements found from the observed and modeled data shows The solutions obtained imply that with the second earthquake the activated fault was subparallel and antithetic to the one activated with the first earthquake. This is a reminder of similar seismotectonics characterizing the strong earthquake sequence that ruptured the eastern Gulf of Corinth on 24 and 25 February and again on 4 March 1981 in Central Greece [40,41] to the south of Thessaly (Figure 1). This case is discussed further later.

The Earthquake of 12 March 2021
The ground deformation caused by the third seismic event of 12 March 2021 is illustrated in wrapped interferograms, in LoS displacement maps as well as in vertical and east-west decomposed displacement maps (Figures 17-19). A clear pattern of three fringes forming a relative small lobe showing subsidence is evident (Figure 17). From the displacement maps the ground subsidence is estimated up to about −10 cm, while from the decomposition displacement the subsidence amplitude is estimated up to −9 cm. The comparison of LoS displacements found from the observed and modeled data shows small residuals at the southern part of the deformation domain obtained from the ascending mode (Figure 20). In the southern part of the lobe of fringes in the ascending pair ( Figure 17) it is also recognized an incomplete fringe corresponding to ground uplift, which from displacement maps is measured up to 2 cm. From the decomposition displacement ( Figure 18) it is measured up to 3 cm of uplift, but no displacement is observed in the west-east direction ( Figure 19). However, the residuals are equal to zero at the descending mode ( Figure 21). For this reason, this LoS displacement is a preferred one and indicates subsidence of~5 cm in the central side of the deformation field and uplift of~2.5 cm in the NE side of the deformed area.
small residuals at the southern part of the deformation domain obtained from the ascending mode (Figure 20). In the southern part of the lobe of fringes in the ascending pair ( Figure 17) it is also recognized an incomplete fringe corresponding to ground uplift, which from displacement maps is measured up to 2 cm. From the decomposition displacement ( Figure 18) it is measured up to 3 cm of uplift, but no displacement is observed in the west-east direction ( Figure 19). However, the residuals are equal to zero at the descending mode ( Figure 21). For this reason, this LoS displacement is a preferred one and indicates subsidence of ~5 cm in the central side of the deformation field and uplift of ~2.5 cm in the NE side of the deformed area.
As observed also in the case of the 4 March 2021 earthquake, the epicenter of the third earthquake falls at the north-eastern margin of the subsidence domain. This is again consistent with the small relocated (h = 4.5 km) focal depth of the earthquake. However, the average centroid depth of ~9 km falls well inside the subsidence domain.  small residuals at the southern part of the deformation domain obtained from the ascending mode (Figure 20). In the southern part of the lobe of fringes in the ascending pair ( Figure 17) it is also recognized an incomplete fringe corresponding to ground uplift, which from displacement maps is measured up to 2 cm. From the decomposition displacement ( Figure 18) it is measured up to 3 cm of uplift, but no displacement is observed in the west-east direction ( Figure 19). However, the residuals are equal to zero at the descending mode ( Figure 21). For this reason, this LoS displacement is a preferred one and indicates subsidence of ~5 cm in the central side of the deformation field and uplift of ~2.5 cm in the NE side of the deformed area. As observed also in the case of the 4 March 2021 earthquake, the epicenter of the third earthquake falls at the north-eastern margin of the subsidence domain. This is again consistent with the small relocated (h = 4.5 km) focal depth of the earthquake. However, the average centroid depth of ~9 km falls well inside the subsidence domain.     The best-fit geodetic source solution for the earthquake of 12 March 2021 (Table 5, Figure 22) provides a seismic fault dipping to about SW with strike ESE-WNW (117°), dip angle of 45° and rake −95°. This solution is similar to the one obtained for the second earthquake of 4 March 2021 as well as to one of the two average NPs found from various seismic fault-plane solutions ( Table 4). The geodetic moment (Table 5) results in magni-   The best-fit geodetic source solution for the earthquake of 12 March 2021 (Table 5, Figure 22) provides a seismic fault dipping to about SW with strike ESE-WNW (117°), dip angle of 45° and rake −95°. This solution is similar to the one obtained for the second earthquake of 4 March 2021 as well as to one of the two average NPs found from various seismic fault-plane solutions ( Table 4). The geodetic moment (Table 5) results in magni-   The best-fit geodetic source solution for the earthquake of 12 March 2021 (Table 5, Figure 22) provides a seismic fault dipping to about SW with strike ESE-WNW (117°), dip angle of 45° and rake −95°. This solution is similar to the one obtained for the second earthquake of 4 March 2021 as well as to one of the two average NPs found from various seismic fault-plane solutions ( Table 4). The geodetic moment (Table 5) results in magni- As observed also in the case of the 4 March 2021 earthquake, the epicenter of the third earthquake falls at the north-eastern margin of the subsidence domain. This is again consistent with the small relocated (h = 4.5 km) focal depth of the earthquake. However, the average centroid depth of~9 km falls well inside the subsidence domain.
The best-fit geodetic source solution for the earthquake of 12 March 2021 (Table 5, Figure 22) provides a seismic fault dipping to about SW with strike ESE-WNW (117 • ), dip angle of 45 • and rake −95 • . This solution is similar to the one obtained for the second earthquake of 4 March 2021 as well as to one of the two average NPs found from various seismic fault-plane solutions ( Table 4). The geodetic moment (Table 5) results in magnitude M w = 5.7, which is consistent with the M w~5 .6 determined from seismic moment tensor solutions (Table S1). Seismic slip on the fault occurred mainly in two distinct patches at depth of~5 km the first and~10 km the second (Figure 22). The maximum slip of~1.0 m, however, occurred within the deeper patch, which is situated to the southwest with respect to the first.
, 311 19 of 24 tude Mw = 5.7, which is consistent with the Mw~5.6 determined from seismic moment tensor solutions (Table S1). Seismic slip on the fault occurred mainly in two distinct patches at depth of ~5 km the first and ~10 km the second ( Figure 22). The maximum slip of ~1.0 m, however, occurred within the deeper patch, which is situated to the southwest with respect to the first.

Discussion
The seismotectonic implications of the results obtained from the source solutions for the first and second earthquakes are important since they imply that the fault activated with the second earthquake is subparallel and antithetic to the one associated with the first earthquake (Figures 23 and 24). This is similar to the seismotectonics characterizing the strong earthquake sequence that ruptured the area of the eastern Gulf of Corinth on 24 February 1981 (Mw = 6.58), 25 February 1981 (Mw = 6.32) and again on 4 March 1981 (Mw = 6.23) (magnitudes taken from [40]) in Central Greece to the south of Thessaly province (Figure 1). During the first two 1981 earthquakes, surface fault-ruptures striking roughly WSW-ENE and dipping to ~NNW were detected, while a roughly parallel but antithetic fault was activated with the third earthquake, e.g., [41,42]. However, during March 2021 the faults activated are blind, unmapped and remained unknown so far. Complex normal fault systems that involve antithetic fault segments have been described in several cases of seismic sequences with multiple earthquake ruptures, e.g., the Kozani-Grevena (NW Greece) earthquake (Mw = 6.5) of 13 May 1995 [43], the 26 December 2003 Ban earthquake (Mw = 6.5) in Iran [44], and the Central Italy earthquakes on 24 August and 26 and 30 October 2016 (Mw = 6.1, 5.9, and 6.5) [45].
The four times less amplitude of the geodetic subsidence displacement, d, found for the second strong earthquake of 4 March (d = −10 cm), with respect to the one found for the first shock of 3 March (d = −40 cm), is better explained by a moment magnitude of 6.0 or 6.1, produced by national seismological institutes for the second earthquake (Table S1), or even by our geodetic estimation of Mw = 6.2, than by the Mw = 6.3 determined by

Discussion
The seismotectonic implications of the results obtained from the source solutions for the first and second earthquakes are important since they imply that the fault activated with the second earthquake is subparallel and antithetic to the one associated with the first earthquake (Figures 23 and 24). This is similar to the seismotectonics characterizing the strong earthquake sequence that ruptured the area of the eastern Gulf of Corinth on 24 February 1981 (M w = 6.58), 25 February 1981 (M w = 6.32) and again on 4 March 1981 (M w = 6.23) (magnitudes taken from [40]) in Central Greece to the south of Thessaly province (Figure 1). During the first two 1981 earthquakes, surface fault-ruptures striking roughly WSW-ENE and dipping to~NNW were detected, while a roughly parallel but antithetic fault was activated with the third earthquake, e.g., [41,42]. However, during March 2021 the faults activated are blind, unmapped and remained unknown so far. Complex normal fault systems that involve antithetic fault segments have been described in several cases of seismic sequences with multiple earthquake ruptures, e.g., the Kozani-Grevena (NW Greece) earthquake (M w = 6.5) of 13 May 1995 [43], the 26 December 2003 Ban earthquake (M w = 6.5) in Iran [44], and the Central Italy earthquakes on 24 August and 26 and 30 October 2016 (M w = 6.1, 5.9, and 6.5) [45]. geodetically (Table 5) and from moment tensor solutions (Table 4) is 42°. A Mw = 5.7 was determined, which is consistent with the Mw = 5.6 found from moment tensor solutions, while the maximum slip of 1.0 m occurred at depth of ~10 km. However, further research is needed in the future, e.g., regarding space-time aftershocks distribution, to verify that our preferred fault geometry is a valid one. Otherwise, we do not rule out that a fault geometry similar to that associated with the first shock may fit better to the case of the 12 March aftershock.  (Tables 4 and 5): 314°/41° for SF1, 123°/44° for SF2, 112°/42° for SF3. Positions of the three faults were adopted by taking into account the fault models, resulted from the inversion of teleseismic P waveforms and of InSAR data, as well as the respective relocated focal depth (Table 3) and the average dip angle for each fault. Lengths of about 23, 17 and 11 km were calculated for SF1, SF2 and SF3, respectively, by taking into account rupture lengths determined from our seismic and geodetic fault models (Figures 8, 11, 14 and 22) and a worldwide empirical relationship between subsurface rupture length and magnitude found for normal earthquakes [46]. Εpicenters (stars) and magnitudes of the strong earthquakes of 3, 4 and 12 March 2021, as well as of the main faults (black lines) in the Thessaly area, as in Figure 1.
The relocated epicenters of both the second and third earthquakes fall at the NE margin of the down-dip projection of the geodetic fault solutions. One may argue that epicentral locations shifted by a few km to SW would fit better the SW-dipping fault planes. However, such an epicentral shift requires focal depths of ~14.5 and ~11.5 km instead of the very shallow relocated hypocenters of 4.0 and 4.5 km found for the 2nd and 3rd earthquakes, respectively. However, the surface projection of the average centroid depth of ~9-10 km of these earthquakes is shifted towards SW and falls well inside the subsidence domains.
The 3 March earthquake was preceded by a short-living foreshock sequence occurring from 28 February to 3 March 2021. Most of the foreshocks occurred within the main  (Tables 4 and 5): 314 • /41 • for SF1, 123 • /44 • for SF2, 112 • /42 • for SF3. Positions of the three faults were adopted by taking into account the fault models, resulted from the inversion of teleseismic P waveforms and of InSAR data, as well as the respective relocated focal depth (Table 3) and the average dip angle for each fault. Lengths of about 23, 17 and 11 km were calculated for SF1, SF2 and SF3, respectively, by taking into account rupture lengths determined from our seismic and geodetic fault models (Figures 8, 11, 14 and 22) and a worldwide empirical relationship between subsurface rupture length and magnitude found for normal earthquakes [46]. Epicenters (stars) and magnitudes of the strong earthquakes of 3, 4 and 12 March 2021, as well as of the main faults (black lines) in the Thessaly area, as in Figure 1.  Our seismotectonic interpretation for the March 2021 earthquake sequence is different from the one that considers that the three earthquakes ruptured at spatially sequential segments of a normal fault system striking SE-NW and dipping NE [16]. Such a fault model requires focal depths of ~15 km to explain the position of the epicenters of the second and third earthquakes at the northeast sides of the respective surface deformation fields, a point already discussed as regards our seismotectonic model. In addition, the NE/SW uplift/subsidence pattern recognized in our geodetic results for the second and third earthquakes (Figures 13, 16 and 21) is better explained by the fault model proposed in this paper. However, our interpretation is consistent with the suggestion [16] that the The four times less amplitude of the geodetic subsidence displacement, d, found for the second strong earthquake of 4 March (d = −10 cm), with respect to the one found for the first shock of 3 March (d = −40 cm), is better explained by a moment magnitude of 6.0 or 6.1, produced by national seismological institutes for the second earthquake (Table S1), or even by our geodetic estimation of M w = 6.2, than by the M w = 6.3 determined by GEOFON/GFZ, which perhaps is an overestimation.
The third earthquake of 12 March ruptured to the NW continuation of the second rupture. The comparison of LoS displacements found from the observed and modeled InSAR data shows small residuals at the southern part of the deformation domain obtained from the ascending mode ( Figure 20). However, the residuals are equal to zero at the descending mode ( Figure 21). For this reason, the last LoS displacement is a tentatively preferred one and indicates subsidence of~5 cm in the central side of the deformation field and uplift of~2.5 cm in the NE side. Our best-fit geodetic source solution for the third earthquake (strike/dip/rake: 117 • /45 • /−95 • ) is quite similar to the one obtained for the second earthquake and consistent with published moment tensor solutions (average values for one nodal plane: 106 • /39 • /−105 • ). The mean of the dip angle determined geodetically (Table 5) and from moment tensor solutions (Table 4) is 42 • . A M w = 5.7 was determined, which is consistent with the M w = 5.6 found from moment tensor solutions, while the maximum slip of 1.0 m occurred at depth of~10 km. However, further research is needed in the future, e.g., regarding space-time aftershocks distribution, to verify that our preferred fault geometry is a valid one. Otherwise, we do not rule out that a fault geometry similar to that associated with the first shock may fit better to the case of the 12 March aftershock.
The relocated epicenters of both the second and third earthquakes fall at the NE margin of the down-dip projection of the geodetic fault solutions. One may argue that epicentral locations shifted by a few km to SW would fit better the SW-dipping fault planes. However, such an epicentral shift requires focal depths of~14.5 and~11.5 km instead of the very shallow relocated hypocenters of 4.0 and 4.5 km found for the 2nd and 3rd earthquakes, respectively. However, the surface projection of the average centroid depth of~9-10 km of these earthquakes is shifted towards SW and falls well inside the subsidence domains.
The 3 March earthquake was preceded by a short-living foreshock sequence occurring from 28 February to 3 March 2021. Most of the foreshocks occurred within the main asperity that ruptured. A different pattern was found regarding the~6-month lasting foreshocks preceding the 25 October 2018 mainshock (M w = 6.8), associated with the thrust-obliqueslip rupture offshore Zakynthos Isl., Ionian Sea [47]. Namely, the foreshocks bounded the mainshock asperity up-dip and at the north of it. However, the two largest imminent foreshocks (M w = 4.1, M w = 4.8) occurred very close to the 2018 mainshock hypocenter as happened also with the imminent foreshocks of the 3 March 2021 shock.
Our seismotectonic interpretation for the March 2021 earthquake sequence is different from the one that considers that the three earthquakes ruptured at spatially sequential segments of a normal fault system striking SE-NW and dipping NE [16]. Such a fault model requires focal depths of~15 km to explain the position of the epicenters of the second and third earthquakes at the northeast sides of the respective surface deformation fields, a point already discussed as regards our seismotectonic model. In addition, the NE/SW uplift/subsidence pattern recognized in our geodetic results for the second and third earthquakes (Figures 13, 16 and 21) is better explained by the fault model proposed in this paper. However, our interpretation is consistent with the suggestion [16] that the activated fault system is clearly belonging to the Tyrnavos Graben that started forming in the middle-late Pleistocene, and whose bordering structures are still in a growing phase.

Conclusions
From our results we concluded that the sequence of strong earthquakes that occurred in Thessaly province, Central Greece, during March 2021 was produced by a quite complicated rupture process. The three strong earthquakes of the 3, 4 and 12 March 2021, were associated with three main blind, unmapped and unknown so far fault segments of nearly pure normal faulting with a small strike-slip component. The entire system of normal fault segments activated forms a graben-like structure. The first earthquake of 4 March (M w = 6.3) was produced by a fault segment striking SE-NW and dipping towards NE.
With the strong earthquake of 4 March (M w = 6.2) the rupture propagated further NW but in an antithetic fault segment dipping~SW and striking~SE-NW. A segment of possibly the same fault, striking ESE-WNW was ruptured with the strong earthquake (M w = 5.7) of 12 March. The new findings shed light in the Thessaly seismotectonics, particularly in the north side of this area, which remained seismically silent since AD 1781. The new knowledge acquired provides important insights for the better assessment of the seismic hazard level in the area.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/geosciences11080311/s1, Table S1: Source parameters of the three strongest earthquakes recorded in Thessaly area during March 2021; Table S2: Geometric features of the two nodal planes determined in moment tensor solutions produced by several seismological centers for the three strong earthquakes; Table S3: Specifications of the Sentinel-1 SAR SLC images used in our InSAR analysis; Table S4: Master and slave scenes forming interferometric pairs used to detect and measure deformation of the three earthquake events.