Geodetic Model of the March 2021 Thessaly Seismic Sequence Inferred from Seismological and InSAR Data

: In this work, we propose a geodetic model for the March 2021 Thessaly seismic sequence (TSS). We used the interferometric synthetic aperture radar (InSAR) technique and exploited a dataset of Sentinel-1 images to successfully detect the surface deformation caused by three major events of the sequence and constrain their kinematics, further strengthened by seismic data analysis. Our geodetic inversions are consistent with the activation of distinct blind faults previously unknown in this region: three belonging to the NE-dipping normal fault associated with the Mw 6.3 and Mw 6.0 events, and one belonging to the SW-dipping normal fault associated with the Mw 5.6, the last TSS major event. We performed a Coulomb stress transfer analysis and a 1D pore pressure diffusivity modeling to investigate the space–time evolution of the sequence; our results indicate that the seismic sequence developed in a sort of domino effect. The combination of the lack of historical records of large earthquakes in this area and the absence of mapped surface features produced by past faulting make seismic hazard estimation difﬁcult for this area: InSAR data analysis and modeling have proven to be an extremely useful tool in helping to constrain the rupture characteristics.


Introduction
The Thessaly seismic sequence (TSS), Central Greece, started on 3 March 2021 (10:16:08 UTC) with a Mw 6.3 event and struck an area about 25 km WNW of the town of Larissa ( Figure 1). Fortunately, it caused no fatalities, while only three people were slightly injured due to the partial collapse of buildings with load-bearing masonry walls in the village of Damasi; moreover, widespread liquefaction phenomena were observed in the fields located close to local rivers [1]. Based on various seismological observatories and institutes (AUTH, NOA, GCMT, USGS, OCA, INGV, GFZ), the first event was characterized by a normal fault mechanism (Figure 1a). Several aftershocks (Mw > 5) concentrated to the southeast during the first few hours after the event within a distance of 10 km (Figure 1a).
In the following days, TSS was characterized by other two major events: a Mw 6.0 on March 4 (18:38:19 UTC), which has been localized about 7 km northwest from the first event; and a Mw 5.6 on March 12 (12:57:50 UTC), localized 12 km further toward the NW. A large number of smaller events have also been recorded until mid-April when the sequence decreased in frequency and magnitude. Historical earthquakes that occurred in Thessaly area are also shown (data by [3]). Recent earthquakes occurred in an area not covered by past events, thus filling a seismic gap.
The Thessaly region is located at the back-arc area of the Aegean microplate and is one of the most seismically active areas of Greece [2,4]. Since 1 January 2008, more than 320 revised events have been detected from the Institute of Geodynamics-National Observatory of Athens database (IG-NOA) in the extended fault area near Tyrnavos: this background seismicity at the end of 2020 had a magnitude value ranging between 0.7 to 3.6. The seismic activity resumed on 9 January 2021, with events with a magnitude ranging between 1.5 to 2. A comprehensive description of the historical and instrumental seismicity of the area is provided in detail in [3,5].
This region is mainly characterized by NW-SE-striking normal faults, testifying an extensional stress regime with low strain rate. Several active fault zones have been detected into the central and northern part of Greece and some of them have been associated with past destructive earthquakes [6] (Figure 1b). In particular, the seismicity of the Thessaly Basin originates mainly along two fault zones: the northern one, concentrated in the Larissa plain, and a southern fault zone [7]. These two zones suggest that there has been a seismic gap in the Larissa plain [8]. Overall, the 2021 TSS sequence affected the area surrounded by past events, thus filling the gap between previous ruptures.
Concerning the structures involved in TSS, a sporadic series of NW-SE striking surface breaks were found on the mountains to the north of Zarko village [1]. Unfortunately, no tectonic ruptures were found corresponding to the Mw 6.0 and Mw 5.6 events. Therefore, the geodetic data play a crucial role in locating the faults and in better constraining their geometry.
Here, we exploit the complete dataset of Sentinel-1 (S1) InSAR measurements from both ascending and descending orbits to investigate the ground displacement field. The S1 short repeat time made it possible to isolate the three major events that characterized the sequence (Figure 2a). In addition, given the rather low magnitude of the aftershocks (Mw ≤ 5.0) compared to the main events, we assumed that they only marginally contributed to the overall surface deformation. We modeled the available InSAR deformation maps to retrieve the parameters characterizing some finite dislocation sources. Then we performed a Coulomb stress transfer and a 1D diffusivity analysis in order to investigate possible fault interactions. Finally, we highlight the key role played by the configuration of the Thessaly Basin characterized by blind faults interconnected at depth, particularly interesting from both the neotectonic and active tectonics points of view.

