The Geological Structure and Tectonic Complexity of Northern Thessaly That Hosted the March 2021 Seismic Crisis

Abstract

Knowing the rich presence of active faults in northern Thessaly and the lack of any significant seismic activity since at least the mid-1940s, the 2021 seismic sequence did not surprise us. What did surprise us was the fact that (i) despite the great knowledge of the neotectonic faults in the area, the causative faults were unknown, or almost unknown; (ii) the direction of the 2021 faulting was different than the expected, and given that the focal mechanisms showed almost pure normal dip-slip motion, the extensional main axis was also different than the one we thought we knew for this area; and (iii) besides the co-seismic ruptures that occurred within the Domeniko-Amouri basin and along the Titarissios River valley, there is evidence of rupturing in the alpine basement of Zarkos mountains. After thoroughly reviewing both the alpine and neotectonic structural setting and all the available literature concerning the seismotectonic data and interpretations of the 2021 sequence, including investigations of our own, we end up in a complex tectonic setting with older alpine structures now operating as inherited faults, and we also suggest the possible occurrence of a roughly N-dipping, low-angle, detachment-type fault. This fault runs below Mt Zarkos, reaching at least the Elassona Basin, with splay faults bifurcating upwards from the main fault zone. Following this complexity, rupture of the first mainshock must have chosen a split route reaching the surface through the gneiss rocks of Zarkos and almost (?) reaching the basinal sediments of the local tectonic depressions. This seismic sequence is a perfect case study to shed some light on the tectonic and rupture processes in the context of both geodynamics and seismic hazard assessment.

1. Introduction

The 2021 seismic crisis in Thessaly (central Greece) consisted of two strong shocks: the first occurred on March 3 with a magnitude of M_w6.3, and the second, one day later, on March 4, with a magnitude of M_w6.0 (Figure 1), followed by a strong (M_w5.6) aftershock on March 12. The first strong shock was not preceded by pre-shocks, but itself caused extensive damage to the buildings in the surrounding area. The second shock worsened the already suffering buildings, but fortunately no fatalities were reported by either the first or the second strong shock. The secondary ground deformation phenomena were plenty: ground fissures and liquefaction, often observed together, were mainly constrained in the basinal part of the epicentral area, whereas rockfalls were scattered mostly on the mountainous areas.

The geological setting of the broader epicentral area (northern Thessaly) has been extensively studied by many researchers, especially during the 1990s and 2000s in the context of alpine [1,2,3,4,5,6,7,8] as well as active tectonics [9,10,11]. Many prevailing tectonically active structures have been identified, although no significant earthquakes have been recently instrumentally recorded, as will be shown afterwards. The March 2021 was the first significant event to be well-captured by the Hellenic Unified Seismic Network. Surprisingly enough, both events did not follow what was already known in the past. The focal mechanisms proposed a different orientation of faulting, and hence extension, since the movement was defined as quasi-pure normal dip-slip. What was also unexpected is the fact that, despite the secondary and possibly primary ruptures found along the southern margins of the Domeniko-Amouri basin and the Titarissios River valley, we also found co-seismic ruptures and alpine faults with geometry and kinematics almost perfectly matching the 2021 focal mechanisms, in the mountainous area of the alpine basement further to the south. Thus, the question that rises is: Do the older, ‘inherited’ structures still define the orientation of faulting or do we have a wrong view on the modern geodynamic setting of this part of the Aegean region? To attempt answering this question, we first need to go deeper to the review of the alpine setting.

2. Geological Setting

2.1. Alpine Structural Evolution

The alpine structural evolution of the Hellenides began during the Early-Middle Mesozoic. Discontinuous westward and eastward orogenic migration, characterized by intense collisional tectonics, related to successive subduction processes of the Tethyan oceanic basins, produced a complete stack of several nappes of Jurassic to Tertiary age. Extensional tectonic took place in between the nappes stacking processes causing tectonic denudation and crustal exhumation [16,17,18,19,20,21,22,23,24,25,26,27,28,29]. Various palaeogeographic models have been proposed arguing about the tectonic evolution of Tethys and the emplacement of the ophiolites, most of which are extensively discussed by Robertson et al. [30]. Towards the final stages of the alpine cycle, the region is marked by the convergence of the Eurasian and Nubian (African) plates, passing from the Mesozoic to the Tertiary, with the consequent closure of the intermediate oceans (Tethys Ocean) and the consequential back-arc extensional regime that led to the orogen collapse through large, normal, low-angle detachment faults often causing the exhumation of underlying units as metamorphic core complexes [3,27,31,32,33,34,35,36].

