A Deeper Look Into the 2021 Tyrnavos Earthquake Sequence (TES) Reveals Coseismic Breaching of an Unrecognized Large‐Scale Fault Relay Zone in Continental Greece

Large magnitude (Mw ∼ ≥6) earthquakes in extensional settings are often associated with simultaneous rupture of multiple normal faults as a result of static and/or dynamic stress transfer. Here, we report details of the coseismic breaching of a previously unrecognized large‐scale fault relay zone in central Greece, through three successive normal fault earthquakes of moderate magnitude (Mw 5.7–6.3) that occurred over a period of ∼10 days in March 2021. Specifically, joint analysis of InSAR, GNSS and seismological data, coupled with detailed field and digital fault mapping, reveals that the Tyrnavos Earthquake Sequence (TES) was accommodated at the northern end of a ∼100 km wide transfer structure, by faults largely unbroken during the Holocene. By contrast, the southern section of this relay zone appears to have accrued significant slip during Holocene. InSAR‐derived displacements agree with the loci of eight subtle, previously undetected, faults that accommodated coseismic and/or syn‐seismic normal fault slip during the TES. Kinematic modeling coupled with fault mapping suggests that all involved faults are interconnected at depth, with their conjugate fault‐intersections acting largely as barriers to coseismic rupture propagation. We also find that the TES mainshocks were characterized by unusually high (>6 MPa) stress‐drop values that scale inversely with rupture length and earthquake magnitude. These findings, collectively suggest that the TES propagated north‐westward to rupture increasingly stronger asperities at fault intersections, transferring slip between the tips of a well‐established, but previously unrecognized, relay structure. Fault relay zones may be prone to high stress‐drop earthquakes and associated elevated seismic hazard.