Seismological Data
In order to constrain the geometry and location of the main fault structures involved during the TSS, we considered 1853 earthquakes that occurred in the area from 28 February 2021 to 26 April 2021 with magnitudes ranging between 0.2 and 6.3 ( Figure 2). The analyzed data were collected by the Hellenic Seismic Network operated by the Institute of Geodynamics of the National Observatory of Athens (NOA-IG). For the absolute earthquake location, we used the phase readings provided by NOA-IG at 116 seismic stations ( Figure S2). Our database is composed of 23,371 P-and 21,455 S-arrival times. Phase arrival times were re-weighted according to the epicentral distance provided by a preliminary location. As a result, most of the data (about 69% of P-waves and about 71% S-wave) used for absolute earthquake locations were recorded by stations located within an epicentral distance of 100 km.
Absolute earthquake locations were calculated according to a probabilistic approach using the NonLinLoc code [9] and a local 1-D velocity model [10] (see Table 1). The VP/VS ratio used for the analysis was 1.78: this choice corresponds to the best value that allows to minimize the rms and location errors obtained after re-locating the data subset of the main earthquakes (M W ≥ 5.0). We obtained the mean error values, vertical, and horizontal equal to 0.67 km, 0.81 km, and a mean rms of 0.26 s, respectively. Then, we applied a double difference earthquake location technique (HypoDD [11]) to compute the earthquake relative locations that minimize the residuals between the observed and theoretical travel time differences (or double differences) for pairs of earthquakes at each station. We obtained the relative locations of 1760 earthquakes with mean errors, vertical, and horizontal, equal to 242 m and 246 m, respectively. The mean rms was 0.08 s ( Figure S3). To achieve proper least squares error estimates, we followed the method proposed in [11]. The double difference earthquake locations showed a clustered seismicity, tighter along the NE-SW direction with clear earthquake alignments in vertical sections, which was not evident in the absolute catalog locations. For more details about the use of double difference technique (see the Supplementary Materials).
Earthquake locations (Figure 1 and Figure S4) suggest a 40 km NW-SE alignment, which is consistent with the regional trend of the mapped surface faults, although the earthquake clusters of 5-10 km length revealed secondary fault structures. This result is supported by the vertical section plots, which clearly show NE-dipping planes, in agreement with the retrieved focal mechanism solutions (Figure 2b). We pointed out that the continuity of these alignments seems to be interrupted laterally. As shown in Sections 2 and 4, earthquake locations seem to show NE-and SW-dipping planes, replicating an antithetic-synthetic fault set, typical of areas involved in normal faulting. In the northern area, the earthquake location highlights a S-dipping fault plane (Section 5). Our results show a segmented fault zone with predominantly NE-trending structures and complex geometries, typical of normal fault systems as those activated by Mw 6.1 L'Aquila and Mw 6.5 Central Italy earthquakes [12,13].