The main feature of the Middle-West Greek orogen is the successive tectonic nappes and slices, which are thrust on top of each other from east to west. The thrusting events started in the Middle-Late Jurassic with the obduction of the ophiolites onto the Pelagonian carbonate platform from one or more Tethys oceans [23,26,27,28,37]. Paleogeographically, the Pelagonian zone is considered to be a large continental shelf, part of the Cimmerian continent, which was detached from Gondwana (i.e., Africa with other smaller continents; [26,38]). During the Tertiary, the tectonic evolution continues with the thrust orientation toward the west-southwest and the thickening of the continental crust. The accumulation of successive nappes and slices, in combination with the development of subduction zones, is associated with the creation of two sub-parallel metamorphic zones of high pressure/low temperature (HP/LT) of Paleocene-Eocene and Oligocene-Miocene age. The extensional regime that followed the compressive tectonics led to the syn- and post-orogenic collapse with large, low angle, normal detachment faults, and the exhumation of the underlying zones in the form of tectonic windows in Olympos-Ossa, Cyclades, and Peloponnese-Crete areas (Figure 2 and Figure 3; [27,31,32,33,34,35,39,40]). Extensional tectonics has played a crucial role from the Tertiary to the present day with extensive shear zones and detachment faults. The reactivation of the older reverse faults (inherited structures) as extensional structures (inverse tectonics) set a study case to be revised in correlation with the recent earthquakes.

2.2. Alpine Tectonostratigraphy

Our study area (Figure 1) lies on the Pelagonian zone (or nappe) which consists of a Palaeozoic or older crystalline basement (gneiss, gneiss-schists and amphibolites) intruded by Upper Palaeozoic granitoids (Carboniferous age). It is overlain by a low metamorphic Permo-Triassic volcano-sedimentary series, composed by intercalations of marble layers, phyllites, sandstones and bimodal volcanic products. Triassic-Jurassic crystallised limestones and dolomites follow at the top of the Pelagonian stratigraphic column. Characteristic is the occurrence of the tectonically emplaced, initially during the Upper Jurassic, Neo-Tethyan ophiolite masses on the Pelagonian Triassic-Jurassic carbonate platform and secondary imbricated during the Cretaceous and Tertiary with the Pelagonian formations. Between the underlying Olympos-Ossa Triassic to Eocene-Oligocene platform carbonate formation and the continental block of the Pelagonian, High-pressure rocks (Blueschists of the Ampelakia Unit) are tectonically emplaced. They were metamorphosed in high pressure conditions during the Paleocene-Eocene, when they subducted under the Pelagonian block. Thus, after the Eocene, the Pelagonian formations were thrust upon the carbonate series probably of the External Hellenides, which were revealed afterwards as a tectonic window (Figure 2; [26,33,41,42,43]).

The study area is dominated by the Palaeozoic basement rocks and their overlain Triassic-Jurassic neritic carbonate nappe, which are tectonically emplaced on a continuous carbonate series (crystalline limestones and marbles) of Triassic-Jurassic and possibly Cretaceous age, also known as the carbonate series of Kranea, which is considered relatively autochthonous exhumed as the tectonic window of Kranea, equivalent to the Olympos-Ossa tectonic window more eastern of the study area (Figure 1 and Figure 2; [33,41,44]).

2.3. Alpine Tectonics

The study area is characterised by a complex alpine tectonic deformation as a result of many tectonic events that have affected the described geological formations of the area [32,45,46]. A complete view is given below (Figure 4).

A predominant S1 schistosity is recognized in the broader area. It is SW-ward dip-directed at the west-southwestern part of the Kranea tectonic window (Study area) and NE-ward dip-directed at the north-northeastern part of the window. On the S1 planes, a L1 lineation is recognized in a consistent NE-SW to NNE-SSW orientation but NE-ward plunging at the north-northeastern part of the Kranea window and SW-ward plunging at its south-southwestern part. Both S1 and L1 are due to a deformation event (D1) associated with the SW-ward thrust of the Pelagonian nappe upon the Kranea carbonate series, as it is shown by shear criteria (e.g., S-C fabric, asymmetric boudins, mica fish, σ and δ clasts etc.). Isoclinal folds and sheath-folds also developed during D1, sup-parallel to parallel to the L1 stretching lineation after rotation and hence parallel to direction of the maximum extension (X-axis). During the last stages of the D1 deformation, extensional SB mylonitic shear zones were created and in the upper tectonic horizons low-angle normal faults, both dipping to the SW and NE. The presence of glaucophane (HP/LT conditions) on the SB planes related to a DSB deformational stage indicates that it initially developed at great depth, as the glaucophane is preserved rotated along the SB or S1 planes and parallel to the L1 lineation. The evolution of the D1 deformation in the extensional DSB stage resulted in the creation of extensive shear zones while the orogen emerges with the simultaneous tectonic escape of the accumulated nappes. The D1 is associated with a greenschist metamorphic facies (pressure up to 6–7 Kb and temperature around 400–480 °C). On the other hand, the HP/LT metamorphism and the development of the glaucophane characterize the period of the subduction of the geological formations below the Pelagonian continental block during the Paleocene-Eocene. An older, residual synmetamorphic deformation is also detected on the Pelagonian series, but strongly reworked by the already described Tertiary deformational events (compression vs. extension). It is possibly of Late Jurassic-Early Cretaceous age related to the first Alpine structural evolutionary stages of the Hellenides. Older Palaeozoic structures cannot be recognized in the studied area because they are totally disappeared by the younger alpine tectonics.

