The 27 September 2021 Earthquake in Central Crete (Greece)—Detailed Analysis of the Earthquake Sequence and Indications for Contemporary Arc-Parallel Extension to the Hellenic Arc

Featured Application: Field validation of combined remote sensing and seismological data after a large earthquake. Abstract: The Arkalochori village in central Crete was hit by a large earthquake ( M w = 6.0) on 27 September 2021, causing casualties, injuries, and severe damage to the infrastructure. Due to the absence of apparent surface rupture and the initial focal mechanism solution of the seismic event, we initiated complementary, multi-disciplinary research by combining seismological and remote sensing data processing, followed by extensive ﬁeld validation. Detailed geological mapping, fault surface measuring accompanied with tectonic analysis, fault photorealistic model creation by unmanned aerial system data processing, post-seismic surface deformation analysis by DInSAR image interpretation coupled with accurately relocated epicenters recorded by locally established seismographs have been carried out. The combination of the results obtained from these techniques led to the determination of the contemporary tectonic stress regime that caused the earthquake in central Crete, which was found compatible with extensional processes parallel to the Hellenic arc.


Introduction
The broader area of central Crete represents a neotectonic structure in the vicinity of the Hellenic Trench, which comprises two post-orogenic basins with trending orientations normal to each other, forming a uniform basin complex [1]. The northernmost part of the latter, the Heraklion Basin (HB), trends approximately N-S, whereas adjacent to its southern margin, the Messara Basin (MB) developed trending E-W [2,3]. We refer to this area as the Heraklion-Messara Basin complex, as it is covered by the same Miocene formations, implying that it shares a common paleo-environmental history, even though there are significant structural differences between the two basins.  [11]). The NNE-pointing vector at the lower-left indicates the direction and convergence rate of the Nubian plate [12]. (b) Seismotectonic map of Crete (red rectangle in panel (a). Black lines are active faults from the NOAfaults v3.0 database [13]. Beachballs indicate focal mechanisms of earthquakes (M w ≥ 4.0) from the compilation of [14] and sources referenced therein (see also Data Availability). Circles (M < 5.0) and stars (M ≥ 5.0) represent past seismicity of years 1000-1899 (SHEEC; [15]) and 1900-2009 [16] for events with focal depths down to 70 km. The maximum (S 1 ) and minimum (S 3 ) principal stress axes of the regional crustal stress field [14] are depicted by blue and red arrows, respectively, with arrow length proportional to the cosine of the plunge angle. The epicenter of the Arkalochori earthquake is presented by a star within MB and HB areas.
MB has been considered an E-W-trending, typical graben structure [17]; however, it is characterized by the absence of classic continuous bounding faults at the north and south margins of the basin. Its northern margin is probably covered with Miocene and/or Pliocene sediments. Therefore, the exact boundaries between MB and HB can hardly be observed and, consequently, several studies, with different arguments [18][19][20], implicate arc-parallel extension dynamics, trending along the local Hellenic trench axis (Figure 1).
Intense historical and recent seismicity reveal the continuous contemporary tectonic activity of the entire Cretan territory (e.g., [21]). Destructive earthquakes have occurred in the broader Crete region since antiquity and during the pre-instrumental era, until 1900 AD. The earthquake of 8 August 1303 [15,22], with an estimated magnitude of 8.3 [23], is one of the largest that has ever been reported in the entire Mediterranean region. This earthquake also caused a large tsunami that was observed at Crete, Rhodes, Peloponnese, Egypt, Syria, the Adriatic Sea, and along the eastern Mediterranean coast. In the vicinity of the study area, two major historical earthquakes have been reported ( Figure 1). The first is the 1 July 1494 (M = 5.4) earthquake [15] that occurred close to Heraklion and caused damage to church towers and private buildings, while large waves were observed in the harbor [22]. One century later, on 26 November 1595, another earthquake (M = 6.4) took place a few kilometers off the Heraklion city, causing severe damage and destruction in the entire Crete island [15,22]. It is worth noticing that no large events (M ≥ 6.0) occurred in the study area during the instrumental period (after 1900), with the nearest one being recorded on 23 May 1994 (M w = 6.1) and that was an intermediate depth (60)(61)(62)(63)(64)(65)(66)(67)(68)(69)(70) km) event ( Figure 1b) located approximately 40 km WNW of Heraklion [16]. Tsapanos [24] assessed the seismic hazard for the three largest cities of Crete and obtained maximum peak ground acceleration (PGA) values for Heraklion between 0.130 g and 0.165 g for rock and soft soil conditions, respectively.
One of the major structures that seem to control the HB-MB complex is the NNE-SSW fault zone that controls its eastern margin, located only a few kilometers east of Arkalochori village ( Figure 2). It bounds the easternmost post-Alpine sediment outcrops of the HB-MB in a morphologically clear way, even though several active structures seem to intersect. It was only recently (27 September 2021) that the lethal "Arkalochori earthquake" occurred in this area, causing one fatality, and resulting in hundreds of damaged buildings, mainly of old stone-built masonries. As the authorities reported, more than 6150 houses have been declared uninhabitable due to excessive damage in the mainly affected municipality of Minoa Pediadas (Figure 1), as a result of the mainshock and the major aftershocks.
The earthquake of 27 September 2021, 06:17:21 UTC, occurred near Arkalochori, Central Crete,~25 km SSE of Heraklion city [25,26]. Its focal mechanism is characterized by an SSW-NNE to SW-NE-trending, nearly dip-slip normal faulting. Considering the solutions proposed by different agencies (Table S1), except for a few outliers, its strike generally ranges N200 • E-N230 • E and its dip angle varies between 40 • and 60 • . The magnitude of the main event was initially determined as M L = 5.8 by the Seismological Laboratory of the National and Kapodistrian University of Athens (SL-NKUA) and the Geodynamic Institute of the National Observatory of Athens (GI-NOA), but later measurements of seismic moment upgraded its value to M w = 5.9-6.0. The largest aftershock (M w = 5.3) occurred on 28 September 04:48:08 UTC, one day after the mainshock, with a similar focal mechanism. According to the report of ITSAK [27], the recorded PGA at station ARK1 in the epicentral area (Arkalochori) was 609 cm s −2 in the horizontal component (N-S) and 806 cm s −2 (~0.82 g) in the vertical one, with the duration of strong ground motion (>0.1 g) being 6 s. The horizontal to vertical spectral ratio reached 8 at the dominant frequencies of 0.65Hz and 1.4 Hz. The accelerometric recordings at station ARK1 far exceeded the provisions of the National Building Code [28] for effective PGA = 0.24 g at the epicentral region in the seismic hazard zone II and for soil classes A and B.
Herein, we present a multi-disciplinary work, based on extensive fieldwork and detailed geological mapping of the affected area at central Crete, combined with results from seismological and geodetic imaging techniques. Our aim is to understand the contemporary seismotectonic regime of a densely populated area close to the active Hellenic arc and provide useful data for a realistic model concerning the neotectonic evolution of Crete. We intend to examine whether arc-parallel extension is a temporal phenomenon that is related to the Arkalochori earthquake or a dominant contemporary stress component.  Figure 1b. Black lines show the mapped tectonic structures. Seismicity of the period from 2008 to May 2021 is from GI-NOA. The Arkalochori earthquake is represented by a star with bold outline. Note the offset beachball linked with a dashed callout line representing its focal mechanism solution from GI-NOA (see Table S1 at Supplementary Materials) and the faults, described in this work, are drawn in red (KF: Kastelli Fault, NF: Nipiditos Fault, LF: Lagouta Fault). Additional focal mechanisms are from [29].