InSAR Measurements
The TSS was imaged by Sentinel-1 (S1) with a regular six day repeat interval along two ascending (tracks 102 and 175) and two descending (tracks 7 and 80) orbits; we considered the interval from 24 February 2021 to 15 March 2021, which spanned the whole seismic sequence and included four acquisitions per track. The complete list of Sentinel-1 images used in this study is reported in Table S1 and are shown graphically in Figure S1, Supplementary Materials. Sentinel-1 data were processed by using our own internally developed InSAR processing chain [14][15][16].
Notably, almost simultaneous acquisitions from ascending and descending orbits were available, which may support the joint inversion of these two components of the displacement vector. A peculiarity of this S1 dataset is that the second image of track 102 was acquired just in between the first two largest earthquakes (Mw 6.3 and Mw 6.0); in this way, although they occurred only 32 h apart, we could separate the two contributions. This was not the case for the other three orbits, of which the interferograms between the first and second acquisition contained no signal and, therefore, were useless for the following analysis.
To ease the discussion, we defined several "timeframes" according to the Sentinel-1 passes and to the events they include (see colored bars in Figures 2a and S1). In particular, we defined T1 and T2 as the two intervals covering the sole Mw 6.3 and Mw 6.0, respectively ( Figures S5 and S6). In addition, each interferogram was uniquely labelled with a letter according to Table S1 and Figure S1.
Furthermore, we defined T3 as the timeframe interval including both the first and second main events; in this case, we had seven interferograms (blue bars in Figures 2a and S1): four "long" (from end of February to 8-9 March) and three "short" interferograms (from 2-3 to 8-9 March); clearly, the short ones did not include track 102 since it covered one event only. Figure 2 shows the four T3 "long" interferograms (wrapped and unwrapped) coupled by the same acquisition times. In particular, track 102 (ascending, IW1) and track 07 (descending, IW3) are shown in the first row, while track 175 (ascending, IW3) and track 80 (descending, IW1) in the second one. These maps show three NNW-SSE striking deformation lobes located in the Larissa Basin; to the south, the two largest lobes, associated with the Mw 6.3, were evident, particularly on the descending data; this displacement pattern reached a maximum value of~−38 cm in LOS (line-of-sight) on ascending data (negative LOS values represent increasing distances from the satellite). To the NW of this pattern, a smaller lobe relevant to the Mw 6.0 was observed. We also note that the T3 ascending pattern was quite consistent with the superimposition of T1 and T2 (Figures S5 and S6) once the different look angles were taken into account. The corresponding descending interferograms presented similar characteristics, considering the effect of opposite side views and possible horizontal components. Furthermore, we defined T4 as the interval covering the Mw 5.6 only, which was imaged once by all four tracks (yellow bars in Figures 2a and S1). The corresponding wrapped and unwrapped interferograms are presented in Figure 2, showing a prevalent negative LOS displacement pattern extending along the EW direction, with a minimum of~−10 cm, consistent with the lower energy of this event with respect to the previous ones.
Finally, we introduce T5 as the overall timeframe that includes all three main events; similarly to T3, this was formed by seven interferograms (magenta bars in Figures 2a and S1), which are presented in Figures S5 and S6. Note that multiple interferograms within a single timeframe (Figures 3, S5 and S6) basically comprise the same displacement signal since the four considered tracks were shifted in time by no more than one day and there was no significant activity between the end of February and the beginning of the sequence. As a further element for the interpretation, we also estimated the vertical and eastwest (EW) components of the displacement [17] for each of the considered timeframes; for brevity, they are reported in Figure S7.