A compressive younger D2 deformational event followed, possibly during the Early Miocene, at very low P/T conditions, characterized by knick-folds or zones, shear zones and thrust faults. At this stage, the orogen has already emerged. The thrusting of the Pelagonian nappe and the HP/LT rocks of Paleocene-Eocene age on the underlain Kranea carbonate formation is assumed to have taken place during the Eocene-Oligocene. From the Late Oligocene to the Early Miocene, the orogen continued its emergence, while at the same time it collapsed (extensional event DSB). The following Neogene-Quaternary extensional tectonics, which is directly related to the alpine one, as will be discussed below, is described in the next section.

3. Post-Alpine Geology

The large Neogene basins in the central and northern Greek mainland are formed by the post-alpine extension of the Late Miocene that followed the final ENE-WSW-trending alpine compressional stage in the Early-Middle Miocene (Aquitanian-Langhian; e.g., [9]), which is still active further to the west, in Epirus and the Ionian Sea [9,47,48,49]. The post-orogenic collapse phase, following the convergence front of the Eurasian and African plates, led to the first NE-SW oriented extension. Given that this phase is responsible for the initial formation of most of the large basins, the major trend of the basins has a NW-SE direction, shaped by NW-SE-striking normal faults.

In our study area, three major basins are met, the intermontane Elassona-Domeniko-Amouri basin, and the much larger Eastern (also known as the Larissa plain) and Western (also known as the Karditsa plain) Thessalian basins (central mainland Greece) (Figure 1). The former is the northwestern one and can be subdivided into two smaller sub-basins, the Elassona to the NE and the Domeniko-Amouri to the SW. These two sub-basins are connected through two narrow valleys (Figure 1). The Larissa plain is located to the east and is separated from the previous basin by a mountainous, marble-consisting mass; there is only a narrow passage that connects the Larissa plain with the Domeniko-Amouri sub-basin, i.e., the Titarissios River, which flows toward the Larissa plain. The West Thessalian basin is located south of the Domeniko-Amouri sub-basin and west of the Larissa plain. It is completely isolated from the former by the Zarkos mountains and connected with the latter via the Penios River, flowing toward the Larissa plain (Figure 1). All basins are separated by the alpine rocks of the Pelagonian zone, which are either schists and/or marbles.

Based on the oldest basinal deposits, all the aforementioned basins started to form in the Late Miocene, right after the last Early-Middle Miocene alpine compressional stress field.

3.1. Post-Alpine Lithostratigraphy

The neotectonic activity, along with the local palaeoclimatic and palaeogeographic conditions, facilitated the formation and shaping of the Elassona-Domeniko-Amouri basin and its fill with the Neogene and Quaternary sediments. The Late Neogene (Pliocene) sediments are unconformably deposited upon the schist-crystalline basement of the Pelagonian zone, they are significantly developed, and they mostly outcrop along the margins of Domeniko-Amouri sub-basin. During a mine project seeking the lignite deposits in the Domeniko-Amouri sub-basin, a dense network of boreholes was drilled, revealing important tectono-stratigraphic information for our study area (Figure 5; [50]). The stratigraphic description of the findings includes (from older to younger; Figure 5):

• Basal formation of the Upper Miocene unconformably lying over the alpine basement, with a maximum thickness of 55 m. It consists of mud, sandy and silty clay, with granule and gravel layers, representing alluvial fan deposits that filled the initial stage of the basin’s formation. mostly of fluvio-torrential origin.
• The gray-green lignite-bearing formation of the Upper Miocene, consisting of lignite beds, mostly at the upper members, within clayey silts, sandy clays, sands, clays, and rarely silts. Fine-grained sand passes into medium- and coarse-grained toward the lower members. Conglomerate intercalations also occur, rarely solid. Fossils lack. The lithofacies represent a fluvial environment developed by a meandering river. Maximum thickness of 150 m.
• The overlying clastic deposits separated into (i) a 50 m-thick, fine-clastic lower member of Upper Miocene, consisting of friable and rarely stiff siltstones, sandy clays, sands, with local intercalations of conglomerates, and (ii) a 30 m-thick, coarse-clastic, unconformably overlying upper member of Upper Pliocene, consisting of unconsolidated to loose breccio-conglomerate with sand and boulders.
• The Pleistocene deposits divided into two members: (i) the Upper Pliocene? to Lower Pleistocene (Villafranchian), 70 m-thick, terrestrial deposits composed by khaki sandy clays with mud, and hard, solid, white limestone intercalations, and (ii) the Upper Pleistocene, 60 m-thick, terraces and talus cone deposits, consisting of brown mud with sandy gravels and conglomerates, occasionally in thin interbeds.
• The Holocene eluvial deposits and river terraces. They can reach 25 m of thickness.