Geological Setting
The study area is located at central Crete ( Figure 1b) and, specifically, at the easternmost part of the supra-detachment MB, which is bounded by active fault systems producing high relief landscape [3]. It is the area where HB shares significant structures with MB and, consequently, it is quite complex to distinguish the margin between these Neogene basins with a single discipline method. In the following sections, we present interpretations on the geological setting of the study area, based on geomorphology, lithology, and seismotectonics.

Geomorphology
The general morphology of the wider area is rather mountainous, with steep slopes which very often exceed the order of 45 • , due to active tectonic structures that divide it into fault blocks that move either independently or in groups, depending on the local dominant tectonic regime [3]. The latter seems to have been changing since Miocene and is either expressed with the change of sediment facies or imprinted on the drainage network patterns, which are developed within a hilly area, with tilted monocline blocks, between higher relief mountains, implying a listric geometry of the involved faults (e.g., [30][31][32]). The highest mountain is Dikti Mt which exceeds an elevation of 2100 m and bounds both HB to the east and MB to the N-NE ( Figure 2). The much lower elevation E-W-trending Asteroussia mountain chain (~1200 m) is the geomorphological structure that acts as the southern margin of MB ( Figure 2).

Lithology
The geology of the area consists of Alpine (mainly Triassic to Late Jurassic limestones/dolomites, flysch, marbles, and ophiolites) and post-Alpine formations (silts, clays, conglomerates, and sandstones of Middle-Upper Miocene). In many places, they are covered by Quaternary alluvial sands and conglomerates of small consistency, as a result of the erosion of the aforementioned formations. More specifically:

•
The Alpine basement formations can be distinguished into two major groups, regarding their tectonic emplacement on the Cretan nappe pile [1]. They are members of the tectono-stratigraphic succession which builds the island of Crete as a stacked pile of 10 geotectonic units with a total thickness of about 10-12 km [33]. The entire nappe pile behaves homogenously till after the Middle Miocene, when the compressional, almost N-S-trending, orogenic stress field is converted into an extensional one due to the African plate slab roll-back [9,[34][35][36], which happened at the edge of the Aegean microplate [37]. Within this context, Crete is a large part of the Mid-Miocene Southern Aegean core complex exhumation that took place between 24 and 15 Ma [9] which hosts several supra-detachment depositional areas [3,4], with MB being one of them.

•
The Cretan nappe pile is segmented by many faults and covered by sedimentary formations representing different paleo-environments, as follows (from deeper to shallower): • Viannos Formation (Mid. Miocene). Sandy/silty deposits of lacustrine origin, coexisting with alluvial conglomerates. The average thickness of the formation is about 400 m. Its outcrops may be found at most of the eastern MB margin, either covering unconformably the basement rocks or in tectonic contact with them through ruptured fault zones. • Skinias Formation. It is a transitional succession consisting of Middle Miocene silts, with a thickness of~200 m, yielding the sea intrusion within the MB area. At the area between Skinias and Martha villages, the Skinias silts were deposited in a deep (over 200 m) marine environment. The base of the deposits consists of Pelites, whereas layers of conglomerates and sand have been observed in the area close to Martha village. • Ampelouzos Formation. It is deposited unconformably above both aforementioned formations. The lower part of the formation includes a variety of sedimentary, well-layered deposits, transitioning gradually from shallow to deeper (fan conglomerates to homogenous continental shelf sandstones). According to fossils found in the continental and shallow-marine deposits, the formation is dated as Middle to Upper Miocene [38]. • Alluvial deposits. Unconsolidated sands and conglomerates with pebbles originated from all the above formations, mainly along riverbeds.

Seismotectonic Setting
The vast majority of the large events affecting Crete, before the 2021 Arkalochori sequence, have occurred offshore, either to the north or to the south, with the latter probably related to the intra-crustal graben system, i.e., the Ptolemy and Pliny trenches [39] ( Figure 1b). The only knowledge regarding the local seismic activity in the study area comes from a temporary seismological network, which consisted of 10 analog stations, installed in the area of Heraklion during the period of September to December 1995 [29]. Data recorded by this network revealed the existence of shallow onshore activity. During that period, the seismicity was mainly concentrated along the eastern margin of the Heraklion Basin, decreasing from east to west, with epicenters also located in the vicinity of Arkalochori. Furthermore, a W-E-trending decrease in the sedimentary cover has been identified offshore, north of Crete [40].
In a broader view frame, the location of Crete Island, north of the Hellenic Arc (Figure 1a), comprises an area of complex geodynamics, with transpressional tectonics dominating the south of Crete, in contrast with normal and oblique-slip faulting onshore, yielding a heterogeneous stress field [14,41]. On the other hand, an extension has been observed to the north, between Santorini and Crete [42][43][44][45]. The latter reveals extension in a general E-W direction, with NNE-SSW shortening observed in western Crete [46]. According to [14], the sub-horizontal S 1 axis is in an SW-NE (in the east) to the N-S direction (in the west), likely affected by the proximity of the region to the subduction interface and compression parallel to the direction of convergence, which becomes more oblique to the arc towards the east. However, the sub-horizontal S 3 axis also reveals E-W (in the west) to WNW-ESE (in the east) extension onshore, even in the region where the 2021 Arkalochori earthquake occurred, in accordance with the observed normal fault scarps, trending in a general N-S direction [19]. The resulting regional shear stress favors SW-NE left-lateral faults, mainly offshore, in agreement with [41,47]. On the other hand, fault sources, determined by [11], indicate mainly E-W normal faulting to the north and south of Crete, incompatible with the resolved stress. It should be taken into account that the aforementioned stress field was determined using offshore earthquakes at crustal depths, both to the north and the south of Crete, as no large onshore events were located onshore until recently. The stress field in the region of Crete was characterized as transpressional by [14] due to the involvement of shallow events related to strike-slip (primary stress state) or reverse faulting (secondary state). In central Crete, due to the absence of focal mechanism data for shallow onshore events, the stress field is largely a result of interpolation from resolved nodes of neighboring regions. The stress shape value of 0.5 indicates instability in both S 1 and S 3 axes. However, the direction of S 3 seems to be generally stable in a WNW-ESE direction in the region from central Crete to the volcanic arc in the north. This complex stress regime is also revealed by the focal mechanisms, as all types (normal, strike-slip, and reverse) can be observed offshore, either to the north or to the south ( Figure 1). Focal depths increase towards the north, due to the African slab rollback phenomenon. Onshore focal mechanisms at shallow depths are few. De Chebalier et al. [48] worked on the western edge of Crete and obtained fault-plane solutions of shallow events, revealing either approximately N-S normal or strike-slip faulting. Delibasis et al. [29] constrained 29 fault-plane solutions close to the 2021 Arkalochori earthquake sequence, which also indicate normal or reverse motion with, in some cases, a significant strikeslip component.

Materials and Methods
We combined data from different sources, performing a multidisciplinary study, to further understand the geotectonic regime that caused the 2021 earthquake activity and possibly reveal its association with the mapped fault pattern. Therefore, our methodology included field mapping accompanied by the necessary laboratory work (e.g., GIS techniques for tectonic geomorphology, photogrammetric processing of unmanned aerial system images), seismological data analysis and interpretation, as well as SAR interferometry processing.

Field Mapping
The fieldwork included detailed mapping of the outcrops, especially the post-Alpine sediments, as well as of the fault surfaces that crosscut the stratigraphic contacts between them. At several places, a drone was used for image acquisition, followed by photogrammetric processing and the generation of high-resolution digital surface models and ortho-photographs.
A quite clear marginal structure, trending NNE-SSW, which delineates the eastern side of the MB-HB, can be identified, even on a medium resolution shaded relief of the study area in central Crete ( Figure 2). This represents the major Kastelli fault zone, which constitutes the western boundary of Dikti Mt that lies to the east [3]. One of its segments, just uphill the new international Kastelli Airport (Figure 2), was captured by the UAS and modeled through structure-from-motion techniques, aiming to map its trace with the highest possible detail, especially at areas where the morphology is relatively smooth. About 320 images were processed and a point cloud of more than 100 million points was generated for the construction of a photorealistic 3D model, with a spatial resolution of 2.5 cm ( Figure 3). The mapped faults were measured along with the kinematic indicators that were found on their surfaces and interpreted within the frame of kinematic analysis. The classification and clustering of the field measurements were performed using algorithms and methods already published and incorporated in software packages, such as M.I.M. v.4.0.1 [49].

Seismological Data and Methods
The 2021 seismic crisis in the wider area of Arkalochori began in the form of an earthquake swarm in early June 2021. The situation, however, changed dramatically after the occurrence of the M w = 6.0 mainshock of 27 September 2021. Although several permanent stations of the Hellenic Unified Seismological Network operate on the island of Crete, at the time of the mainshock, the closest station to the epicentral area was KNSS of the Hellenic Mediterranean University (HC network at formerly Technological Educational Institute of Crete [50]), at a distance of about 12-22 km ( Figure S1). As a result, the seismicity of the first stage of the sequence (foreshock swarm) could not be determined with high precision, due to the lack of data from stations at local distances, which would help constrain the focal depths. Resolving hypocentral locations with data from more distant stations is further problematic due to the inhomogeneity of the velocity structure, because of the proximity of the area to the Hellenic subduction zone. The day after the mainshock, four temporary seismological stations (CRE1-4) were deployed by GI-NOA (HL Network [51]) around the aftershock zone ( Figure S1b). The inclusion of data from these stations drastically improved the hypocentral locations, even with the use of a generic velocity model. We collected the available catalogues, including manually analyzed phase data from SL-NKUA and GI-NOA at all operational stations, and compiled a dataset of 2500 events for the period between 1 June and 18 October 2021. We divided the data to be processed in two periods, one before (Period A) and one after (Period B) 28 September, when the temporary local network was deployed.
We used the HypoInverse code [52] to locate the sequence, initially with a generic model used by routine analysis at SL-NKUA, as well as with the local velocity model of [29] for Central Crete. The latter model was preferred, as it yielded lower root mean square (RMS) travel-time residuals and locations errors, while providing a better general image regarding the distribution of epicenters, in terms of scattering, especially between events that form spatial clusters. Another issue concerns the distance weighting parameters for the two periods, i.e., the maximum epicentral distance of stations to be used. After several tests, the maximum distance was set to~280 km for the first period and~120 km for the second one. A V P /V S ratio of 1.78 was measured using the Chatelain [53] diagram for the first period, whereas a V P /V S = 1.76 was preferred for the second period with the shorter distance limit. The application of station corrections, especially for the first period, plays an important role to reduce artifacts of the initial epicentral distribution, where spurious epicenters of small magnitude events seemed to spread northwards of the main seismicity cloud.
To further improve the hypocentral distribution, with emphasis on the second period, we applied a relocation procedure using the HypoDD code [54]. The algorithm works by minimizing the double difference between calculated and observed travel times in combinations of event pairs with neighboring hypocenters, at inter-event distances much smaller than their respective hypocentral distance from a station. This method enhances the hypocentral distribution by minimizing relative location uncertainties between correlated events, caused by discrepancies between the velocity model and the real structure. Furthermore, waveform cross-correlation data are used to minimize uncertainties due to arrival-time reading inconsistencies in groups of events with similar waveforms, called multiplets.
We used station KNSS ( Figure S1b) as the main reference station for the identification of multiplets, as this is the closest station that was operational during both periods, with the exception of its N-S component that was damaged after the occurrence of the 27 September mainshock. Recordings of the temporary station CRE1 ( Figure S1b) were also used to enhance the reference cross-correlation measurements, particularly those of smaller events, during the second period. The waveforms were band-pass filtered in the range 2-8 Hz, for KNSS, and 2-15 Hz, for CRE1, before cross-correlation was calculated for the full P and S waveform signals in all combinations of event pairs with available data. The crosscorrelation maximum (XC max ) of every pair was registered in a matrix for every component. The matrices produced for different reference stations and components were then combined, retaining the RMS value of the XC max determined for each pair. The nearest neighborhood linkage was applied and a C th = 0.81 optimal threshold was selected, maximizing the difference between the size of the largest multiplet and the sum of clustered events [55]. This procedure created 234 multiplets containing a total of 1725 events. The P-and S-wave windows for each pair of events within the same multiplet were then cross-correlated for all stations with available picks, measuring the XC max and its respective time lag to be used as input for HypoDD.
The relocation procedure was applied separately for the two periods, acknowledging issues due to the different network geometry, so as not to degrade the data quality of the second period. On the other hand, we chose to relocate the mainshock and its early aftershocks (roughly a hundred), which occurred before the deployment of the temporary local stations, together with events that were located with data from the local network. This manages to improve the relative locations of the first aftershocks by exploiting their links with the rest of the events of the second period. Each relocation procedure was mainly divided into two types of sets, one with stronger a priori weights to the catalogue data and the other with stronger weights to the cross-correlation data.

SAR Data and Interferometric Processing
Space-born synthetic aperture radar (SAR) imagery is routinely utilized for measuring co-seismic surface displacements based on the DInSAR technique. DInSAR allows mapping the co-seismic motion by using two satellite images taken before and after an earthquake [56][57][58]. The systematic availability of SAR data from the Copernicus Sentinel-1 mission, with a revisiting time of 6 days over Europe (12 days globally), as well as the wide swath coverage of 250 km, enables the near-real-time response to seismic events [58][59][60][61][62][63].
The contribution of platform-based solutions in the early response phase of an earthquake has been well demonstrated [64,65]. The Geohazards Exploitation Platform (GEP; https://geohazards-tep.eu, accessed on 23 January 2022) focuses on mapping hazard-prone land surfaces and monitoring terrain motion, providing access to several Earth Observation (EO) missions and a broad range of relevant online services and development on cloud processing resources [66].
The present work utilized Sentinel-1 data acquired during the period from 23 September to 1 October 2021, with different viewing geometries to optimize the displacement map generation. For the interferometric processing, the GEP DIAPASON service was exploited, considering precise orbit state vectors and SRTM Digital Elevation Model 3arcsecond (≈90 m) data to remove the topographic phase. Interferometric measurements correspond to movements detected along the Line-of-Sight (LoS) of the satellite, i.e., the oblique direction between the satellite and the Earth's surface. Positive LoS values indicate uplift or motion towards the satellite, whereas negative ones correspond to subsidence or motion away from the sensor. The DInSAR processing scheme implemented herein is described in detail in [67,68]. Four differential interferograms were generated for the descending 036 (D036) and 109 (D109) orbital tracks, as well as for the ascending 102 (A102) and 029 (A029) ones, with temporal separation of 6 to 12 days and spatial resolution of approx. 45 m. Figure 4 and Figure S2 in the Supplementary Materials show the obtained displacement fields, while Table 1 details the interferometric pairs considered for the various geometries. Based on their temporal coverage, it is important to note that each co-seismic displacement map is affected by a different contribution of post-seismic motion. The interferograms for tracks A102 and D036 represent those with the smallest and largest post-seismic contribution, respectively (Figure 4). October 2021) (c,d) and the corresponding wrapped differential interferograms (e,f). The spatial extend of Sentinel-1 acquisition frames for both geometries are also indicated (a,b). The difference in the temporal span following the mainshock (27 September 2021) for the two opposite geometries is worth noting. While A102 spans approx. two days, D036 extends over five days after the earthquake, denoting larger contribution of post-seismic motion. Black lines correspond to fault zones (after [3]). The relocated epicenter of the 27 September 2021 mainshock is marked with a star. Furthermore, due to the opposite geometries, independent LoS measurements were decomposed into vertical (Up) and E-W-motion components, to facilitate the optimal interpretation of motion patterns. The horizontal and vertical components were calculated using both Sentinel-1 geometries in a comprehensive decomposition scheme [69], considering the variability of the incidence angles within the scenes. In order to simplify the geometrical assumptions, only measurements of similar incidence angles were considered in the decomposition (i.e., orbital tracks A102-D036 and A029-D109) (see Table 1).

Data Analysis and Results
The combination of data from various and heterogeneous sources is a complex task that needs to be evaluated and prioritized to examine their compatibility grade. The field data analysis provides us with strong arguments for establishing plausible planes of movement, which host the hypocenters of an earthquake. Therefore, it offers a tool for filtering them with higher relevance to the contemporary tectonic regime. In addition, the field data were used to validate the surface deformation detected through DInSAR data and attribute it to specific geological structures (e.g., active faulting).

Field Data Analysis
The length of the Kastelli fault zone exceeds 30 km since it can be traced from the northern coast of Crete to Asteroussia Mt to the south, within the recent sediments of MB. It was mapped and measured at several places, especially where the Alpine basement outcrops hosted the fault plane. It was quite frequent that striations along with other kinematic indicators were overprinted on the calcite coating of the fault plane ( Figure 5). The fault slip data recorded in the field were used as inputs within the stress inversion method TRM [70] and a stress tensor with S 1 : 90/079, S 2 : 00/349, and S 3 : 00/079 (R = 0.25) was defined. The resolved stress tensor, although calculated with fault planes that do not deviate strongly from each other, is rather compatible with the E-W-trending extension, which in turn agrees with the Arkalochori earthquake focal mechanism.
During the fieldwork, it was rather obvious that the morphological discontinuity, which the Kastelli fault zone has generated through its activity, is gradually eliminated and finally disappeared within the early Miocene lacustrine sediments of the Viannos, Skinias, and Ampelouzos units [71] and, hence, the fault trace which comprises its southern segment (Lagouta fault, according to [3]) is not very clearly delineated.
Despite this, a significant number of smaller fault traces that were mapped around Lagouta and Kassani villages (Figure 6a) yield that the Kastelli fault zone continues to the south and ends up to an E-W-trending fault, almost parallel and antithetic to the south Cretan detachment [3]. In addition, a pair of conjugate syn-sedimentary faults was located, measured, and analyzed ( Figure 6b) from a kinematic and dynamic perspective. The measurements of the two faults (p1 88/114, p2 72/150) and the sliding directions recorded on them (s 1 27/025, s 2 19/066), respectively, were found to be compatible-after tectonic analysis-with left-lateral strike-slip faulting and, specifically, with a stress field which is defined by the axes S 1 : 24/046, S 2 : 63/200 and S 3 : 11/311 [3]. The left-lateral component of this movement is in agreement with the en-echelon formation of numerous faults found along the previously activated E-W-trending, MB marginal fault zone, at the southern edge of this segment at Asteroussia Mt [3]. The monoclinal stratigraphy of the area yields that the same west-dipping fault zone has also a strong normal component. This conclusion is based on the fact that the outcrops of the older post-alpine sediments (Viannos formation) are found beneath more recent successions that have been mapped in detail on the western fault block and simultaneously at higher altitudes at the eastern one, yielding that the latter has been relatively uplifted (Figure 7).

Seismological Results
Regarding the initial locations of earthquake hypocenters (Figure 8a) of the first period, before 28 September 2021, focal depths remain uncertain, due to the lack of local stations, mainly ranging between 3 and 20 km ( Figure S3c), with the majority of events being concentrated near a discontinuity of the velocity model at~6 km [29]. On the other hand, more constrained focal depth values, between 5 and 13 km ( Figure S3d), are determined for the second period. Location statistics for the two periods are presented in Table 2 and Figures S3 and S4. Average vertical location errors are nearly halved during Period B and horizontal location errors are reduced by~26%, as a result of the deployment of stations at local distances, which also reduced the azimuthal gap ( Figure S4).  In order to assess the 2021 Arkalochori earthquake sequence in more detail and to identify the activated tectonic structures, a relocation procedure was applied, as previously described in Section 3.2, to both periods (Figure 8b). The events distribution was divided into four groups. Group #1 (red) is a temporal cluster that concerns the first period of the foreshock swarm activity, whereas groups #2, #3, and #4 (cyan, blue, and green, respectively) were determined by spatial clustering of the events of the second period. This was achieved by applying Ward's linkage to the matrix of inter-event epicentral distances of the relocated hypocenters [72]. This procedure of agglomerative hierarchical clustering, also known as the "minimum variance method", fuses two clusters into a new one, provided that the resulting increase in the sum of squares (i.e., its objective function) between the objects of the new cluster is minimized, compared to the outcome of any other potential combinations of cluster fusions. After creating the clustering hierarchy with Ward's linkage by setting an appropriate threshold to the fusion level, we divided the spatial distribution into three clusters (i.e., groups #2, #3, and #4) to aid the description. The evolution of the earthquake sequence is examined through the spatiotemporal projection of Figure 9, along the SW-NE-oriented, 16 km-long profile A-B of Figure 8b. The foreshock swarm began in June 2021 with a few small events (M L ≤ 3.0) and an M L = 4.2 event on 4 June which did not, however, cause any significant outbreak. Minor bursts occurred on 8 June and then on 18 July including three M L = 4.1-4.2 events. During that period, a slow migration of seismicity towards the southwest can be observed in Figure 9. Then on 24 July 02:07:37 UTC, an M L = 4.8 occurred, triggering a subsequence. Afterwards, apart from a few more minor bursts, the seismicity rate gradually diminished, notably to the point where only 17 events were located in the period between 16 and 25 September, with a complete absence of detectable events for about 41 h preceding the M w = 6.0 mainshock of 27 September 2021, 06:17:21 UTC. The whole foreshock group #1 is roughly distributed in an area of 3.0 km × 1.7 km, elongated in an N-S direction, with its centroid~2.5 km WNW of Arkalochori. The main event's epicenter is located adjacent to the foreshock swarm, near its eastern edge.
The occurrence of the mainshock, near the middle of profile A-B, caused seismicity to spread rapidly nearly through the whole length of the aftershock zone, both towards the southwest (group #4) and the northeast (group #3). The aftershocks appear to have little overlap with the foreshock zone, specifically with the sparser middle group #2 (cyan) and the northern margin of the southern group #4 (green). The two major groups, #4 south-west of the foreshocks and group #3 (blue) towards the NE, the latter with epicenters near the Kastelli airport, are separated by the less populated group #2 (cyan). This deficit of aftershocks between groups #3 and #4 is likely associated with an area of maximum coseismic slip during the mainshock. The rupture of a large asperity in that region apparently caused redistribution of stress towards the northern and southern margins of the main fault, while the central part was relaxed due to stress drop. The largest aftershock was an M w = 5.3 event that occurred on 28 September 04:48:08 UTC, about 4 km west of the mainshock. The epicenter of the main event is located at the NE edge of group #4, whereas the largest aftershock occurred towards its western part, near the southwestern tip of the foreshock group. The major aftershocks include 11 more events with M L ≥ 4.0 between 27 and 29 September, with only smaller events being detected through the rest of the sequence during the period of study (up to 18 October 2021). Regarding the period after 18 October, routine monitoring of seismicity at SL-NKUA detected five more events with 4.0 ≤ M L ≤ 4.5 between 20 and 22 of October, two associated with cluster #3 and three with cluster #4, while another M L = 4.2 event occurred on 29 December associated with the latter cluster. However, no significant changes in the overall spreading of the aftershocks' distribution were observed, besides a small, isolated cluster that was detected~10 km NNW of Kastelli, activated between 16 and 18 January 2022, with the largest event having M L = 3.6.
The three groups of the relocated aftershocks, if considered as a whole, appear to be well distributed on a plane, striking~N216 • E and dipping~45 • (Figure 10, red-dashed line), as determined by the least-squares fit on the hypocenters. This is quite compatible with focal mechanism solutions for the mainshock reported by most agencies, taking into account the WNW-dipping nodal plane (Table S1). Regarding specific spatial groups of the aftershock sequence, the best-fit plane for the northern cluster #3 strikes N181 • E and dips 49 • westwards ( Figure 10, blue dashed line in profile c 1 -c 2 ), whereas the mid-southern part (groups #2 and #4) has a least-squares plane that trends N216 • E and dips 53 • WNW ( Figure 10, green-dashed line in profiles a 1 -a 2 and b 1 -b 2 ). The more north-southwards orientation of the northern cluster's plane makes it more compatible with the Kastelli fault. The fault is steeply dipping at the surface (74 • ), but the aftershocks' distribution indicates a lower dip angle at the hypocentral depths of 6-9 km, implying listricity of the fault. Cluster #3 is fairly more concentrated than group #4 in the south, which appears to include several smaller sub-clusters, some of which may be related to the activation of antithetic structures (e.g., Figure 10, black-dashed lines with a question mark in profile a 1 -a 2 ). Normal faulting is indicated for the mainshock, the major aftershock, and several other aftershocks and foreshocks. However, oblique and strike-slip faulting is also observed for several events, including some foreshocks ( Figure 10). Although the mainshock's focal mechanism is in line with the rest of the aftershocks, it is noted that its hypocenter, along with those of aftershocks that occurred before the deployment of the local network, is not well constrained at depth. Hence, it is reasonable to assume that its true focal depth should be deeper, i.e., between 6 and 14 km, which is an estimated range for its centroid (Table S1).  Figure 7b. The cross-section profiles are drawn in an N103 • E direction, roughly perpendicular to the strike of the mapped Kastelli and Lagouta faults, whose downdip extension is presented in the cross-sections with gray-dashed lines. The red dashed line is the section of the least-squares plane fitted to the relocated hypocenters of period B, whereas the green (a 1 -a 2 , b 1 -b 2 ) and blue (c 1 -c 2 ) dashed lines are the planes fitted to groups #2 and #4 (green), and group #3 (blue), respectively. The possible antithetic structures are marked with black-dashed lines (a 1 -a 2 ). Focal mechanisms of the mainshocks and other major events of the sequence, including foreshocks (red beachballs, on the map only) are from the database of GI-NOA. The location of Arkalochori village is marked with a black square on the map.

DInSAR Analysis
In terms of satellite observations, the systematic acquisition strategy of the Copernicus Sentinel-1 mission ensured the availability of SAR data to investigate the earthquake sequence. Interferometric coherence levels for 6-to 12-days pairs are maintained, allowing for the extraction of reliable DInSAR ground motion measurements.
Independent results from four different satellite tracks (Figures 4 and S2) are consistent, indicating co-seismic downthrow of the epicentral area, extending over an area of about 38 km 2 to the east of the Lagouta and Kastelli faults. The ground displacement patterns are comparable, retaining an elliptical shape of the co-seismic fringes within the epicentral area. The increased ellipticity of the fringes for pairs covering a larger time span after the mainshock (tracks D036, A029, and D109) implies the contribution of post-seismic motion, compared to the more concentric fringes of the A102 track, containing only 2 days of postseismic activity. In addition, the direction of the major axis of the elliptical fringe pattern is following the aftershock distribution, retaining roughly an NE-SW trend. The spatial continuity of the interferometric fringes yields that the rupture did not reach the surface. The higher density of fringes at the eastern part of the epicentral area implies higher ground displacement gradients, and thus, proximity to the upper part of the activated rupture zone (i.e., NNE striking, west-dipping fault).
However, examining in more detail the DInSAR displacement fields for the various independent measurements, variability arises both in terms of magnitude and location of ground motion maxima (Figure 4). Specifically, for the pair with the shortest temporal extent (2 days) after the mainshock (A102, 23 September 2021-29 September 2021), LoS ground displacement is located in the vicinity of the relocated mainshock epicenter, with values reaching −18 cm. On the contrary, for the pair with the largest temporal extent (track D036, 5 days), the location maximum has been shifted approx. 2.5 km towards ENE and the observed ground displacement increased to −20 cm (Figure 4). The other two displacement pairs (A029 and D109), spanning 3 days within the post-seismic period, show comparable motion (downlift of about −20 cm) with the 5-days D036 pair ( Figure S2). The shift of the displacement maxima can be attributed to the amount of post-seismic motion contained in the DInSAR measurements (see Table 1), and it is consistent with the seismological observations indicating migration of aftershocks towards the northern margins (cluster #3) of the mainshock epicentral area. It can also be assumed that no surface deformation occurred, or at least with a magnitude detectable by DInSAR, after three days of post-seismic activity.
In fact, a few centimeters of variation of LoS displacement measurements (still within the error budget of the conventional DInSAR technique) might also be related to the difference in viewing the angle between the ascending and descending satellite observations, and the contribution of the horizontal (mainly E-W) motion component. Therefore, the interpretation of the observed LoS displacements is in advance favored by the decomposition of the ascending and descending interferograms (Figures 11 and S5). The estimated vertical motion indicates down throw reaching −22 cm (Figure 11a), concentrated within the mainshock epicentral region, whereas a more complex pattern, with a maximum of +8 cm, is presented by the E-W horizontal component (Figure 11b). Figure 11. Vertical (a) and E-W (b) ground displacement maps as decomposed using Sentinel-1 DInSAR LoS observations from ascending A102 and descending D036 tracks (see Table 1). Black lines correspond to fault zones (after [3]). The relocated epicenter of the 27 September 2021 mainshock is marked with a star.

Discussion and Conclusions
The recent Arkalochori earthquake sequence, which was recorded before and after the mainshock (27 September 2021), seems to be a quite prominent opportunity to discuss the contemporary tectonic regime of MB. The latter has a very significant geotectonic placement in the context of the African plate subduction that takes place under the Aegean microplate, as it is an onshore post-orogenic basin in the vicinity of the active Hellenic Arc [4,73]. Crete is at a crucial location right above the subducting slab and it is more than certain that the impact of the plate convergence is imprinted on the island, at all scales [74].
The fault pattern that was mapped in detail revealed the different generations of fault-block movements, as a result of the stress field changes since the integration of the Alpine nappe pile buildup [1]. Several direction changes of the main stress axes have been imprinted on the fault structures and the surface morphology [3]. It seems that the easternmost area of MB had been under N-S-trending compression during the orogen buildup, followed by an N-S-trending extension during the slab rollback and the orogen collapse through the detachments' operation (either north or south dipping) [1,8,10,75].
From an evolutionary point of view, we argue that after the generation of the nappe pile on Crete and the south-dipping detachment fault, several almost E-W-trending fault zones were activated (as segments of the south-dipping detachment) since they crosscut rocks that belong to different Alpine units. These fault zones remained active until after the end of late Serravallian, when the deposition of the lacustrine Viannos sediments was completed, as large fault blocks comprising Serravallian age successions are found tilted towards NNE.
A transtensional geotectonic regime, which seems to be dominating the entire Aegean microplate [76], is expressed on the Cretan territory by nearly NNE-SSW-trending strike-slip fault zones [20,77] and one of them crosscuts the central area of the island and simultaneously was acting as the marginal structure of the HB-MB system [2]. It should be mentioned that the main fault zones that were found along the eastern margin of MB are striking NNE-SSW, NW-SE, and E-W. The kinematic analysis of these data showed that most of the faults that were activated after the Early Miocene were either normal or strike-slip faults with left-lateral movement components [3]. According to this analysis, the NNE-SSWtrending shear zone was also affecting the previously generated E-W-trending fault zones, which were inherited by the Early Miocene N-S extensional period, when the south Cretan detachment was still active [1]. The impact on these E-W-trending structures is crucial as they seem to have been segmented and some of the segments have been reactivated during or shortly after the Pliocene (e.g., Nipiditos fault).
Specifically, during Tortonian, while the sediments of the Ampelouzos formation were still depositing, the dominating stress field had a strong left-lateral component, which is evidenced by the syn-sedimentary structure analysis (Figure 6b). The stress field seems to be active during the Tortonian since the syn-sedimentary faults were deforming the Ampelouzos unit strata (see Section 2.2 for lithology description), which were deposited before 10 Ma in a coastal to shallow depth marine environment [78]. A strong case scenario, considering both tectonic and stratigraphic data, is that the same fault zone seems to have a larger left-lateral component during Tortonian and the stress field gradually changed during Messinian, when the normal component dominates, increasing the compatibility of the stress field with arc-parallel extension. The latter seems to be active up to the contemporary era and is tectonically responsible for the Arkalochori earthquake.
The tectonic analysis of the collected measurements during the fieldwork provided strong arguments of stress field changes along the NNE-SSW-trending fault zones, as two groups of structures were found on fault surfaces revealing lateral movement and extension along the NW-SE direction. As explained above, the lateral movement on the same fault surfaces is older than the normal one, which in turn coincides with arc-parallel extension.
Finally, we refer to the Kastelli fault zone and in particular its southernmost segment (Lagouta fault), which is a west-dipping structure hosting a normal displacement movement with a strong left-lateral oblique component inherited by the dominating stress field during Tortonian. It was exactly this structure that was activated during the Arkalochori earthquake, as the hypocenters of the sequence clearly delineate the above-mentioned fault surface. It is obvious that the two hypocentral clusters that have been described, form two discrete segments that belong to the same fault zone. It is the WNW-ESE-trending Nipiditos fault that seems to be the structural reason for this segmentation. The northern cluster is formed on the surface of the Kastelli fault and the southern cluster on the Lagouta fault ( Figure 10). Even though all events belong to the same sequence, it is rather clear that the northern cluster is located at more shallow depths compared to the southern cluster.
Concerning the seismological data, the 27 September M w = 6.0 mainshock was a reminder that strong earthquakes do occur onshore Crete. It was preceded by a rich sequence of over 700 foreshocks that was considered a swarm, due to its largest event occurring halfway through 24 July. However, it did not cause much concern at the time, as similar swarms have been known to occur elsewhere in Greece, without leading to a strong event of M w ≈ 6.0, unlike the case of the Arkalochori earthquake. The herein analysis highlights the importance of local seismological networks to the reduction in location uncertainties. The applied double-difference relocation has managed to improve the relative locations of the foreshock epicenters, which are concentrated in the vicinity of the 27 September mainshock, although their hypocenters could not be adequately constrained. The mainshock apparently broke a large asperity of a west-dipping normal fault and distributed stresses towards its northern and southern edges, triggering aftershocks mainly at two large groups, separated by a spatial gap, where the asperity was located. Similar cases have been previously reported in other significant earthquakes on normal faults in Greece, including the 1999 Athens [79,80], the 2017 Kos [81,82], and the 2020 Samos [83][84][85][86] earthquakes. The herein presented relocation analysis of the sequence during its second period, i.e., with data from the temporary local network that was deployed by GI-NOA, permitted the detailed delineation of the main activated structures which could be related to mapped faults on the surface, allowing for a small degree of listricity, as well as contingent smaller antithetic, conjugate faults at depth. The spatiotemporal evolution of the sequence indicated triggering of seismicity throughout most of the aftershock zone soon after the mainshock, attributed to coseismic stress transfer, followed by slower migration towards its outer edges, indicating possible after slip.
Based on the post-earthquake fieldwork, both activated fault segments (Kastelli and Lagouta) did not show any signs of surface raptures along their traces, which agrees with the interferometric findings. A straightforward interpretation of horizontal motion patterns derived from the DInSAR analysis points out that the observed migration of motion during the post-seismic period (interferograms of different temporal spans) is mainly in agreement with the aftershock distribution and clustering, and to a lesser extent linked to actual co-seismic motion patterns. This is more evident by the collocation of the horizontal motion lobes to the aftershock clusters #3 and #4, as well as the absence of motion in the central part of the epicentral area where the mainshock nucleated. The vertical down-throw at the Arkalochori village area is also in accordance with the local stress field, being attributed to the hanging wall of the normal westward dipping Lagouta and Kastelli faults ( Figure 6).
The activation of the specific fault zone is also implied by the increased displacement gradients at the easternmost part of the epicentral area. However, the −22 cm of vertical motion, based on the decomposition of opposite DInSAR geometries (Figure 10), appears to be contaminated by post-seismic deformation. A more conservative estimate would attribute −18 cm of LoS motion to the main earthquake event, as calculated by the interferometric pair of the shortest temporal span after the mainshock. Additional displacement detected by other DInSAR pairs, also present in the decomposition, is attributed to the migration of motion towards NE, in accordance with the clustering and the shallower aftershocks' depth nearby the northern Kastelli fault segment.
In conclusion, the 2021 Arkalochori earthquake is a characteristic implication of the contemporary existence of a local arc-parallel extensional regime adjacent to major arcnormal strike-slip zones on the forearc region of the Aegean microplate. The consistency between seismological, geodetic, and field observations has been well demonstrated, highlighting the complementarity of multi-disciplinary approaches. The availability of dense seismological recordings and the contribution of Earth Observation platform-based solutions, validated, when possible, with field observations ensure proper mapping and tectonic interpretation of seismic events.