Analytical Modeling
We jointly modeled the LOS displacements retrieved from both ascending and descending InSAR data with a finite dislocation fault in an elastic and homogeneous halfspace [18,19]. Our source modeling approach was based on a consolidated two-step approach, also implemented in SARscape ® software, capable of computing the distribution of the (non-uniform) slip over the fault plan. The first step was implemented by fitting (via a non-linear inversion) a uniform slip model to the data to constrain the geometry of the faulting plane; in particular, it found the east, north, and depth position of the plane center, along with its strike (azimuth) and dip angles.
Clearly, this first step also provides us with an average slip value (and its rake) and with the fault plane dimensions (length and width); however, these parameters are refined afterward. In fact, the second step computes (via a linear inversion) the distribution of the slip values over the defined fault plane by partitioning it into a number of small patches; at this point, the distribution of the patches showing a significant nonzero slip defines the extension of the actual faulting portion of the plane.
We started with a non-linear inversion scheme [20,21] to identify the fault parameters and mechanism (strike, dip, rake, slip, fault location, length, depth, width) using some of the results derived from the seismic waveform analysis (e.g., strike direction or fault depth) for some inversion cases (see Figures S8, S9, S11 and S12). The optimization starts from a random fault configuration within the given parameter ranges and keeps minimizing the cost function Φ, based on the weighted squares of the residuals between the observed and the predicted data: σ i where d i,obs and d i,mod are the observed and modeled displacements of the i-th data point, respectively; σ i is the weighting parameter (here we assume σ i to be the same for all points); and N is the total number of points used for the inversion. The downhill algorithm is implemented with multiple restarts to guarantee the convergence to the global cost function absolute minimum.
In order to refer the source depth to the actual ground surface, we also applied a topography compensation [22] and assessed possible residual offsets and ramps in the InSAR measurements caused, for instance, by residual orbital signals [23]. Moreover, to reduce the computation load of the inversion process, the InSAR data were preliminarily down-sampled to a 140 m regular grid.
After defining the fault geometry, we applied a linear inversion to obtain the slip distribution, extending the fault length and width to let the slip vanish to zero and subdividing the fault plane into sub-faults of 0.5 × 0.5 km 2 . To avoid back-slip, we used a trial-and-error approach to define the system damping, which is the empirical parameter balancing the data fit and the slip distribution roughness. The adopted linear inversion scheme is described by the equation [24,25]: where d is the InSAR data vector; G is the Green's functions matrix; and ∇ 2 is the Laplacian operator, which is tuned with the damping factor ε obtained by trial-and-error [26] to obtain a reliable slip distribution, described with the m vector of parameters. A further constrain of parameter positivity was adopted to prevent back-slip. The two timeframes that minimize the number of included events while allowing us to study the TSS deformation jointly from both ascending and descending orbits are clearly T3 and T4. We inverted for the two events (Mw 6.3 and Mw 6.0) included in T3 and, separately, for the one (Mw 5.6) included in T4 (Figures S8-S10).
The best-fitting fault model with uniform slip for T3 consists of three distinct~58 • , 33 • , and~42 • NE-dipping fault segments located beneath the three observed lobes of deformation shown in Figure 4 (say Figures 1a,b and 2, respectively). They both extend to a maximum depth of 10 km, the F1a showing a small, left lateral strike-slip component (rake: −118 • ). Although we cannot exclude some effects on the first event, the F2 fault is interpreted as associated with the second event, as confirmed also by modeling the second event (T2) alone (see Figure S11); this showed a 13 × 3.7 km 2 source plain. Conversely, data in T4 were best-fitted by a uniform slip model with a single~47 • SE-dipping fault (say F3). We performed uniform slip modeling for all the timeframes defined earlier; for sake of completeness, we report all these results in the Supplementary Materials in Figures S8-S12 (including inversion parameters and statistics). Our final result was achieved by constraining a small number of the model parameters either to a preliminary geodetic model or to seismic data, as explained in the following.
First, to explore the sources responsible for Mw 6.3, we chose to invert the InSAR data for the shortest allowed timeframe (i.e., T1). We noted that a single rectangular fault cannot reproduce the observed bi-lobed fringe pattern associated with this event; as a test, we performed a free inversion of T1 and found that a good fit with the data was reached only when using a two-fault model ( Figure S11).
Second, we note that, while the LOS displacement profiles for the Mw 6.3 and Mw 5.6 events were very abrupt, thus clearly indicating the corresponding dipping directions, this was not the case for the Mw 6.0 event, as shown in Figure 5. It is clear, therefore, that the use of the displacement information alone makes it difficult to resolve for the F2 strike value. To overcome this problem, we thus inverted T2, constraining its depth and limiting the strike search interval according to the seismic analysis ( Figure S11). Among the dislocation model parameters, the seismological analysis can provide strike and dip values, along with information about the depth. Here, we only used seismic-derived strike and depth information to guide the inversion of T2, while all the other inversions were kept free; therefore, seismic analysis data were used to independently validate the location and mechanism of the retrieved faults other than F2 ( Figure 2b).
As an additional analysis, we also inverted for T5, which contained the whole sequence, by applying the same constraints used for T3 and T4; basically, we obtained the same source parameters when using a different combination of all the available interferograms covering the period ( Figure S13). Subsequently, we selected the geometric settings derived from our preferred non-linear inversions ( Figure S8) for the slip distribution calculation. To this end, we subdivided the fault planes into small patches (0.5 × 0.5 km 2 ) and solved for the slip values on each patch ( Figure 6). For the Mw 6.3, the slip was concentrated between 2 km and 4 km depth: the maximum slip for F1a and F1b faults was about 1.1 m at 2 km and 4 km depths, respectively. Both fault plains present small portions of slip at shallow depths, in the area where ruptures are observed on the ground (Figure 1a). The Mw 6.0 event ruptured a western~15 km-long segment of the F2 fault system with a maximum slip value reaching~1.0 m at 5 km depth. Finally, the M W 5.6 ruptured the~10 km-long segment that had remained unbroken after the previous events. We found a main patch of slip (with maximum of 0.57 m) located at the center of the fault plane at a 2.5 km depth. The overall geodetic moment of the F1 fault system was 2.88·10 18 Nm corresponding to a Mw 6.3 earthquake, while for F2 and F3, we obtained 1.1·10 18 Nm and 3.1·10 17 Nm, corresponding to Mw 6.0 and 5.6 events, respectively. Looking at Figure 6, most of the earthquakes occurred outside or at the borders of the slipped areas on each of the fault planes. This feature suggests that some of the shear stress was redistributed outside the main patches of the slip after each event.