According to the IGME’s [12,13] geological maps, the oldest post-alpine sediments within Eastern Thessalian basin are deposited unconformably upon the alpine basement in the Upper Miocene; they comprise of transgressive, polygenetic, compact conglomerates, passing upwards to thick-bedded, micro-brecciated limestone. Lacustrine deposits of marls and marly sandstones follow in unconformity, and above them fluvio-terrestrial formations of sandy-clay material and loam, both of Plio-Pleistocene age. Pleistocene fluvio-lacustrine deposits of clays and sands with variously thick intercalations of coarse-grained material occupy a large extent of the northern margin, and on the top, Holocene alluvial deposits have covered the greatest part of the plain, consisting of unconsolidated clays, sand, angular and rounded pebbles, and fluvio-torrential-lacustrine material.

The Western Thessalian basin is partially imprinted on the older molassic formations of the Meso-Hellenic Trough, a (piggy-back or pull-apart) basin developed during the last alpine stages, from the Mid-Late Eocene to the Mid-Late Miocene [6,7,8,52,53,54,55]. As in the other basins, the older post-alpine sedimentation started in the Late Miocene with the transgressive, polygenetic, compact conglomerates, which was followed by terrestrial clastic deposits of the Pliocene, fluvio-lacustrine deposits of the Pleistocene, and the alluvial deposits of the Holocene [14,15,18].

River terraces along the local fluvial systems, and slope debris with talus cones along the margins, also occur in all basins.

3.2. Post-Alpine and Active Tectonical Setting

The post-orogenic collapse that followed the Eurasian-Nubian plates convergence led to the first NE-SW-trending extensional (σ₃) phase which lasted from Late Miocene until Pliocene (P2, Figure 6 and Figure 7) [9]. This phase was responsible for the formation of the large basins in Thessaly following the main alpine structures. The next extensional phase (P3, Figure 6 and Figure 7) [9] has been active since Middle Pleistocene and remains active until Today. It has a general N-S direction (roughly N10°E of the σ₃ axis) and affects most of the broader Aegean region (Figure 7, bottom-right) [9]. The stress field is reflected in the crustal deformation of the area which is depicted by GPS/GNSS analyses. Indeed, the upper crust seems to be dominated by dilatation, but towards the western margin of the Western Thessalian basin, i.e., toward the Pindos mountain range, shear locally occurs with the presence of contraction (Figure 8) [56].

The result of the active stress field is the creation of a new generation of normal faults preferably perpendicular to the extension (i.e., E-W direction, N100° to N110° systematically), which overprint their geomorphological imprint along the margins of the Thessalian basin (Figure 6 and Figure 7d). The faults of northern Thessaly have been thoroughly studied, especially since the 1990s, resulting in the recognition of many active faults. The ones that affect our area of interest are the Tyrnavos and Larissa faults.

The Tyrnavos fault is a typical normal E-W- to WNW-ESE-striking, N-dipping structure considered as active, similar to the many others in the broader area [9,11,61]. It is clearly seen cutting the Mesozoic marble and it crosses diagonally the Titarissios valley.

It also shows typical morphological characteristics affecting younger sediments. It enters the Larissa plain from the west as proven from geophysical and palaeoseismological investigations near the basinal margin [61,62]. Its surficial length is estimated to 10–12 km, although its basinal length remains unknown. Further details for this fault can be found in the GreDaSS [63,64].

The Larissa fault is a sub-parallel (striking more to the WNW-ESE) normal fault that crosses most of the Larissa plain with an inferred length of ca. 20 km reaching the alpine basement to the west [62,63,65]. Once again, further details for this fault can be found in the GreDaSS [63,64]. However, based on the 1:50,000 scale geological map of IGME [13,15], the alpine basement is fractured by a fault system of the same direction, parallel to the passage that hosts the Penios River, which seems to reach the marginal fault of Titarissios valley (yellow arrows in Figure 1a), implying a possible connection.

The earthquake activity in northern Thessaly lacked any significant earthquakes before the 2021 sequence. Strong to major historically and instrumentally recorded events [57,58,60] mostly occurred in the Eastern and Western Thessalian basins with the last one occurring in 1941 (Figure 7b). The very recent activity (since 2011, when the completeness of the earthquake catalogue [58] was more improved) shows a large gap in Eastern Thessalian Basin, although large tectonically active faults occur (Figure 7d), and a small, clustered, diachronic activity in the epicentral area (Figure 7c). The Western Thessalian Basin shows much more frequent activity during those last 10 years.

The Domeniko-Amouri Basin Tectonic Setting