The Seismotectonics of the TES
The TES is a multi-fault sequence of three moderate magnitude (M w 6.3-5.7) normal-slip earthquakes that occurred in central Greece between 3 and 12 March 2021 (http://bbnet.gein.noa.gr/HL/index.php) ( Figure 2; Table S1 in Supporting Information S1). The earthquake sequence occurred in Thessaly, a region in continental Greece which is characterized by N-S extension and associated normal faulting ( Figure 2) (Caputo & Pavlides, 1993;Galanakis, 1997;Goldsworthy et al., 2002). The study area includes two main NW-SE trending Quaternary alluvial basins (the Karditsa and Larissa basins; shaded light yellow in Figure 3), separated by similarly aligned bedrock ridges (the Central Hills and Pelion Hills; Figure 3), forming a horst-and-graben type topography. A comparable submarine bathymetry, with successive horsts and grabens, is revealed though EMODnet bathymetry (https://www.emodnet-bathymetry.eu/) in the North Aegean Sea (Figure 2) (Papanikolaou et al., 2002). The study area is bounded in the west by the Pindos Range, that is underlain by ophiolites, flysch, dolomitic and pelagic limestone, and is separated from the Central Hills by the Quaternary Karditsa Basin (Figure 3). The Central Hills are underlain by augen gneiss, Permo-Triassic and Triassic-Jurassic limestone while the Pelion Hills, east of the Larissa Basin, comprise ophiolites, schist, gneiss, marbles and dolomites (Bornovas & Rondogianni-Tsiambaou, 1983).
The broader area of Thessaly and central Macedonia in Greece has hosted numerous large-magnitude normal faulting events in the recent past, including the 1941 M6.1 Larissa Earthquake, the 1954 M7 Sofades Earthquake, the 1957 M6.8 Velestino Earthquake, the 1980 M6.5 Almyros Earthquake and the 1995 M6.6 Grevena-Kozani Earthquake (Figure 2), mostly on faults with significant geomorphic expression (Ambraseys & Jackson, 1990). Interestingly, the 2021 TES was not accommodated on any of the previously known active normal faults ( Figure  S1 in Supporting Information S1), which are mostly characterized by spectacular scarps of Holocene age (Caputo & Helly, 2005;Caputo et al., 2004;Tsodoulos et al., 2016) (Figure 3). Instead, all three TES mainshocks (M w 6.3 Event 1 on 3 March, M w 6 Event 2 on 4 March and M w 5.7 Event 3 on 12 March; Figure 2) ruptured previously unrecognized faults that extend beneath hilly country (e.g., Chatzipetros et al., 2021;Ganas et al., 2021;Koukouvelas et al., 2021;Papadopoulos et al., 2021;. The epicentral region of the TES is well monitored instrumentally by an abundance of permanent/temporary seismic and GNSS stations (Figure 2), while InSAR allowed the quantification of the coseismic displacements associated with each of the main events (i.e., Ganas et al., 2021;Papadopoulos et al., 2021). Using a combination of seismological, geodetic and InSAR data, numerous studies attempted to kinematically model the rupture sequence, proposing either a two (Karakostas et al., 2021), a three (Ganas et al., 2021;Kassaras et al., 2022;Papadopoulos et al., 2021) or a four (De Novellis et al., 2021;Yang et al., 2022) fault source model. Despite the differences, all proposed models agree on the geometry and approximate location of the fault that produced the large (M w 6.3) March 3rd earthquake. This normal faulting event was unambiguously modeled by all above studies to have ruptured a shallow (∼40°) NE dipping plane of NW-SE strike for a length of ∼15 km. As the fault that accommodated the M w 6.3 event was characterized "blind" (e.g., Chatzipetros et al., 2021;Ganas et al., 2021;Papadopoulos et al., 2021;, field studies (i.e., Galanakis (Figure 1b; see Section 3.3 for details). The majority of the published studies also seem to agree that the TES ruptured overall immature faults (e.g., Karakostas et al., 2021;Kassaras et al., 2022;Koukouvelas et al., 2021) due to a domino effect of stress-transfer (De Novellis et al., 2021;Ganas et al., 2021;Kassaras et al., 2022) that underpins the westward propagation of the faults of the Larissa Basin (Chatzipetros et al., 2021;Koukouvelas et al., 2021).
Here, we (a) revisit all Quaternary faults in Thessaly and map a previously unrecognized large-scale (∼100 km) relay zone that connects normal faults in Volos region (south) to the normal faults in the region of Deskati (north) ( Figure 3); (b) present geomorphic evidence that the normal faults that ruptured during the TES, including the fault that produced the M w 6.3 earthquake, were not blind (Figure 4a) but elements of this well-established  Map-view (upper panels) and schematic cross-sections (lower panels) of normal faults accommodating coseismic slip on primary and secondary slip surfaces which appear interconnected at depth (a) or concurrent coseismic and syn-seismic slip (i.e., triggered slip on nearby faults which are not connected to the seismogenic zone) (b). (c) Examples of slip accumulation within an active rift during a notional earthquake sequence (Events 1-3) that spans different timeframes (hours = black; days = green; months = orange; see color-coded references for relevant examples). Pink shading implies distributed off-fault deformation. 1 = Bernard and Zollo (1989); 2 = Wallace (1980); 3 = Beanland et al. (1989); 4 = Delano et al. (2022); 5 = this study; 6 = Crone and Hailer (1991) Papanikolaou et al. (2002) and Sakellariou and Tsampouraki-Kraounaki (2019) using EMODENET bathymetric data. DF = Domokos Fault Zone (Palyvos et al., 2010; pink trace indicates demonstrable Holocene activity). Inset: the study area is located in continental Greece.

Figure 3.
A revised active fault map of Thessaly-Central Greece. Map superimposed on a hillshade illustrating the NW-SE-trending large-scale fault relay zone that traverses Thessaly to connect two E-W trending fault systems that extend between the areas of Volos (south) and Deskati (north). The faults are color-coded as per Figure 2, while the magenta color indicates trenched faults with dated Holocene activity (Caputo & Helly, 2005;Caputo et al., 2004;Tsodoulos et al., 2016). Numbers next to each fault/fault system correspond to fault IDs presented in Table 1 Table 1 for references). The main cities are indicated by black filled circles. The locations of the maps illustrated in Figure 4 are indicated by black boxes. relay zone; (c) combine instrumental (seismic, GNSS and InSAR) data sets with field investigations to propose a kinematic model in which the TES breached the northern ∼40 km of this relay zone, transferring normal displacements northwestward; and (d) show that as the rupture propagated from Event 1 to Event 3 (Figure 2), the crust accommodated increasingly higher stress drop to break progressively stronger asperities at normal fault intersections.

Data and Methodology
To study the TES and better understand its geotectonic context, we used data sets and methodologies that collectively span various temporal and spatial scales. Specifically, we undertook fieldwork and analysis of high resolution DEMs to map active faults in Thessaly (Section 2.1), analyzed ∼70 days of seismicity associated with the TES (3 March 2021 until 12 May 2021) (Section 2.2.), analyzed all available InSAR and GNSS data sets that recorded TES displacements (Section 2.3 and 2.4) and kinematically modeled the deformation produced during (or soon after) the three mainshocks (Section 2.5). Subsequently, the TES is discussed within the broader geotectonic setting revealed following our revised normal fault mapping (Sections 3 and 4).

Quaternary Normal Faulting: Field and Digital Mapping in Thessaly
The previously recognized and/or studied active faults within the TES epicentral region are the Tyrnavos Fault (TF), the Gyrtoni Fault (GF), the Rodia Fault (RF), the Elassona Fault (EF), the Larissa Fault (LF) and the Vlachogianni-Mesochori Faults (VM) (Caputo, 1990(Caputo, , 1995Koukouvelas et al., 2021) (Figure 3). Their published fault length ranges from ∼10 to 15 km while youthful landscape features offset by these faults indicate that they were active during the Holocene. Indeed, fault trenching revealed that at least eight paleoearthquakes were accommodated on three of these faults (TF, RF, and GF) during the last ∼10 ka, with average recurrence intervals ranging between 2 and 4.7 ka, while the most recent event on each fault is older than 2 ka (Caputo & Helly, 2005;Caputo et al., 2004;Tsodoulos et al., 2016). These faults in Thessaly, and associated published information, provide the basis (and inspiration) for the work presented here; nevertheless, the vast majority of these faults are only mapped partly (i.e., traces or segments are missing) or approximately (they are represented by straight lines), while numerous other faults in the area remained to date unmapped ( Figure S1 in Supporting Information S1). Revisiting, thus, the faults in Thessaly appeared to be a pressing need, not only for better understanding this and future strong earthquakes, but also for better quantifying the seismic hazard in Thessaly.
Here, we use a 5 m DEM (from the Greek Cadastre Agency) with a vertical resolution of ∼2 m (fault scarps ≥2 m are here resolved), field mapping and QGIS software (version 3.16), to map active fault traces in Thessaly (Figures 2 and 3). Our work allows the refinement of the surface trace of ∼10 previously known active normal faults (Table 1 and Figure S1 in Supporting Information S1) and the identification of ∼40 new active faults and/ or fault systems, including those that ruptured during the TES ( Figure 3 and Table 1). In contrast, five faults within the Larissa Basin which were previously classified as active (the Dimitra, Asmaki, Kastri, Mavrommati and Mitropoli faults; Caputo, 1990Caputo, , 1995 Figure S1 in Supporting Information S1) have not been included in the map of Figure 3 as their active traces could not be confirmed using geomorphological criteria. Figure S1 in Supporting Information S1 provides a useful comparison between the previously mapped active normal faults (upper) or lineaments (lower) with the active fault traces mapped by this study.
Lineaments mappable on 5 m DEM within the epicentral area are commonly represented in the field by active faults. The normal faults presented here are classified according to their relative activity as "active faults with Holocene scarps" and "active faults with pre-Holocene scarps," broadly following the classification provided by Nicol et al. (2020) and Faure Walker et al. (2021). Direct dating of offset Holocene deposits is available only for 4 faults in Thessaly (RF, GF, TF, and Domokos-Ekkara Fault; Caputo & Helly, 2005;Caputo et al., 2004;Palyvos et al., 2010;Tsodoulos et al., 2016), three of which lie in the immediate vicinity of the epicenter (see magenta lines in Figure 3). For the remaining faults there has been no direct dating of offset landforms, and their relative (pre-Holocene/Holocene) activity is based on the following geomorphic criteria. Those characterized as "active faults with pre-Holocene scarps" (the green traces in Figures 2 and 3) meet at least two, and often more, of the following characteristics (see Table 1 for nomenclature): have significant offset of materials considered late Quaternary (ST), are characterized by the presence of facetted spurs (FS) or stream displacement (LS), form linear range fronts (LR), form an elongate lineament clearly distinct from basement texture (L), significantly displace or truncate ridge lines (RT), often with linear scarps close to ridge crests, and/or display elongated changes in slope (SC) that cross-cut bedrock fabric. In addition to any two of the above characteristics, the "active faults with Holocene scarps" (the red traces in Figures 2 and 3) must also meet at least one of the following criteria: have sharp surface traces (SS) on materials of inferred Holocene age and/or lie at the margins of areas of active deposition (DC). For the former faults (i.e., the green faults), even if recent activity cannot be demonstrated, they should be considered active, particularly when their orientation and dip direction are comparable to faults that are demonstrably active. These 10 geomorphic criteria are used to characterize the activity on each fault trace in Table 1 (see column "Geomorphic Criteria").
For the purposes of this study, we assign individual traces to faults and/or fault systems ( Figure S2 in Supporting Information S1). The criteria for this classification are purely geometric: individual traces that appear to be segments of a larger structure are treated as single faults and their length corresponds to the total length of their along-strike extent (i.e., Fault IDs: 4, 5, 10, 11, 26, 28, 40, 41; Figure S2 in Supporting Information S1). Neighboring normal fault strands, with mostly antithetic dips which are located in close proximity to one another (<5 km) and possibly intersect at depth, are also considered here parts of a single system (i.e., Fault IDs: 3, 12, 18, 21, 37, 39; Figure S2 in Supporting Information S1). In a few cases, where the association of a fault with other neighboring faults is unclear (i.e., Fault IDs 13, 19, 20, 23, 33; Figure S2 in Supporting Information S1), faults appear as isolated structures. Although this classification helps in the remaining discussion, especially when we examine the TES within the context of the larger transfer structure that hosted the sequence, the precise inclusion (exclusion) of a normal fault trace in (from) a system, does not impact on the results/conclusion drawn here.

Seismological Data
The bulk of the seismological data are recorded by the Hellenic Unified Seismic Network (HUSN; http://bbnet. gein.noa.gr/HL/), which is the permanent network that has been monitoring seismic activity in Greece since 2008. Data from a total of 49 HUSN stations that are located within a distance of 200 km from the epicentral TES were included in our analysis, as well as six temporary stations (TYR1-TYR6) that were installed shortly after the occurrence of the two first main events ( Figure 2). The earthquakes analyzed here span the period between 3 March and 12 May 2021. The recorded waveforms have been manually picked for P-and S-phase arrival times by the staff of the National Observatory of Athens, Institute of Geodynamics (1997), however, we also inspected and adjusted the picks whenever necessary.  To improve the location of the earthquake sequence, we derived a minimum 1D velocity model with station delays for our study area by using the established Velest software package (Kissling, 1995). Details on the constraints used to derive this model as well as a number of sensitivity tests conducted to test its robustness, may be found in Supporting Information S1 (Figures S3-S5 and Text S1 in Supporting Information S1). The newly derived minimum 1D velocity model and the nonlinear probabilistic algorithm NonLinLoc (Lomax et al., 2009) were utilized to obtain absolute locations for the events of the TES. The OctTree search algorithm (Lomax & Curtis, 2001) was employed for reconstructing the posterior probability function for the purpose of sampling the solution space of each individual event. We further constrained the locations by also utilizing the equal differential time likelihood function (Font et al., 2004) that is calculated as the difference of residuals at station pairs. Uncertainties for each absolute location are the variances provided by NonLinLoc (Maleki et al., 2013). Figure S7 in Supporting Information S1 shows the distribution of the horizontal and vertical uncertainties where their mean values are 0.85 km (±0.33 km) and 1.24 km (±0.55 km), respectively. These small uncertainties reflect the fact that the majority of the events in the sequence were recorded either by the nearby permanent station TYRN, or by the temporary stations installed after the mainshock that surrounded the rupture zone.
In an effort to test whether the absolute locations could be further improved, we also performed a relocation using the double-difference algorithm HypoDD (Waldhauser & Ellsworth, 2000). Details on the relocation process may be found in Text S1 and Table S2 of Supporting Information S1. Figure S6a in Supporting Information S1 shows a map of the relative locations obtained by HypoDD which can be compared with the distribution of the absolute locations ( Figure 5a). As the NonLinLoc locations are already well-constrained and do not differ significantly from the relative locations, and also because they represent the total number of aftershocks within the study period rather than a fraction relocated by HypoDD (2870 vs. 2388), we prefer to discuss and interpret the spatiotemporal earthquake patterns produced by the absolute locations ( Figure 5). Nevertheless, in Supporting Information S1 we also include depth cross-sections of the relative locations ( Figure S9 in Supporting Information S1).

InSAR Data
The earthquake sequence was covered by 6-day repeat Sentinel-1 SAR measurements (C-band, 5.6 cm radar wavelength) from four different viewing directions: two down-ENE looking (ascending) looks of different incidence and two down-WNW looking (descending) looks with different incidence angles (Table S3 in Supporting Information S1). The scheduled acquisition time of one ascending orbit realized a SAR measurement between the first and the second mainshocks on 3 and 4 March 2021, such that with these interferograms we can separately analyze the displacements of both earthquakes (Table S3 and Figure S14 in Supporting Information S1). The 3 March acquisition took place approximately 6 hr after the first large earthquake of the sequence. All other SAR acquisitions only allow the measuring of the cumulative displacement of the first two large earthquakes. The surface displacements caused by the third larger shock on 12 March 2021 can be measured separately from all covering orbits ( Figure S14 in Supporting Information S1). We use the line-of-sight displacement measurement in the kinematic fault modeling (see Section 2.5).
With acquisitions from down-ENE and down-WNW look directions we can extract estimates of the eastward and vertical surface motions (method details provided in Supporting Information S1) of the combined signal of the first two earthquakes and the third large earthquake ( Figure 6). The resulting displacement signals are dominated by smooth and long wavelength displacements. To enhance small scale and sudden changes in the displacement measurements, we derive also maps of differential displacements for the eastward and vertical components of motion, calculating the difference of a pixel value from its second-to-the-east neighbor pixel (Figures 6c and 6d).

Geodetic Data
Data from six permanent GNSS stations distributed within and around the epicentral area of the TES were analyzed ( Figure 2). This data set primarily covers far-field effects at distances ∼15-55 km from the first mainshock (Figures 2 and 6). Data were collected on a 24 hr basis, with a sampling interval of 30 s between 1 January 2021 and 15 March 2021, except for the ELAS station, for which data were available only until 5 March 2021 (Figures S11-S13 and Table S4 in Supporting Information S1). Kinematic (30-s sampled) and daily coordinates for the horizontal and vertical components were calculated using the Bernese GNSS Software v5.2 (Dach et al., 2015). For each station, coordinates were estimated with respect to the GNSS station PAT0 (∼150 km to  southwest; EUREF network) whose coordinates were fixed to its ITRF2014 coordinate (Altamimi et al., 2016). The operational final CODE products (Dach et al., 2020) were introduced as reference products. This relative approach allowed resolution of the carrier phase ambiguities. GNSS coordinate time-series, after 3-sigma outliers were removed, are shown in Figures S11-S13 of Supporting Information S1, in which no significant postseismic signal is evident after the first event.
Coseismic displacements of a few centimeters are recorded on stations ELAS and KLOK for the first two main events (3 and 4 March), while there is no resolvable offset during the third mainshock (12 March). Therefore, coseismic offsets and associated uncertainties were calculated only for the first two mainshocks. Coseismic displacements were derived from the kinematic GNSS coordinates, using step (Heaviside) functions, and taking into account the exact time of the earthquakes within the corresponding day. Figures S11-S13 in Supporting Information S1 indicate consistency between the displacements derived from 30-s recordings and daily coordinates. For both events displacements up to 40 mm are recorded in the north component at stations ELAS and KLOK while the southeastern tip of the first rupture is recorded by station LARI (Figures 2 and 6 and Table S4 in Supporting Information S1). Overall, the displacement pattern indicates NE-SW extension, consistent with the focal mechanisms of the mainshocks and the strike of the mapped faults ( Figure 2). Vertical coseismic displacements, as expected for far-field stations, are minimal and noisy.

Kinematic Earthquake Source Modeling
For the kinematic finite-fault modeling of the three larger earthquakes in the TES, we combine coseismic InSAR and GNSS surface displacements with teleseismic waveforms in various schemes. For example, only for the first earthquake all three different data sets are included in the analysis. During the 2nd event (4 March), a strong (M w 7) earthquake along the Kermadec Ridge in the south Pacific (USGS, 2022) masks the waveforms of the ∼M w 6 TES event. Thus, for the source analysis of this earthquake, we only use InSAR and GNSS data. For the first and the second earthquake source analysis, we used only partial interferograms from the descending look direction to mask out the signal of the second and first earthquake, respectively. Regarding the 3rd event (12 March), its epicenter was far off the available GNSS stations such that they showed no coseismic displacements (see Section 2.4) and have, thus, not been used in the source analysis of the 3rd event.
We downsample the high-resolution InSAR data for efficient modeling with the irregular quadtree algorithm (Jónsson, 2002) and account for data errors through variance-covariance weighting as detailed in Sudhaus and Sigurjón (2009). The available seismic data from the IRIS and Geofon data services (USGS, GEOFON) are used, but spatially decimated for sufficient density and sufficient coverage in azimuth and distance around the epicenters (Figures S15f and S17f in Supporting Information S1). We conduct full-waveform fits of low-frequency P-and S-waves (Figures S15e and S17e in Supporting Information S1).  The preferred earthquake slip representation is a low-parametric rectangular plane with uniform slip. We optimize simultaneously for the location of the rupture plane, its depth and dimensions, the rupture mechanisms (strike, dip, rake) and the seismic moment released (Table 2). For the first (3 March) mainshock only, we have also predicted, using waveform data, a nucleation time and point from which the rupture initiates at a constant velocity (2,000 m/s). The source is embedded in a horizontally layered elastic medium ( Figure S3 in Supporting Information S1). To optimize the source for the best fit to weighted data, we employ a randomized direct search algorithm with Bayesian bootstrapping. We use here the implementations of the Pyrocko software project (Heimann et al., , 2018(Heimann et al., , 2019Isken et al., 2017). Details of the kinematic modeling and associated uncertainties are extensively discussed in Text S3 of Supporting Information S1.

The TES Ruptures the Tip of a Previously Unrecognized Large-Scale Fault Relay Zone
A total of 49 faults and/or fault systems have been mapped (Table 1; Figure 3 and Figure S2 in Supporting Information S1). Twenty-two (22) appear to have ruptured during the Holocene (red-colored faults in Figure 3), twenty-two (22) are characterized by pre-Holocene scarps (green-colored faults in Figure 3), while another five (5) have segments that ruptured both since and prior to the Holocene (Table 1). We increased the available fault trace information for most of the known active faults ( Figure S1 in Supporting Information S1) and identified hundreds of new active normal fault traces that we grouped into ∼40 new fault systems (Table 1 and Figures S1 and S2 in Supporting Information S1).
An example of a "revised" fault length is illustrated in Figure 4b, where the published length of the Gyrtoni Fault of ∼13 km (Tsodoulos et al., 2016) is now revised to ∼18 km. Another striking example of a revised fault length is that of a major seismogenic source within the Larissa Basin, the Tyrnavos Fault, with an assigned length of ∼12 km (Caputo, 1990;Caputo et al., 2004) that has now been revised to 21 km ( Figure S1 in Supporting Information S1; Figure 3 and Table 1). Detailed examples of newly identified active faults with Holocene (red) or pre-Holocene (green) scarps, are illustrated in Figure 4. Displacements measured across faults that traverse youthful landscape prove that these lineaments are active ( Figure S10 in Supporting Information S1). The aggregated maximum displacement on each fault system is presented in Table 1 (see Figure S2 in Supporting Information S1 for the locality of each displacement profile). Collectively, mapped fault lengths range from 5 to 95 km while maximum cumulative displacements range from ∼10 to >800 m ( Table 1). The reviewed fault lengths/ displacements presented here should be still considered as minima, because (a) some faults intersect the coastline and extend offshore, (b) some others are partially (Domeniko FS; ID = 36) or entirely (i.e., the Larissa Fault; ID = 19) concealed beneath basin infill materials; (c) fault scarps of <2 m are not resolvable on the available DEM and (d) because displacements of any size may have been reduced due to burial/sedimentation.
Revised fault mapping reveals a previously unrecognized, large-scale (>100 km) extensional transfer zone that extends from Volos, in the south, to Deskati, in the north (Figure 3). A transfer zone commonly includes the master faults and a number of intervening relay structures (Fossen & Rotevatn, 2016), which are typically characterized by a wider range of orientations (Mercuri et al., 2020;Peacock & Sanderson, 1994). All 49 mapped fault systems (which herein we will refer to as "faults") are interpreted to be elements of this large-scale transfer structure (Figure 3). The southern and northern ends of the transfer zone are characterized by large (>25 km) E-W trending normal faults that in their vast majority appear to have accommodated slip during the Holocene (Faults 2-6 in the south and 25, 40-41, 44 in the north). These E-W faults, which are the master strands of the system, are hard-linked via a series of smaller-sized (<25 km) normal faults that extend for ∼70 km with a NW-SE strike, transferring displacement between the master E-W striking faults (Figures 3 and 7a). Displacement transfer along the relay zone is indeed suggested by four system-perpendicular transects in Figure 7a (for locations see Figure  S2 in Supporting Information S1), where uniform displacement of ∼800 ± 200 m have been accommodated since the inception of faulting along the entire extent of the system. The NW-SE trending faults, that herein we will refer to as "relay faults," are mostly characterized by pre-Holocene scarps (Faults 21-24, 26, 36-37), have subtle but measurable geomorphic signature ( Figure S10 in Supporting Information S1 and Table 1) and have variable orientations, typical for faults in relay zones (Mercuri et al., 2020). Interestingly, although the central section of this relay zone includes faults with pre-Holocene scarps only ("green faults"), its northern and southern tips faults appear to have ruptured during the Holocene (Faults 10.11,22,26), similarly to the E-W trending master faults, pointing to a well-established transfer zone, a property that will be discussed further in Section 4.2.

Coseismic and Syn-Seismic Displacements During the TES
Our geological mapping reveals numerous NW-SE trending faults in the northern domain of the relay zone where the three TES mainshocks occurred (Figures 2 and 3). The Zarkos Fault System (ID = 21; Table 1 and Figure 3), for example, is the most likely candidate for the M w 6.3 March 3rd event. Field observations along its southern section indicate fault gouge and tensional cracks at the ground surface (Figure 4a), in broad agreement with findings in Ganas et al. (2021). The following two mainshocks did not produce slip that could be resolved conventionally in the field; however, a number of candidate structures with subtle but clear geomorphic signature extend above the inferred rupture planes (see faults with IDs 22, 36, 37, 25, and 40 in Figures 3  and 4). To explore which of these faults may have ruptured during these two events in Thessaly, and validate rupture of the Zarkos FS, we analyzed the available InSAR and GNSS data sets that, collectively, sample surface displacements down to a few millimeters ( Figure 6).
For the first mainshock on March 3rd, the InSAR-derived surface displacement maps show significant elliptical subsidence with associated components of horizontal motion ( Figure S14a in Supporting Information S1). The subsidence has a northwest-southeast orientation and extends for ∼25 km between the Kalyvia-Farkadona FS (ID 22), in the southwest, and the Tyrnavos Fault (ID = 28), in the northeast ( Figure S14a in Supporting Information S1). The maximum recorded subsidence is about 40 cm and it is observed directly above the Zarkos FS, with asymmetric displacement gradients perpendicular to the NW-SE striking structures (Figures 6c and 6d). Southwest of the subsidence signal, the displacement gradient is the highest and indicates that the earthquake here coseismically ruptured one or more segments of the Zarkos FS (Figures 6a-6d). The exact location of the surface rupture is not easy to discern in the wrapped interferograms (i.e., Figure S14 in Supporting Information S1) because of the generally high phase gradient causing spatially high density of wrapped phase cycles and partly phase decorrelation (for details see Text S2 in Supporting Information S1). The surface rupture may be more clearly inferred by the localized, probably post-seismic, motion captured in a subsequent InSAR measurement that also includes displacements of the second mainshock on 4 March ( Figure S14b in Supporting Infor mation S1). There, coseismic rupture of the southern section of the fault is recorded, while its northern extent appears to have ruptured in the subsurface only. This finding, together with our field observations and mapping, clearly suggests that the strongest of the three TES mainshocks was accommodated by the Zarkos FS.
The InSAR displacement signal due to the second mainshock (M w 6) on 4 March occurred about 16 km northwest of the first large earthquake ( Figure S14b in Supporting Information S1). Here, surface displacements are also broadly consistent with subsidence, but with displacement values being approximately three times lower compared to the first mainshock. The displacement gradients of this second mainshock are rather symmetric to the signal centers and show no indication of surface rupture (Figures S14b-S14d in Supporting Information S1). Similarly, the third TES mainshock (12 March) occurred ∼9 km to the northwest of the second (4 March) and produced about 50% less subsidence signal compared to the second mainshock, with small and symmetric displacement gradients (Figure 6e and 6f and Figures S14c, S14e, and S14f in Supporting Information S1).
Our InSAR analysis, in addition to the above large-scale well-resolved displacement signals, reveals smaller, step-like, slip gradients that form km-long lineaments along existing normal faults (Figures 6c and 6d and associated profiles). These displacements are significantly smaller (<5 cm) compared to those produced coseismically (up to 40 cm), and are mostly formed at an angle to the main displacement gradient. Such small-scale displacements are illustrated in the differential-pixel maps of Figures 6c and 6d and measured on strike-perpendicular  (Table S5 in Supporting Information S1). (b) Tyrnavos Earthquake Sequence (TES) versus Global Data Sets: Displacement (D)-length (L) data from the relay zone in Thessaly (red circles), global displacement-length data for ancient faults (gray circles; Nicol et al., 2017 and references therein) and single-event displacement rupture lengths for historical earthquakes (orange circles; Wesnousky, 2008). The red filled circles with the black outline are the faults that accommodated the three TES mainshocks in March 2021 while the three magenta circles are the faults in Thessaly with demonstrable repeated Holocene activity (see also magenta lines in Figure 3). Least squares lines of best fit and R 2 values are also indicated for each data set.
profiles in the following cases: (a) Along the Kalyvia-Farkadona FS (ID 22) where a >5 km long section of the fault records ∼1 cm slip down-to-the-east (Profile F in Figures 6c and 6d). (b) At the southern tip of the Zarkos FS (ID = 21), where a >3 km long section slips also down-to-the-east by ∼1 cm in both components (Profile X in Figures 6c and 6d). (c) Along the Vlachogianni-Pretorio FS (ID = 26) where an >8 km fault section, accommodates shallow slip with alternating kinematics (i.e., throw directions vary from down-to-the-east to down-to-the-west along the Profiles V1 to V3 in Figures 6c and 6d). Indeed, careful examination of these profiles (V1-V3) reveals that the northwestern-end of the activated fault section slipped (for ∼3 km) with an east-down component of slip, as it has also been recorded by Koukouvelas et al. (2021) using UAV measurements, while further southeast, the slip directions are completely reversed, with a steady east-up component of slip (up to >1 cm in both the vertical and east components). (d) On a strand of the Zarkos FS which is sub-parallel to the one that slipped coseismically, we record, for about 1 km, a small down-to-the-east displacement that extends across where it splays into multiple strands (Profile Z in Figures 6a-6d). (e) Along the junction of the Elassona-Mavrelli F.S. (ID = 40) and the Vlachogianni-Pretorio FS (ID = 26) where small down-to-the-south displacements were recorded in response to the third mainshock (Figures 6e and 6f).
Because these displacements are small and localized along shallow depths only, and also because they form at an angle to the main displacement gradient, we do not attribute them to the earthquake rupture process itself but to syn-seismic (or sympathetic) rupture ( Figure 1b) in response to ground shaking (e.g., gravitational collapse). It is very likely that these features, which were all formed within the first 6 hr of the first earthquake (Text S2 in Supporting Information S1), were produced during the 3 March mainshock without, however, being genetically related to its seismogenic zone.

Comparison Between Geological and Remotely Sensed Displacements (InSAR)
It is intriguing to observe that nearly all ground surface displacements inferred by InSAR align well with active faults, chiefly mapped by this study (Figure 6). Specifically, here we find that the only coseismic displacement recorded by InSAR to have ruptured the ground surface during the TES (see solid blue line in Figure 6 and in the wrapped interferogram of Figures S14a, S14c, S14d in Supporting Information S1), extends along a section of the Zarkos FS that has clear geomorphic signature (Figure 4a). Thus, the fault on which the mainshock took place is not blind-as many have suggested (e.g., Chatzipetros et al., 2021;De Novellis et al., 2021;Ganas et al., 2021;Papadopoulos et al., 2021; but, instead, includes traces with a clear history of preceding surface ruptures prior to the March 3rd event (Figure 4a).
Places where InSAR suggests rupture at depth (see dashed blue lines in Figure 6) align well with mapped fault traces at the ground surface (i.e., faults with IDs 21, 36, and 40). In circumstances where InSAR suggests syn-seismic displacements (purple lines in Figure 6), in all cases but one, they align well with mapped fault traces at the ground surface (IDs 21, 22, and 26). The one case where InSAR-derived displacements do not extend along mapped fault traces relates to the 3rd event, where syn-seismic displacements extend, with a curved map-view projection, across the intersection of the Elassona-Mavreli and Vlachogianni-Pretorio fault systems (IDs 40 and 26; Figures 6e and 6f). Interestingly, this displacement pattern mimics the map-view geometry of the northern-end of the Kalyvia-Farkadona FS (ID 22), likely indicating the presence of an undetected fault ( Figure 6).
The overall very good agreement between remotely recorded ground displacements and mapped active faults reveals the high potential of InSAR techniques in resolving small-scale coseismic and syn-seismic displacements along individual active faults but also in revealing transient kinematics (i.e., reversal of slip along slip surfaces like the one recorded along Fault 26 and which has also been documented by InSAR during the 2019 Ridgecrest earthquakes; Xu et al., 2020). In reverse, where InSAR suggests slip in areas with no mapped faults (i.e., see junction between faults 40 and 26 in Figure 6), there is implication of a concealed/undetected active structure, highlighting the importance of InSAR in detecting previously unrecognized active faults with subtle or concealed traces.

Seismicity Patterns During the TES
The characteristics of the three normal-fault mainshocks are presented in Table S1 of Supporting Information S1 and are discussed in detail in various preceding publications (i.e., Karakostas et al., 2021;Kassaras et al., 2022).
Here, we focus on the characteristics of the aftershock sequence and the inferred relationships between the distri-bution of aftershocks and mapped faults. Figure 5a shows a map of the distribution of 2870 located events, forming a ∼40 km long NW-SE trending zone of seismicity, in accordance with the strike of focal mechanisms of the mainshocks (NOA) (Figure 2). The hypocentral depth distribution shows that most of these events occurred at depths ranging between 5 and 10 km (average ∼7 km), with their number diminishing significantly at greater depths ( Figure 5a). The earthquake magnitude of the aftershock sequence ranges from 0.1 to 5.2.
To investigate the spatiotemporal distribution of the seismicity patterns within the TES, we constructed depth cross-sections (Figure 5b and Figure S6b in Supporting Information S1) along selected profiles that extend across the broader epicentral region (Figure 5a). Examination of the depth profiles, which are color-coded according to their timing since the first mainshock on March 3rd (Figure 5b), reveals a broad migration of the epicenters toward the NW (i.e., Profiles G, F, E, D are dominated by red color while profiles A-C by blue). This NW migration of the aftershocks, which is also clearly illustrated along the H-H′ Profile, broadly reflects the spatial and temporal distribution of the mainshocks. On the same profile, we also observe a 10 km long section of the crust immediately above the 3 March epicenter to be depleted from aftershocks (see gray box in Figure 5b). Microseismicity outlining sections of a fault-plane has been found to characterize both, locked asperities that may rupture during future events (Chlieh et al., 2011;Moreno et al., 2010;Saltogianni et al., 2020;Sippl et al., 2021) and/or fault plane patches that have recently accommodated one (or more) large-magnitude earthquakes (Andinisari et al., 2020;Konstantinou, 2017Konstantinou, , 2018. The relationship between mapped faults and aftershocks, as depicted solely by the distribution of microearthquakes, is less clear ( Figure 5). Aftershocks appear to be widely distributed in between fault strands, possibly due to complex fault networks mapped within the epicentral area ( Figure 3). Indeed, during the 3 March event significant slip was accommodated by the southern tip of the Zarkos FS, where at least five fault strands extend across a strike distance of <2.5 km (Figures 3, 4a, and 6). The Profile F-F′, that extends proximal to the 3 March epicenter, reveals an overall NE-dipping plane with most of the aftershock activity occurring below the inferred plane (Figure 5b). Similarly, the Profiles D-D′ and B-B′ that extend proximal to the 4 and 12 March epicenters, reveal weak southwestern and southern dipping planes, respectively (Figure 5b). The majority of the aftershocks of the 4 March event extend also below the inferred rupture plane (Figure 5b). This is in support of geological and geodetic data (see Sections 3.1-3.3) which, collectively, indicate rupture of the SW-dipping Fault 36 and a south-dipping strand of Fault 40 during the last two TES mainshocks (Figures 3 and 6; Table 1). The remaining cross-sections (Profiles A, C, E, and G in Figure S6b of Supporting Information S1) are inconclusive.

Kinematic Model for the TES
After considering various geometric configuration and testing different kinematic scenarios, we derive our preferred kinematic model for each of the three mainshocks (Figure 8a, Text S4 and Figures S15-17 in Supporting Information S1). All three events are modeled as predominantly normal faults with a minor left-lateral component (−66°, −66° and −76° for the 1st, 2nd, and 3rd TES mainshocks, respectively), in agreement with the documented N-S extension in the region and the overall NW-SE trend of the mapped faults (Caputo & Pavlides, 1993;Konstantinou et al., 2017; this study). The first mainshock (3 March) resulted from slip on a shallow dipping plane (34°NE), which extends for ∼14 km along-strike and for a down-dip width of ∼10 km (Figures 8a and 8b and Table 2). The rupture plane has a NW strike (29°W) and aligns well with the mapped Zarkos FS (ID = 21). The small mismatch (<1 km) observed between the surface projection of the modeled rupture plane and the mapped Zarkos FS (Figure 8), indicates steepening of the fault-dip toward the surface (in the upper <2 km). InSAR displacement data (Figures 6a and 6b) show that significant slip (∼25 cm) reached the ground-surface along the southern section of the rupture plane, however, our low-parametric rupture model could not resolve this slip. In the northern half of the modeled rupture, significant slip occurs at depths greater than 1.5 km, which is the upper extent of our optimal modeled rupture, while the average slip of this event is modeled at 72 cm ( Table 2).
The 2nd TES mainshock is modeled to be blind and to extend at crustal depths between 8 and 5 km, significantly deeper than the 1st mainshock (Figures 8a and 8b and Table 2). The best fit model is compatible with uniform normal slip of 90 cm along a 9 km long (strike of 131°), southwest-dipping (∼46°) rupture plane. The strike of this second rupture plane is rotated slightly counterclockwise with respect to first rupture plane (Figure 8a), aligning well with the strike of mapped faults (Figure 3). Indeed, a linear upward projection of the rupture plane aligns well, at its southeastern end, with the mapped Domeniko FS (ID = 36 in Figure 8a). The extension of the rupture plane toward the northwest suggests that a section of the Domeniko FS may be concealed beneath recent basinal sediments. A northeast-dipping rupture plane, as proposed by Kassaras et al. (2022) may also be possible, but it would require a more NW-SE strike ( Figure S19 in Supporting Information S1), which is incompatible  Table 1). Fault sections are color-coded according to rupture style (coseismic, syn-seismic, etc). A-A′ and B-B′ mark the locations of depth profiles across key faults (and their intersections), schematically illustrating the geometric relationships of faults/ruptures associated with the first two mainshocks (b) and the latter two mainshocks (c), respectively. Fault and rupture plane geometries/ dimensions are tailored to the values of the preferred kinematic models (a) and to field measurements (i.e., for faults with IDs 28,26,26a,22). Syn-seismic slip information on the Mesochori Fault (ID = 26a) comes from Koukouvelas et al. (2021). Topography is derived from the DEM. Color coding on faults of both profiles follow that of (a). Close-up view of normal fault intersection in (b) is modified from McCormack and McClay (2018).
with the strike of the mapped active faults in the region (it would cut across and violate the grain of the mapped structures; Figure 3).
The rupture model of the 3rd mainshock shows further counterclockwise rotation of the fault strike, in agreement with the mapped structures of the Elassona-Mavreli Fault System (ID = 40 in Table 2 and Figure 8a). The comparatively steep dip of this plane (∼58°SE) brings its surface projection to align well with geologically mapped faults as well as with traces mapped by InSAR only, as small syn-seismic steps in the displacement maps ( Figures 6 and 8a and Figures S14e and S14f in Supporting Information S1). The modeled rupture occurred at depths between 5 and 2 km (Figures 8a and 8c and Table 2), which is slightly shallower than the reported hypocenter depth by NOA (http://bbnet.gein.noa.gr/HL/index.php).

Fault Intersections and Their Role During the TES
In order to visualize the subsurface geometries of the faults involved in the TES, we used geometric constraints from the optimum kinematic modeling (Section 4.1; Figure 8a) and the geology (Section 3.1) to derive two depth profiles across the inferred intersections of the faults that ruptured (Figures 8b and 8c). Profile A traverses the Zarkos (ID = 21) and Domeniko (ID = 36) fault intersection, while Profile B is drawn across the intersection between the Domeniko Fault (ID = 36) and the strand of the Elassona Fault (ID = 40) that ruptured during the 3rd mainshock. From these profiles, we derive the following observations: (a) The M w 6.3 3 March event nucleates on the Zarkos Fault at depths (∼7 km) clearly greater than its antithetic intersection with the Domeniko Fault (∼4 km) and extends to the ground surface. (b) Slip on the M w 6 event that follows a few hours later (4 March), is restricted to depths beneath this intersection which, in this case, appears to have acted as barrier to rupture propagation ( Figure 8b). The microseismicity that precedes this 2nd mainshock, and is mostly associated with the 1st event, clusters at the inferred intersection of the two rupture zones, while the early aftershocks of the 2nd mainshock align favorably with the location of its rupture plane at depth ( Figure S18 in Supporting Information S1). has not accommodated resolvable coseismic slip during the 3 March event although its northern section intersects the Zarkos Fault proximal to where this event nucleated. This may reflect low stress being currently stored on the Tyrnavos Fault, which is known to have accommodated three M > 6 earthquakes in Holocene (Caputo et al., 2004) and may, therefore, experience currently a quieter phase in its seismic cycle (Nicol et al., 2009). Although intriguing, exploring further these rupture relationships using Coulomb stress modeling, may be problematic as we lack significant information on the pre-stress conditions (Mildon et al., 2019) on all remaining faults. Collectively, the above points (a-e) suggest that both synthetic and antithetic normal fault intersections in the broader TES region, may have acted as barriers to rupture propagation but also as loci to rupture initiation (Aki, 1989;Walters et al., 2018;Zhang et al., 1991).
To explore further the likely control of fault intersections on the size and timing of the TES mainshocks, here we combine the rupture area (A) obtained from our kinematic model for the three mainshocks (Table S2 in Supporting Information S1) with the Global CMT seismic moment (M 0 ) of each event to infer stress drop (Figure 9; for details see Text S4 in Supporting Information S1). Together with the TES stress drop values, in the graph of Figure 9 we have also plotted rupture areas and seismic moments of 53 other events (Mw 4.5-7.5) that occurred in the Mediterranean between 1976 and 2013 (for more details see Konstantinou, 2014). Most of these events exhibit stress drop values between 1 and 6 MPa, with the exceptions of the 1995 Kozani (Greece) and 1997 Umbria (Italy) earthquakes whose stress drop was notably higher (indicated in Figure 9). When compared to this pattern, all TES events  Konstantinou [2014]). Gray circles indicate dip-slip events and gray squares strike-slip ones (note that the largest of these earthquakes is the 1999 Mw 7.6 Izmit, Turkey, earthquake). Colored circles indicate dip-slip (normal faulting) events with stress drops higher than 6 MPa, while filled symbols represent the three largest events of the 2021 Tyrnavos Earthquake Sequence analyzed in this work.
are associated with high stress drop values. The sequence initiated with the March 3rd event that exhibited stress drop close to the maximum value of 6 MPa, while the 4th and 12th March events had increasingly higher stress drop values, which were similar to (or higher than) the Kozani and Umbria earthquakes (Figure 9).
Earthquake stress drop is thought to inversely relate to the notion of structural maturity, where maturity is defined as clear surface expression, large rupture length and faulting age in the order of million years (i.e., Manighetti et al., 2007). Thus, immature (or strong) faults would tend to produce earthquakes with smaller moment magnitude but higher stress drop, contrary to the more mature (or weak) faults which are usually associated with events of low stress drop. To test whether the faults involved in the TES, as well as the faults that form the relay zone, are structurally immature, we compare their displacement-length (D-L) relations with displacement-length data from a global data set of inactive faults and a compilation of global normal fault earthquakes (see gray and orange circles, respectively, in Figure 7b) (Nicol et al., 2017;Wesnousky, 2008). Examination of Figure 7b suggests that the faults in Thessaly show a clear positive D-L relationship with a slope of ∼1.1 (R 2 = 0.52) and plot in a comparable part (although mostly at the lower-end) of the log-log graph of the global data set of inactive faults, which typically have D-L slopes that range from 1 to 1.5 (e.g., Nicol et al., 2017 and references therein). The similarity in the slopes of these two fault populations, as well as the disparity between the D-L slope for Thessaly faults and the single-event normal fault global compilation (1.1 vs. 0.3), indicates that the Thessaly faults are neither structurally immature (i.e., Karakostas et al., 2021;Kassaras et al., 2022) nor newly formed (i.e., Koukouvelas et al., 2021). By contrast, the TES mainshocks (red circles with black outline in Figure 7b) appear to occupy a similar part of the graph with faults in Thessaly of demonstrable Holocene activity (magenta circles in Figure 7b), as also suggested by recognition of geomorphic indicators for past rupture on these faults (Section 3.3), and be part of a well-established large-scale relay fault zone (Figure 7a). A consequence of the above is that high stress drop values during the TES more likely reflect high-strength asperities formed at fault intersections and less so structurally immature faults.

Progressive Rupture Across a Relay Zone
The anticorrelation between stress-drop and earthquake magnitude (and fault rupture length) during the TES mainshocks suggests that, as the rupture sequence propagated to the northwest with a gradual counterclockwise rotation of its strike (Figure 8a), it broke increasingly stronger asperities formed at fault intersections ( Figure 8b). Such dependencies are supported by global field and laboratory studies (Brodsky et al., 2016;Goebel et al., 2015;Kemna et al., 2021;Walters et al., 2018), but also numerical simulations (Zielke et al., 2017), that relate high stress drop to strong fault asperities at fault intersections. Further, laboratory experiments have also shown that strong asperities on fault planes may often result in longer earthquake recurrence on these faults (Beeler et al., 2001). This is in good agreement with our geological mapping that suggests that the structures that ruptured during the TES were in their vast majority faults with subtle topographic signature which were largely unbroken for thousands of years (see "green" faults in Figure 3).
From the above we conclude that the 2021 TES propagated north-westward, breaking increasingly stronger asperities at fault intersections, transferring slip at the northern tip of a well-established relay zone that remained largely unbroken in the Holocene. Our findings, thus, provide a rare snapshot of the growth of a large-scale relay zone during sequential normal fault earthquakes and highlight the significance of normal fault intersections to earthquake rupture propagation. Fault plane heterogeneities appear to be responsible for earthquakes at the higher end (≥6 MPa) of the expected stress drop values globally. This may have significant repercussions to seismic hazard, as asperities due to geometric fault complexities have historically produced earthquakes with high stress drop and significantly enhanced ground motion amplitude that affected populated areas (i.e., the 2010-2011 Canterbury Earthquake Sequence in New Zealand, the 2016-2017 Central Italy Sequence, etc.) (Kemna et al., 2021;Oth & Kaiser, 2014;Walters et al., 2018). Thus, we conclude that fault relay zones, but also normal fault systems that include elements in close proximity to one another, are prone to grow by accumulating slip due to high stress drop earthquakes, posing locally elevated seismic hazard.

Conclusions
We report details on the coseismic breaching of a previously unrecognized large-scale relay zone in central Greece, through three successive normal fault earthquakes of moderate magnitude (M w 5.7-6.3) in March 2021. We use field and digital fault mapping, together with analysis of InSAR, GNSS and seismological data, to show that the Tyrnavos Earthquake Sequence (TES) ruptured the northern-end of a ∼100 km wide transfer structure. Indeed, geological fault displacements, aided by the displacements recorded by InSAR, show that the three TES mainshocks ruptured faults with subtle but clear geomorphic signatures that were previously unbroken in the Holocene, as opposed to the southern section of the relay zone that has accrued slip during the Holocene. InSAR-derived slip spectacularly agrees with the location of eight subtle, previously undetected active normal faults that, during the TES, accommodated coseismic and syn-seismic fault slip. Kinematic modeling coupled with fault mapping suggests that all involved faults are interconnected at depth, with their fault-intersections acting as barriers to coseismic rupture propagation. As a result, the three TES mainshocks were characterized by stress-drop values at the high-end of those recorded by Mediterranean earthquakes during the last ∼40 years (≥6 MPa). These findings, collectively suggest that the TES propagated north-westward to rupture increasingly stronger asperities and transfer slip between the tips of a well-established, but previously unrecognized, large-scale relay structure within a normal fault system. Thus, fault relay zones are prone to produce high stress drop earthquakes, posing locally elevated seismic hazard. 000GRC), for providing the GNSS data analyzed here. To generate Figures 2-4 we used QGIS (version 3.16). For the InSAR data processing and analysis we used the open-source GMTSAR software (Sandwell et al., 2011;Xu et al., 2017). For the kinematic modeling we used the modular open-source software of the Pyrocko project . Andy Nicol is thanked for useful comments. Two anonymous reviewers are thanked for their constructive comments that helped improve the content of this article.