Coulomb Stress Transfer Analysis
To study the interaction between the three main events of the sequence, we calculated the change in the stress field on each fault plane in terms of Coulomb failure function (∆CFF) by using an approach similar to that proposed in [27]. We calculated the co-and postseismic deformation and the associated Coulomb stress change ∆CFF = ∆τ + µ(∆σ + ∆p), where ∆τ and ∆σ are the shear and normal stress change, respectively; µ is the friction coefficient; and ∆p is the pore pressure change. We used a computer code based on the viscoelastic-gravitational dislocation theory [28] and assumed that ∆p = 0, corresponding to drained conditions, and used µ = 0.4. The method allowed for the use of finite source fault models with heterogeneous slip distribution. It uses the standard linear solid rheology defined by three parameters: the unrelaxed shear modulus µ 0 , the viscosity η, and the parameter α, which is the ratio of the fully relaxed modulus to the unrelaxed modulus. As a point of difference with usual analyses that consider the location of the nucleation of the following earthquakes relative to the induced stress variation, here, we investigated the heterogeneity of the stress field on the whole fault plane and the time evolution of the earthquake preparatory process.
As undertaken for the earthquake location, we also used here the layered model proposed in [10] (see Table 1), while the detailed slip distributions and the fault geometry were those inferred from the InSAR modeling. We computed Green functions for a 10-day time window to include the effect of three main events contributing to the stress field. We used 60 equally spaced horizontal points at a distance range of 0-60 km, and 40 points in depth, ranging between 0 and 20 km. Stress field variation was computed on 16 different layers with depths ranging between 0 and 15 km.
We first considered the effect of the Mw 6.3 event on the fault planes of Mw 6.0 and Mw 5.6. Next, we considered the additional effect of the Mw 6.0 event on the fault plane of the Mw 5.6 earthquake. The first interaction is shown in Figure 7a,b while the second is shown in Figure 7c. The results reported in Figure 7a indicate that the Mw 6.3 earthquake strongly affected the Mw 6.0 by producing a significant ∆CFF increase both east and west of the hypocenter and a slighter increase in the upper right corner that likely promoted its next breaking. Similarly, the Mw 6.3 also increased the ∆CFF in the lower part of the fault plane of the Mw 5.6. The occurrence of the Mw 6.0 increased the ∆CFF in the central part of the fault, leading to the following breaking; notably, the earthquakes that occurred before the Mw 6.0 (reported as square in Figure 7a) seemed to encircle the patches that will slip. Moreover, all the aftershocks (reported as crosses in Figure 7a) occurred in the area that was initially inhibited by the ∆CFF and where remaining stress was released after the main event. Although less evident, a similar picture is depicted in Figure 7c for the Mw 5.6.

Diffusivity Analysis
The results reported in the previous section suggest a clear static stress contribution to the evolution of the seismic sequence. However, in order to further investigate the delayed triggering of the main events, we tested the hypothesis of pore pressure diffusion from the region enclosing the TSS. To this aim, we modeled the pore pressure perturbation caused by a point source in an isotropic fluid saturated medium. We considered the Mw 6.3 as the source and retrieved the best 1D isotropic diffusivity value (Diso) following the approach proposed in [29]. According to [30], the triggering front was modeled as r = √ 4πDt, where D is the scalar diffusivity corresponding to Diso/4. We only analyzed events with magnitude larger than the minimum magnitude of completeness that corresponds to Mc 1.6 ± 0.2 and was obtained by using the maximum curvature approach [31].
The resulting r − t plot is shown in Figure 8 and the inferred best Diso was 60 ± 37 m 2 /s. A similar value has been obtained for the 2009 L'Aquila sequence in Italy [32]. We note that, except for the very early days after the main event where aftershocks are expected to occur mainly due a redistribution of the stress on the fault planes, the evolution of the whole sequence is compatible with a pore pressure diffusion mechanism likely due to the volume deformation produced by the two main events. In particular, the Mw 6.0 and the Mw 5.6 occurred in correspondence with a triggering front characterized by the lower limit of the inferred diffusivity values (Figure 8).