The previously mentioned mine project for the lignite deposits in the Domeniko-Amouri sub-basin with the borehole network [50] revealed tectonic structures that are buried under the Holocene alluvial deposits. After correlating the borehole logs, vertically displaced (downthrown) horizons and the basin’s base were observed indicating the presence of normal Faults (Figure 5b,d). Given that the layer succession is not always cut up to the top, the detected faults are interpreted of having variable ages: sometimes they displace only the deeper formations of the Late Miocene, characterized as “Neogene” faults, whereas other times they reach and displace the base of the Quaternary (Villafranchian), characterized as “post-Neogene” faults [50,51]. The latter, which are of our interest, have mainly a NW-SE direction forming bookshelf or graben-style patterns (Figure 5b,d). Along the southwestern margin of Domeniko-Amouri sub-basin, a fault of this direction, dipping to the NE, demonstrates a downthrow of several tens of metres (fault f11; Figure 5). As we will see below, this fault lies along the coseismic ground ruptures with the liquefaction phenomena, mapped southwest of the Titarissios River and is possibly related with the second strongest shock (March 4). The Titarissios fault (yellow arrows in Figure 1a), which lies SE-wards along the southern margin of the Titarissios valley, had been previously mapped [15,65], without knowing how accurately and by which method, but never elaborately studied. This fault forms a prominent linear slope bringing to contact the alpine basement (schist-gneiss-marbles) with the valley’s alluvial deposits (Figure 1a).

4. The 2021 Seismic Sequence

In March 2021, two strong earthquakes were recorded close to the towns of Tyrnavos and Elassona (Figure 7d and Figure 9): the first one occurred on March 3 with a magnitude of M_w6.3; the second one occurred one day later, on March 4, few kilometres WNW of the first one, with a magnitude of M_w6.0. Both events and the rich aftershock activity that followed them were shallow, with depths rarely exceeding 12–15 km [51,67,68] (Figure 9). According to Michas et al. [69], “the aftershock sequence exhibited scaling properties consistent with well-known empirical relationships for aftershocks. In particular, the frequency of aftershock magnitudes follows the Gutenberg-Richter scaling law for a b-value of 0.90 ± 0.03”. Published focal mechanisms from various institutes revealed NW-SE-striking normal faulting, and hence, NE-SW oriented extension (σ₃), which is quite different to the ca. N-S estimated from neotectonic investigations in previous studies. This faulting direction is also delineated by the horizontal distribution of the aftershock sequence [51,67,68].

The two strong shocks with the few kilometres distance between them strongly suggest the occurrence and reactivation of two adjacent and roughly aligned fault segments. This scenario can be also supported by the spatiotemporal evolution of the sequence (Figure 9). Indeed, the aftershocks that followed the first main event, and exactly until the second main event (Figure 9a), were constrained to the west as far as the second mainshock lies. Right after the second main event, most of the aftershock activity migrates toward the WNW on a parallel direction, delineating the neighbouring sub-parallel segment. The overall picture shows a WNW-ESE- to NW-SE-trending horizontal distribution of the hypocentres extending from the western margin of the Larissa plain to the east, until the eastern slope of the Antichasia Mountains (toward the Domeniko-Amouri basin) to the West.

4.1. Ground Deformation Phenomena

The earthquake sequence produced a wide range of ground deformation phenomena, including (primary and secondary) coseismic ruptures, liquefaction, rockfalls, and landslides.

Coseismic ground ruptures abounded in the Domeniko-Amouri basin and the Titarissios valley. Almost persistently, they remained in a NW-SE direction, following the shape of the valley, parallel to the strike of the main faulting. They usually crossed farmed fields and fluvial deposits near the Titarissios riverbanks, as well as manmade constructions, such as asphalt roads and local irrigation network channels made of reenforced concrete. In fact, whatever the constructions direction was, the ruptures continued almost unhindered. However, not all of these ruptures are considered primary. Taking into account the alluvial-covered ‘Post-Neogene’ faults found by Dimitriou and Giakoupis [50], we suggest that the ruptures above these faults are probably primary. The ruptures were often observed to be combined with liquefaction. Ruptures with liquefaction were observed further SE, in the Titarissios valley. These ruptures are very surficial and considered gravitational and they were probably created by the sympathetic slip of the Titarissios fault as it will be explained in Discussion section.

As we will later see in the InSAR results, unexpected was the fact that (near-to-the) surface rupturing is indicated further south-west of the Titarissios valley, within the gneiss-schist alpine basement, where few coseismic ruptures were found in a thin soil cover (Figure 10f). Moreover, considering the focal mechanisms, as well as the aftershock distribution hypocenters, we searched for matching faults in this particular area. Indeed, slickenlines of normal dip-slip movement were observed on N160°E-striking, moderately dipping (~50° average), polished surfaces (Figure 10d,e), as well as cataclasite and micro-fractures, implying the neotectonic reactivation of a pre-existing shear zone in brittle conditions. These slickensides are directly comparable to the proposed geometry of the focal mechanisms. Moreover, these ruptures partially coincide with the broader fault pattern of ‘unknown’ activity that was recognised in IGME’s [13,14] geological maps (Figure 1a), which implies the tectonic connection of the Titarissios valley with the Larissa plain through discontinuous NW-SE-striking faults.