Discussion
The TSS represents the largest seismic sequence affecting a continental extensional domain in Greece that has been widely monitored by modern geodetic techniques; thanks to the short satellite revisit time, InSAR measurements allowed for isolating the contribution from the major earthquakes of the sequence to study their interaction. Available geological data indicate that the northern sector of Thessaly represents a large seismic gap. This may be a direct consequence of the limited size of the faults (less than 20 km) and their intrinsic capability to originate earthquakes of small-to-moderate magnitude only [8]. TSS, which finally filled the gap, confirmed this hypothesis.
Our model shows that the seismic activity distributed along distinct segments of a previously unknown fault system composed of three fault plains (F1a, F1b, and F2) dipping NE for the Mw 6.3 and Mw 6.0 earthquakes, and of a small fault (F3) dipping SW for the Mw 5.6 ( Figure 6). However, the activation of some portion of the known Tyrnavos and Larissa fault during the aftershocks cannot be completely ruled out. The orientation of the retrieved fault plains is in good agreement with the published focal mechanisms by Greek and other international Institutes (Figure 1a), and is in accordance with the extensional tectonic regime of the area.
For the first two events, the vertical deformation component highlights strong subsidence phenomena reaching more than~40 cm. Moreover, to make the whole picture more complicated, the E-W component showed significant westward motion (~20 cm) ( Figure S7). This latter motion can be related to the movement of the F1a fault, which features a slight left lateral strike-slip component.
Extensive and deep ground cracks generated by the Mw 6.3 were observed close to river beds along with ejected liquefied material; other cracks were parallel to the Vlachogianni fault. Some field reports showed scattered surface alignments only in the Zarko region [1], which can hardly be associated with faults. However, the inferred slip distribution showed that the highest values were confined in the S-E area (Vlachogianni, Damasi, Zarko villages) and this could explain the observed damage distribution (Figures 1a and 6).
Our ∆CFF results (Figure 7) indicate that the seismic sequence developed in a sort of domino effect. Regarding the temporal evolution of the sequence, the delayed triggering of the Mw 6.0 earthquake can be explained by the distribution of the events that occurred earlier; indeed, that distribution envelops the patches of the fault that will break a feature, described in [33] as the encircling maneuver. A discrimination between the three main conceptual models-cascade-up, pre-slip and progressive localization [34,35]-for earthquake initiation would require a specific spatial-temporal analysis of the earthquakes sequence and stress perturbations between different events, which is beyond the scope of the present manuscript. We thus suggest that the most plausible initiation mechanism could be the recently proposed mixed stress loading process involving the cascade-triggering and aseismic-slip transient mechanisms [36].
As for the Mw 5.6 earthquake, the results depicted in Figure 7b,c show that the most relevant effect, in terms of static stress contribution, was due to the Mw 6.0 event. Furthermore, the delay in the triggering (about eight days from the Mw 6.3 and Mw 6.0 events) suggests that mechanisms other than the static stress change-in addition to a dynamic triggering contribution-took place. This is the case, for example, of pore fluid pressure diffusion as observed in other seismic sequences [37][38][39]. Indeed, a pore pressure increase can reduce the normal stress promoting the rupture. Moreover, fluids can also contribute to decreasing the fault friction [40]. The performed 1D diffusivity analysis indicates that almost all the sequence, and in particular the occurrence of the Mw 6.0 and Mw 5.6, are indeed compatible with a mechanism of pore pressure diffusion that was triggered by the deformation pattern produced by the previous events.
The combination of the lack of historical records of large earthquakes in this area, and the absence of mapped surface features produced by past faulting, make seismic hazard estimation difficult. Therefore, the identification of blind faults may significantly contribute to improve the seismic hazard analysis, particularly when dealing with moderate dip angles. In this context, the InSAR products have proven to be an extremely useful tool in helping to constrain the rupture characteristics.
Considering the possible presence of a multi-fault rupturing structure and the complexity of the tectonic framework of the studied area, we believe that the study of the seismogenic mechanism of the blind faulting deserves further investigation through additional data.

Conclusions
In this work, we studied the seismic sequence occurred in Thessaly (Greece) during the early days of March 2021. To this end, we used seismic and InSAR data to set up a geodetic model of the causative faulting system. Moreover, we performed Coulomb stress transfer and 1D pore pressure diffusivity analyses to investigate the space-time evolution of the sequence.
Our main findings are summarized as follows: • the Thessaly seismic sequence nucleated at shallow depths (<12 km) and is related to the activation of several blind, previously unknown faults; • the seismic sequence developed in a sort of domino effect involving a complex interaction among the normal faults within the activated crustal volume; • InSAR data and modeling are also extremely useful to constrain the rupture characteristics in the case of blind faults; and • the used approach can help improve our knowledge of the seismic potential of the Thessaly region and refine the associated seismic hazard.