Extensive liquefaction phenomena were manifested not only in the Domeniko-Amouri basin and the Titarissios valley (Figure 11 and Figure 12), but in Penios valley (Peniada and Zarko) as well [51,70]. They were mostly continuous and aligned along ground fissures.

Secondary gravitational effects other than liquefaction, such as rockfalls and landslides, were observed throughout the area (Figure 11 and Figure 12). Their distribution did not follow a specific pattern, but was rather dependant on the local site conditions, i.e., rock formation attributes and slope angle. Quite common was lateral spreading observed next to riverbanks, such as in Damasi, Mesochori, and other nearby areas.

The extent and intensity of ground deformation are revealed from InSAR (Sentinel-1) images (Figure 13) [60]. Utilizing a pair of images just before and after the first mainshock, the uplift and subsidence were estimated, also delineating the surficial location of the seismic fault in the alpine basement of the Zarko Mt. Similar results have been published by other research teams as well [67,71,72,73,74,75,76,77]. The results of Yang et al. [77], based on InSAR analysis and relocated aftershocks, suggest that at least four unmapped low-angle normal faults were activated. They also suggest that afterslip propagated along a shallow steep fault while coseismic slip occurred on a deep, moderately dipping fault, revealing a ramp-flat structure.

4.2. The Seismic Sources

The determination of the seismic sources that produced the two mainshocks has been the main issue in this work, as well as in many already published papers. The principal evidence on which all research teams based their models are (i) the hypocentral locations, especially of the two mainshocks, (ii) the co-seismic ground ruptures, and (iii) the InSAR images [51,60,69,71,76,78]. Some researchers focused only on the Domeniko-Amouri basin and the Titarissios valley where most of the ground ruptures occurred [51,72]. However, few paid attention to the fresh slickensides southwards, within the alpine gneiss-schist basement [60,71]. Taking into account both sets of primary coseismic ruptures, the InSAR fringes, and the geometry suggested by all published focal mechanisms, we suggest that two roughly parallel faults were ruptured during the first strong event: the southern one, the Zarkos fault (orange arrows in Figure 1a and Figure 14), accommodated the main slip, whereas the northern one, the Titarissios fault (yellow arrows in Figure 1a and Figure 14), was partially ruptured generating the eastmost ground fissures along the Titarissios valley. The rupture of the northwestern adjacent fault segment, i.e., the one that produced the next day’s (March 4) strong shock probably produced the westmost ground fissures between the villages of Vlachogianni and Amouri and is probably related to the blind ‘post-Neogene’ faults of Dimitriou and Giakoupis [50] (Figure 11). The exact location and length of the seismic fault, however, were not possible to obtain and for this reason we did not present this information on a map.

As shown in a previous work [60], we suggest that the March 3, M_w6.3 seismic source is best represented by a 310°-striking, 44°-dipping, with −89° rake fault model, as the focal mechanism of the GFZ GEOFON program suggests (Figure 14). Given that the basement’s primary ground ruptures are in good agreement with the InSAR fringes, we matched the seismic source with the Zarkos fault. Based on the dimensions calculated by magnitude vs. rupture dimensions empirical relationships [79] and the mainshock’s hypocentral location, we estimated a length of 19 km and a down-dip width of 11 km. The spatial aftershock distribution of Figure 9 and the hypocentral relocation of the mainshock [67] were also taken seriously into account for the geometry and location estimation.

Based on our fault model, we calculated the theoretical ground vertical displacement of the first mainshock (March 3) using the Okada formulae [80], which demonstrates a maximum subsidence of 29.3 cm (Figure 15; [60]). The interpolated pattern shows great similarity with the InSAR findings with a comparable maximum value (35 cm from the InSAR).

There is no doubt that the two mainshocks were produced by two adjacent and aligned faults or fault segments, with a small strike difference. The ground ruptures in the Domeniko-Amouri basin and along the Titarissios valley strongly support the surficial expression of the seismic faults, but they do not coincide to InSAR results. If we also take into account the hypocentral location and the proposed focal mechanisms of the first event, and the upward theoretical fault plane continuation, it becomes possible the occurrence of a low-angle, detachment-type normal fault with upward strands, similar to the dominant theory of the Corinth Gulf tectonic setting [81,82,83,84]. This means that the rupture of the first major event bifurcated along two strands, (the northern) one reaching the Titarissios valley, and the other (southern one) cropping out in the alpine basement on Zarkos Mountains. A corresponding situation failed to be observed for the second major event. The detachment hypothesis is also referred by Koukouvelas et al. [72], noting the complex geometry caused by the possible occurrence of a low-angle, detachment fault.

The study of the static stress change patterns due to fault slip is crucial in order to understand the effect of the seismic source on the surrounding receiver-faults; this study is achieved with the Coulomb static stress modelling. Plenty of the published works about the currently studied earthquake sequence present stress transfer models which diversify due to the input parameters used in each case (e.g., fault location and geometry) [60,67,68,69,74,76,77]. Chatzipetros et al. [60] calculated the Coulomb static stress changes at a depth of 8 km for receiver faults similar with the seismic source (Figure 16). Their source fault model was based on the published focal mechanisms and the hypocentral spatial distribution of the first mainshock (March 3). Their results show bilateral stress load along fault strike which can explain the triggering of the adjacent fault (NW-ern ‘segment’ of the Zarkos fault zone) that produced the second mainshock on March 4. In the study of Kassaras et al. [68], several stress models clearly show the different stress patterns of each one of the main events corresponding to the relocated aftershocks cloud.

5. Discussion

The March 2021 earthquake sequence with its two strong mainshocks rose new queries concerning what we knew about the previously studied active faults and the extensional stress field they belong to. Given the roughly N-S extensional stress field, we thought it was the prevailing direction in the epicentral area, which turned out to be NE-SW according to the focal mechanisms and the observed ground ruptures. This direction is not unfamiliar in the broader area as it used to be the dominant direction during the previous tectonic phase (P2, Figure 7; [9]). Moreover, it is close to the late alpine, compressional P1 tectonic phase (Figure 7; [9]) that caused the NNW-SSE-striking thrusts in the region. Then, the following questions arise: Do the previous (inherited) tectonic structures still act as active in this area? Can formerly thrust alpine faults be reactivated as normal during the afterward extension?

If we also consider (i) that the epicentral area lies next to the spine of mainland Greece, i.e., the NNW-SSE-trending Pindus mountain chain, a region where instrumentally recorded seismicity is rather low, and (ii) that further to the west (Epirus coast and northern Ionian Sea) the stress field has already turned into compressional, then perhaps this σ₃ axis rotation is part of the gradual transition between the N-S extension to the east and the ca. E-W compression to the west. Kassaras et al. [68] showed a rotation of the main extensional axis along a W-E direction in the broader region of central Greece, i.e., from NW-SE direction in Pindos to roughly NE-SW in the 2021 epicentral area. The compressional axis significantly diminishes along the same direction, almost disappearing much western than the epicentral area. Consequently, this question arises: does the poorly formed Elassona-Domeniko-Amouri basin (compared to the Eastern and Western Thessalian basins) represent the westward migration of extension like Koukouvelas et al. [72] also suggest?

Besides our suggested fault model for the first strongest shock [60], various seismic sources have been modelled by many researchers for the strongest shocks of the Thessalian seismic sequence; these are briefly described below.

The geodetic inversions of De Novellis et al. [74] 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 M_w6.3 (March 3) and M_w6.0 (March 4) events, and one belonging to the SW-dipping normal fault associated with the M_w5.6 (March 12) event. The respective fault models of Ganas et al. [71] are of similar geometry: blind NE-dipping fault planes for the first two mainshocks and a blind antithetic, along a straight line with the others, to the northwesternmost part of the sequence (almost perpendicular to the E-dipping slope of Antichasia Mountains). This fault setting was also adopted by Kassaras et al. [68]. Kontoes et al. [75], based on the InSAR results and the proposed nodal planes, modelled all three previous events dipping to the NE, with the fault plane lying a couple of kilometres (1.4–3.1 km) below the surface. Koukouvelas et al. [72] insist that there is an apparent mismatch between field and interferometric data implying that the responsible fault for the first mainshock (March 3) is the SSW-dipping “Mesochori” fault (Figure 17), i.e., a fault that bounds the northern flank of the Titarissios River. Given that our suggestion for the first mainshock’s fault is in accordance with the majority of the published studies about this seismic sequence, we would interpret this fault as an accommodation structure of the Titarissios fault or a sympathetic fault (Figure 17). However, the above authors do not directly connect the three strongest shocks with their causative faults, only suggesting that at least three fault segments were reactivated by the major shocks of the March seismic sequence following the Titarissios river graben as NW extension of the Tyrnavos fault.

Song et al. [76], using InSAR techniques to model the fault plane, showed that all three successively occurring events (March 3, 4, and 12) were characterised by normal dip-slip kinematics on previously unknown faults confined in the upper 15 km of the crust. In particular, the March 3, M_w6.3 earthquake ruptured a NE-dipping fault (301° strike and 46° dip), producing a maximum slip of ca. 2.2 m. The slip motion of the March 4, M_w5.9 aftershock showed a similar fault geometry (297° strike and 42° dip), but with a considerably smaller dip-slip component (~0.8 m).

The seismological and geodetical investigation (hypocentral relocation, teleseismic P-waveforms inversion, moment tensor calculation and InSAR) of Papadopoulos et al. [73] result in three seismic sources corresponding to the three strongest shocks of March 3, 4, and 12 respectively: the first shock (M_w6.3) was produced from a fault segment striking 314° and dipping NE41°, causing surface subsidence of −40 cm and seismic slip of 1.2–1.5 m at the depth of ~10 km; this model is in good agreement to our results. In contrast with other studies, the second earthquake (M_w6.2) occurred to the NW, on an antithetic subparallel fault segment (strike 123°, dip SW44°) producing seismic slip of 1.2 m at the depth of ~7 km, and surface subsidence of −10 cm. Their third source for the March 12 (M_w5.7) shock is considered as the north-western continuation of the previous fault (strike 112°, dip SSW42°) which caused ground subsidence −5 cm and seismic slip of 1.0 m at the depth of ~10 km.

In terms of seismic hazard assessment (SHA), the 2021 earthquake sequence showed us that more research is still needed to better determine the seismic sources in the broader area of Thessaly. The seismic gap proposed by Caputo [11] was a fact indeed, putting us in wait for the next earthquake to occur. However, in the 2021 sequence, faults that were either unknown, such as the Zarkos blind fault, or faults which were never expected to be active (e.g., the Titarissios fault marked by the yellow arrows in Figure 1a), came to our knowledge, no matter how well this area was studied before. The Titarissios fault could also be considered as the prolongation of the Larissa fault, although no sufficient evidence exists. The complex geometries and rupture process add new considerations in SHA. The Coulomb stress change estimations indicate a triggering effect between the two mainshocks, as would have been expected for tangentially aligned normal faults (stress rises bilaterally of the fault plane’s tips).

6. Conclusions and Open Questions

The experience of the March 2021 Tyrnavos–Elassona seismic sequence presented new problems that require rational answers, such as: When and how do old faults, non-preferably oriented to the modern stress field, activate? How does rupture propagate, both bilaterally (horizontally), through different segments, and upwards (vertically) following old, inherited structures? How is the morphology affected, i.e., when do normal faults occur along the margins of obvious geomorphological depressions (e.g., basins, valleys, etc.), and when do they occur in unexpected locations (e.g., mountainous areas)? Which are the best stand-alone or combined methods for detecting and recognising faults that are hidden, either by natural processes or human interventions? Any answer to any of the above questions can crucially contribute to a twofold hazard estimation: the ground motion simulation as deduced from deterministic seismic hazard assessment (DSHA), and the surface faulting hazard or surface fault rupture hazard, a crucial assessment for building and infrastructure design considering that a possible fault displacement could damage the foundations of any technical construction [78].

The most important conclusion that is drawn from the 2021 earthquake sequence is that nature can always surprise us. In this review paper, the complex alpine tectonics of the area is sufficiently described, which forms the background (upper crust) of the seismogenic volume and probably contributes to the active tectonics. The neotectonic regimes well as the surface morphotectonic structures are presented. Although typical active faults have been very well studied in the area, the seismic structures of March 2021 manifest as blind unknown faults in the shale bedrock of the Pelagonian geological zone.

• The focal mechanisms of the two mainshocks revealed NW-SE-striking fault rather than WNW-ESE. Both show almost pure normal dip-slip motion. Consequently, a NE-SW oriented extension is indicated instead of N-S.
• The faulting and σ₃ directions of the two 2021 mainshocks better match the corresponding directions of the previous P2 extensional phase.
• The Domeniko-Amouri basin and the Titarissios valley were full of ground ruptures and liquefaction phenomena, and while where there was a relief, landslides and rockfalls also occurred. Most of the time, ground ruptures followed the general trend suggested by the focal mechanisms.
• The InSAR images, along with the observation of neotectonically formed slickensides on old shear zones of the gneiss and schist comprising alpine basement, suggest a main rupturing path during the first major event with a possible secondary and steeper one aiming at the Titarissios valley (through the Titarissios fault). If we also add the hypocentral location and the nodal plane geometry, a bifurcated rupture can be inferred, possibly by a low-angle, detachment-type normal fault and its upward branch.
• Our inferred tectonic setting includes the occurrence of a low-angle, detachment-type fault dipping to the north, with the steeper parallel to sub-parallel normal faults in its hanging-wall corresponding to branches bifurcating from the detachment’s surface. The Zarkos blind fault can be considered as the almost emerged seismogenic source, while two of its NNW-SSE-trending and ENE-dipping segments were activated during the double event of March 3 and 4.
• Special attention needs to be given (i) tp the role of inherited structures in seismogenesis deviating from standard rules and known criteria, especially blind faults in mountains without any morphotectonic feature, and (ii) to the unknown, hidden, active faults and their role in SHA, especially close or under the modern urban areas and along lifelines. New methodologies and scientific tools are needed to identify the weak zones and the associated earthquake patterns.