The Arkalochori Mw = 5.9 Earthquake of 27 September 2021 Inside the Heraklion Basin: A Shallow, Blind Rupture Event Highlighting the Orthogonal Extension of Central Crete

Abstract

A strong, shallow earthquake occurred near Heraklion (Crete, Greece) on 27 September 2021. The earthquake produced significant ground deformation in the vicinity of Arkalochori village but without any evidence for surface ruptures of primary origin. We used geodetic (InSAR and GNSS) data to map motions of the Earth’s surface that occurred during and shortly after the earthquake. A 14 cm subsidence of the GNSS station ARKL and a maximum of 19 cm distance from the SAR satellite were recorded. The measured surface displacements were used to constrain the rupture geometry and slip distribution at depth. Our best-fitting inversion model suggests that the rupture occurred on a 13 km-long planar normal fault striking N195° E dipping 55° to the northwest, with major slip occurring to the east and updip of the hypocentre. The fault tip is located 1.2 km beneath the surface. The maximum coseismic slip occurred in the uppermost crust, in the depth interval of 4–6 km. A decrease in the fault offsets toward the Earth’s surface is likely caused by an increased frictional resistance of the shallow layers to rapid coseismic slip. Satellite observations made in the first month after the earthquake detected no post-seismic deformation (i.e., below one fringe or 2.8 cm). The seismic fault may be identified with the Avli (Lagouta) segment of the NNE-SSW striking, west-dipping, 23 km-long neotectonic Kastelli Fault Zone (KFZ). Part of the rupture occurred along the Kastelli segment, indicating a fault segment linkage and a history of overlapping ruptures along KFZ. Based on geological data and footwall topography we estimate an average slip rate between 0.17–0.26 mm/yr for the KFZ. The Arkalochori earthquake is a paradigm example for the on-going extension of Heraklion basin (central Crete) in the WNW-ESE direction, which is almost orthogonal to the E-W Messara graben and other active faults along the south coast of Crete.

1. Introduction

A strong earthquake with magnitude ML_(NOA) = 5.8 (Mw_NOA = 5.9; USGS M_w = 6.0) occurred on 27 September 2021 06:17 UTC near the village of Arkalochori about 25 km SSE of Heraklion, Crete, (Figure 1; [1,2,3]). The event was shallow (the slip centroid of the NOA moment tensor solution was at 6.4 km; https://bbnet.gein.noa.gr/HL/seismicity/mts Event ID: noa2021sybnk accessed on 10 March 2022) and generated large vertical ground acceleration (0.8 g; [4]), that was recorded by a strong motion instrument in Arkalochori (Figure 1). One person was killed and several were wounded. Information on the impact of this event is available at https://en.wikipedia.org/wiki/2021_Arkalochori_earthquake (accessed on 1 December 2021). The earthquake occurred in a location of central Crete where information on previous seismicity is very limited [5]. Strong previous events within a distance of 20–30 km from Arkalochori were reported by Papadopoulos [5,6] to have occurred in 1508, 1779 and 1810.

The 27 September 2021 M5.9 event was preceded by a sequence of foreshocks that started to gain attention on 4 June 2021, when an event of magnitude M = 4.6 occurred (see Table S1; [1,2]). In fact, microseismic activity around Arkalochori had been recorded about one year before the mainshock (this study; [1]). Thus, the Mw = 5.9 27 September 2021, Arkalochori event occurred as a culmination of the foreshock sequence. The European Mediterranean Seismological Centre, EMSC, https://www.emsc-csem.org (accessed on 1 December 2021)) reported testimonies on intensities for several foreshocks (see Figure S1 and Table S1). In Table S1 and Figure S1 we present the testimonies, their number and location for the seven of them of magnitude larger than M = 4.0, the mainshock and the main aftershock of ML(NOA) = 5.3 that occurred one day after the mainshock on September 28, 2021 at 04:48 UTC. The distribution of severe damage was around the village of Arkalochori (Figure 1 and Figure S2).

Table 1. Moment tensor solutions calculated by various agencies. In the last line, the values of our best fitting geodetic model (see Section 4).

Agency	LAT	LON	Strike	Dip	Rake	Strike	Dip	Rake	Depth	Mw	Moment
	°	°	°	°	°	°	°	°	km		×1018 Nm
GCMT	35.12	25.18	40	43	−85	214	47	−94	12	5.9	1.01
USGS	35.25	25.27	55	57	−79	215	35	−106	12	6.0	1.08
ERD	35.39	25.35	35	47	−97	225	44	−83	10	5.7	0.76
INGV	35.10	25.18	22	53	−108	230	41	−68	10	6.0	1.2
GFZ	35.18	25.25	30	32	−96	218	57	−85	16	6/0	1.1
CPPT	35.24	25.37	20	68	−107	239	28	−55	16	6.0	1.23
IPGP	35.25	25.26	24	38	−111	229	55	−75	6	6.0	1.16
KOERI	35.12	25.21	44	41	−76	207	50	−101	10	6.0	1.13
OCA	35.13	25.26	11	45	−35	245	60	−56	10	5.9
NOA	35.14	25.26	25	33	−81	195	57	−95	6.4	5.9	0.86
UOA	35.14	25.27	21	25	−75	185	66	−97	14	5.9	0.91
Median			25	43	−96	218	50	−85	10		1.09
Mean			30	44	−95	218	49	−83	11		1.04
r.m.s. scatter			12	14	13	13	12	16	3.7		0.18
Model (this study)	35.15	25.28				195	55	−90	5.1	5.96	1.1

reports the various available moment tensor solutions for the mainshock. The mainshock hypocentre depth range reported in Table 1 is between 6–16 km. This difference appears as depth can be estimated from moment tensor inversions by fitting the spatial variation in surface or body waveform amplitudes. Then, the moment tensor inversion is performed for a range of focal depths and the best solution is assumed to be the one with the lowest misfit. However, as different institutes use different stations and different seismic frequencies it is evident that the depths will vary. We note that the NOA solution indicates a centroid depth of 6.4 km, which is close to the depth range of the main slip patch as determined in [2].

The location and depth of the 27 September 2021 M5.9 earthquake onshore Crete highlights the ongoing extensional deformation of the overriding Aegean Plate (e.g., [7,8,9,10,11,12,13,14,15,16]) in tandem with ongoing compression along the subduction front [17,18,19,20,21,22]. According to the latest GNSS results by Briole et al. ([23]; see their Figure 9) there is extension in the ~NW-SE direction of the order of ~0.5–1 mm/yr near Heraklion. This process is also manifested in the orientation of the local active faults and the topography of this region [9,13,16,24,25,26].

In this study, we present an analysis of seismological, geological, and geodetic data that constrain the location and geometry of the activated fault that ruptured with the M5.9 earthquake of 27 September 2021. In Section 3.1 we describe the method of relocation of hundreds of aftershocks. Due to the shallow depth of the earthquake rupture (hypocentre depths range from 6 to 16 km; see Table 1) and the good phase coherence of the area, it was possible to accurately map the surface deformation using InSAR (Section 3.3). The interferograms showed one main lobe of subsidence, elongated with an NNE-SSW orientation. The geodetic displacement data (InSAR & GNSS) were jointly inverted (Section 4) to model the dislocation source, assuming a homogeneous elastic half space. Our field survey showed that the earthquake generated secondary phenomena, with several areas on alluvial deposits exhibiting opening cracks. No tectonic surface ruptures were found in the field, in agreement with the modelled rupture. Our study highlights the WNW-ESE orientation of the stress field in this key region of central Crete and the particular type of strain release, which is by occurrence of earthquakes along normal faults that ruptured a large patch (10+ km long segment) of a longer (20+ km) neotectonic structure.

2. Tectonic Setting

The Heraklion basin is located in the central part of the Hellenic Arc and is characterised by extensional faulting, mainly along N-S and NW-SE directions [9,25,26]. Extensional features (intrabasinal highs, axial drainage etc.) are also manifested near the active faults [9,13,16,27,29] and the present-day topography of this region. The basin is found to the north of the Messara Basin, between Mts Psiloritis and Dikti in the west and east, respectively and Mt Asteroussia in the south (Figure 2). Both basins (Heraklion and Messara) are related to the back-arc extensional tectonic regime established in the Hellenic Arc since the Lower Miocene. The sedimentary fill in these basins includes a thick succession of 600–800 m of sediments which represent successive stages in their post-Lower Miocene evolution. The older basin sediments suggest lacustrine conditions in the Middle Miocene (Ano Viannos formation). Marine conditions were established in the Middle-Upper Miocene (Profitis Ilias, Ambelouzos and Agia Varvara formations); this marine environment continued in the Lower–Middle Pliocene (Foinikia fm). In the Pleistocene, the depositional environments are diversified, with marine bioclastic limestones and sandstones in the Heraklion Basin (Iraklio fm.) and fluvial–fluviolacustrine sediments in the Messara Basin [9,10,24]. The evolution of the Heraklion and Messara basins is thus characterized by diverse depositional environments, ranging from marine to fluviolacustrine and fluvial, since the Middle Miocene [30]. This alternation of depositional environments can be attributed both to climate-driven (eustatic) and tectonic factors, the interplay of which has led to the formation of numerous intra-basinal grabens and horsts.

The structural fabric within the Heraklion and Messara basins follows a pattern that reflects their evolution since their inception in the Lower-Middle Miocene; the former is dominated by NNE-SSW and NE-SW normal or oblique-normal faults, which appear to crosscut and overprint N-S and E-W faults that are believed to represent earlier deformation phases [9], including the south-dipping Psiloritis detachment [37]. The structures that control the geometry of the Messara basin, on the other hand, are mainly E-W and ENE-WSW faults. However, there appears to be considerable interaction and structural overprint between the two basins, at the eastern edge of the Messara basin (Kastelliana—Ano Viannos; Figure 2). This complex deformation pattern was also manifested in the 2005 Asimi and 2013 Houdetsi earthquake sequences, respectively, which confirmed the occurrence of both NE-SW normal and N-S strike-slip intra-basinal faulting, compatible with overall NW-SE crustal extension [25].

The village of Arkalochori is built on a topographic intra-basinal high, which is the result of both the occurrence of the erosion-resistant Agia Varvara bioclastic limestones and the fact that it is part of the footwall block of the N-S, west-throwing Arkalochori fault (see ALF in Figure 2), whose trace runs along the western edge of the hill. The cumulative throw of the Arkalochori fault is in the order of a few to several tens of meters, as is also the case of another minor fault in the east (Roussohori; see RF in Figure 2; this is fault GR2485 in the NOAFaults database http://0-doi-org.brum.beds.ac.uk/10.5281/zenodo.4304613 with a length of ~4 km). We have not studied these faults in detail, but certainly the geological data imply that their ages could be around 3–5 Ma, and assuming a preliminary throw of 100-m (implying a dip-slip of 140 m) the average slip rate is about 0.02–0.04 mm/yr. This value is expected for minor structures that do not accommodate most of the deformation. Both faults are characteristic of the ~N-S structural grain in this part of the basin and run sub-parallel to the marginal Kastelli normal fault (Figure 2).

The meizoseismal area of the 27 Sep. 2021, earthquake includes the villages Arkalochori, Avli, Nipiditos and Kastelli (Figure 1) including the aforementioned intra-basinal high, which corresponds to an eastward-tilted fault block, controlled by the intra-basinal N-S Arkalochori fault. The latter is one of a series of N-S to NNE-SSW normal faults that develop in the hanging wall of the major Kastelli Fault, which marks the eastern boundary of the Heraklion basin (KF in Figure 2). The SSW-ward segment of the Kastelli fault is the Lagouta Fault (LF in Figure 2), whose morphological expression is subdued, compared to the Kastelli fault. According to Vassilakis [27], this is due to the nature of the juxtaposed formations, as the Lagouta Fault juxtaposes the erodible footwall lacustrine marls (Ano Viannos fm.) against the Miocene marine sediments (Ambelouzos fm.) in the hanging wall. Nonetheless, this fault is characterized by a possible earlier strike–slip phase, in the Late Miocene, while later it acts as a pure extensional fault, in accordance with the suggested NW-SE extension in the broader area [25].

In terms of fault nomenclature, we note that in this paper, the neotectonic fault zone bounding the Kastelli Neogene basin is named as the Kastelli Fault zone (KFZ; following [9,26]) comprising two segments, the Kastelli Fault (KF in Figure 1) and the Avli-Lagouta Fault (AF in Figure 1). The normal fault in the footwall of KFZ is named the Geraki Fault Zone (GFZ in Figure 2). Two minor faults are also examined: the west-dipping Arkalochori Fault (ALF in Figure 2) and the west-dipping Roussohoria Fault (RF in Figure 2), respectively. In the literature, the above faults are named as follows: (1) in the NOAFAULTs v4.0 database https://arcg.is/04Haer and in Ganas et al. [25], the two segments of KFZ are named as Agnos (code GR2467), and Avli (GR2465), while the normal faults in the footwall of KFZ as Kastelli (GR2468; west dipping) and Panagia (GR2464; south dipping) (2) in Caputo et al. [13,29] the normal fault in the footwall of KFZ is identified as Kastelli Fault (3) in Nicol et al. [16] the KFZ is reported as Kastelli Fault, the normal faults in the footwall of KFZ as Geraki Fault, and Nipiditos Fault (4) in Vassilakis et al. [3] the two segments of KFZ are reported as Kastelli Fault and Lagouta Fault, the NW-SE normal fault at the footwall of KFZ as Nipiditos Fault and (5) in Triantafyllou et al. [2] the two segments of KFZ as North Kastelli and South Kastelli Fault; the normal faults in the footwall of KFZ as Geraki Fault zone (GFZ) and Nipiditos Fault.

3. Data & Methods

3.1. The 2021 Seismic Sequence of Arkalohori

In the present study, we have collected catalogue data (arrival times of events and seismic phases) from the online database of Institute of Geodynamics (National Observatory of Athens—NOA) for earthquakes that occurred in the study area during the period between November 2020 and October 2021. The mainshocks and aftershocks of the sequence were recorded by the broad-band and strong motion stations of the Hellenic Unified Seismological Network (HUSN; http://eida.gein.noa.gr/ (accessed on 1 December 2021); [38]); moreover, four (4) stations, named as CRE1-4, were installed on 28 September 2021, in the epicentral area by NOA (HL network). The available manually picked, P- and S-wave arrival-time data form a dataset of 1062 events with magnitude ML between 0.5 and 5.8. The local velocity model and the value 1.78 for the Vp/Vs ratio of Delibasis et al. [39] were adopted. The initial location was performed using the HypoInverse algorithm [40]. The applied distance weighting was set to nullify the influence of stations at epicentral distances further than ~100 km. The average absolute location errors reported by HypoInverse were: mean RMS travel-time error 0.14 s, ERH = 0.74 km and ERZ = 0.98 km for horizontal and vertical uncertainties, respectively.

Subsequently, to improve the initial events locations, the double-difference relocation HYPODD [41] procedure was performed. HYPODD determines relative locations within clusters, using a double difference algorithm, developed by [42]. It improves relative location accuracy by strongly reducing the influence of the velocity structure on locations and uncertainties caused by arrival-time reading errors. The algorithm minimizes the double-difference between observed and calculated travel times through iterative weighted least-squares using the conjugate gradients method. The relocation procedure takes into account both waveform cross-correlation and catalogue differential travel-time data to provide an enhanced picture of the seismic activity.

Multiplets were detected, for the closest stations, followed by the calculation of the cross-correlation maximum and time-delay measurements for P-waves, which were used as differential travel-time data input for the double-difference relocation algorithm, HYPODD [41]. Full waveforms were cropped in windows containing both P and S waves, filtered between 0.5 and 6 Hz, and a threshold of 70% was set concerning the waveform coherency. The process mainly implicated two groups and the results were then imported, along with catalogue travel-time data, into the double-difference for the relocation process. The relocation procedure is mainly divided in two stages, the first with stronger a priori weights to the catalogue data and the latter with higher weight to the cross-correlation data, as suggested by [41]. Distance thresholds were applied to ensure that no links were formed between unrelated events within the bigger groups. The velocity model used in the relocation was the same with the model used in the initial location process [39] and 19 stations within 100 km distance were used.

In total, 1056 events were relocated and clustered in two groups (Figure 3); the main cluster (1), oriented NNE-SSW) with 951 events that occurred at the vicinity of the mainshock and a smaller cluster (2), with 95 events, located towards the NE of the cluster (1). The HYPODD final results include the 98% of the initial dataset (absolute cross-correlation root-mean-square: 0.3860 s, absolute catalogue rms: 0.335 s). The relocated hypocentres have an estimated mean rms residual of 4 milliseconds (ms) and the mean location formal uncertainties x, y, z and t were 17.4 m, 19.3 m, 23.4 m and 5.6 ms, respectively. The clustering of the aftershocks is a robust feature of the aftershock sequence, activated in diverse times [2,3].

We made four cross-sections, oriented NE-SW and NW-SE in order to examine the geometry of the distribution of events. A NW-SE cross-section through the main cluster shows that the distribution of the hypocentres indicates a NW-dipping structure (cross-Section A1-A2, Figure 3 and Figure 4; as also imaged by [3]). We can infer the position of the fault plane dipping towards northwest between depths 5–15 km. The inferred dip direction is similar to the dip direction of the northwest-dipping nodal plane of the moment tensor solutions of NOA and GFZ (see Table 1; the median value of the northwest plane is 49°). The northern cluster (cross section C1-C2, Figure 3 and Figure 4) was activated at shallower depths (6–10 km;) however, the distribution of the hypocentres is not conclusive on the dip direction of the seismic fault.

Sparse events had already taken place at the same area since November 2020, while the activity seems to become denser by the end of May 2021 and more intense by the end of June 2021 (Figure 5); a majority of the foreshock events were restricted within the vicinity of the mainshock (±2 km; Figure 5 middle panel). In particular, the May–July 2021 events occurred immediately to the north (± 30°) of the mainshock within a distance of 4 km (Figure 5). The foreshock activity was intense; two earthquakes of magnitude Mw >= 4.0 took place while the largest foreshock (Mw = 4.6) occurred on 24 July 2021. A noticeable pause of seismicity appeared on September 26, 2021, followed by the mainshock during the 27th of September. The aftershock sequence was intense, while twelve of the aftershocks were of magnitude Mw > 4.0; the strongest one (Mw_(NOA) = 5.1; ML = 5.3) took place on 28 September 2021, in shallow depth, as well. Another characteristic of the aftershocks’ distribution was the bilateral occurrence of seismic events with respect to the mainshock hypocentre, already during the 27th of September 2021, including a tight cluster 5–10 km towards the NE of the mainshock (Figure 5). During the next few days, it was evident (Figure 5) that the northern cluster remained an isolated feature of the aftershock sequence. Its occurrence is due to either aftershocks occurring at the termination of the rupture or off-fault aftershocks occurring on a neighbouring structure (or even structures), indicating the strong interaction between the faults of the area close to Arkalochori.

3.2. GNSS Data Processing

Dual-frequency data from three (3) permanent (continuously recording) Global Navigation Satellite System (GNSS) stations were processed. We used GNSS observations covering several days before and after the earthquake (Figure 6) in order to obtain position time series and their static offsets due to the earthquake. The sampling interval was 30 s, and the data comprise 24 hr daily files. All station records were complete (rejected epochs 0.00%) and included in our processing, thus providing substantial observations for mapping the co-seismic displacement field. The stations belong to the HxGN SmartNET and Uranus private networks of Greece. The azimuthal distribution of the GNSS stations is very good, extending uniformly around the epicentral area (Figure 1). However, most of the GNSS stations are located more than 20 km away from the epicentre and therefore captured only small displacements of the order of a few mm.

The GNSS data of local continuous stations ARKL, HERA and MOI1 together with the data of four regional continuous stations (BSHM, DYNG, IZMI, MERS) were processed using the GAMIT/GLOBK/TRACK software version 10.71 [43]. The TRACK module is used to process kinematic GPS data. Our TRACK processing uses 30s resolution data from local ARKL, HERA and MOI1 stations. We obtained co-seismic and early post-seismic deformation (hours) only in the nearby ARKL station (Figure 6). The displacements of HERA and MOI1 stations (Figure 1) during the earthquakes were negligibly small. The presented data in Figure 6A–C show ARKL displacements relative to MOI1 station using the TRACK module. We also obtained a daily solution using the GAMIT/GLOBK moduli. GAMIT uses a priori data of satellite and station locations and obtains daily position and baseline solutions. GLOBK combines all data and creates a multi-day time-series solution for the position of stations. Our processing procedure included runs of both GAMIT and GLOBK. Then, outliers were detected and removed from the data. The final sets of locations are obtained with respect to the International Terrestrial Reference Frame 2014 (ITRF2014; [44]).

The only GNSS station located in the near field and showing large displacements is the station ARKL located just above the fault plane (Figure 1). The antenna of this station is located on top of an industrial building south of the village (Figure S3). Figure 6D–F shows the relative displacements of ARKL station days after the earthquake. The offsets indicate cm-size motion along the N, E and Up-down directions. Station ARKL (Figure 6) moved 7 cm to northward and 5 cm eastward. In terms of Up-Down (vertical) displacement, we obtained clear trends of co-seismic motion with amplitude −14.3 cm. The vertical downwards displacement at the GNSS station enhances the robustness of the west-dipping fault model as they allow discrimination between the two conjugate focal planes (see also forward model of Okada-type dislocations; [45]) in Figure S4. Regarding the large aftershock on 28 September. 2021, 04:48 UTC (Mw = 5.1; located 4 km to the W of Arkalochori, Figure 3), we could not find any evidence of this aftershock in the GNSS time series. Furthermore, the GNSS time series of the ARKL station show the existence of a small post-seismic motion (~2 cm of subsidence) following the mainshock with a relaxation period of 75 days (Figure S5). The processing of high-rate 1 s data (Figure S7; Brendan Crowell personal communication) show some motion preceding the co-seismic offset.

3.3. Geodetic Data from Sentinel-1 SAR Interferometry

Synthetic Aperture Radar (SAR) data from different directions are used to explore surface displacements due to the 27 September 2021 central Crete earthquake. SAR data that tightly bracket the earthquake dates were acquired using the Sentinel-1A and 1B satellites (see

Table 2. Characteristics of the interferometric pairs of Sentinel-1 data that were used in this study.

	Interferometric Pairs (yyyy/mm/dd)	Satellites	Relative	Incidence Mid Angle	Comments
			Orbit Number
Ascending	2021/09/17–2021/09/29	Sentinel-1A	102	44	Co-seismic (Figure 7)
Descending	2021/09/18–2021/09/30	Sentinel-1A	109	33.93	Co-seismic (Figure 7)
Descending	2021/10/01–2021/10/13	Sentinel-1B	36	43.87	Post-seismic (Figure S6)

) from the European Space Agency (ESA). Sentinel-1A/B satellites operate using C-band (radar wavelength of 56 mm) in the Terrain Observation by Progressive Scan (TOPS) mode, which allows for wide swaths and short repeat intervals. Each TOPS image consists of three sub-swaths, with a total swath width of ∼250 km, providing a complete coverage of the earthquake-struck region from two opposite directions. Image pairs that span the earthquake dates are available from both the ascending and descending satellite orbits. The Sentinel-1 radar remote sensing data was downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/ (accessed on 1 December 2021)) and processed using the Generic Mapping Tools Synthetic Aperture Radar (GMTSAR) package [46]. The topography contribution to the radar phase was calculated and removed using digital elevation data from the Shuttle Radar Topography Mission (SRTM) with 30 m resolution https://www2.jpl.nasa.gov/srtm (accessed on 1 December 2021). Accurate orbital data of the satellites are available via https://qc.sentinel1.eo.esa.int/aux_poeorb/ (accessed on 1 December 2021).

The differential interferograms (unwrapped) are shown in Figure 7A–D, which indicate clear coseismic displacements in line of sight (LOS) caused by the 2021 Arkalochori event. The negative LOS displacement indicates the ground motion away from the satellite. The main deformation zone depicts an elliptical shape with a size of 10 km × 7 km along the NNE-SSW direction (Figure 7A,D). The interferograms indicate a similar main feature of surface subsidence, although their maximum values are different because of their opposing satellite views. The fact that in both LOS imaging geometries the negative vertical component displacement is dominant, this implies a normal-slip event. The rupture may be related to the west-dipping normal fault on the eastern side of the Kastelli graben (see blue line in Figure 7).

4. Inversion for Coseismic Slip Model

The observed geodetic deformation can be modelled by slip along a rectangular dislocation plane in a homogeneous elastic half-space [45]. To assess the fault slip distribution during the earthquake we performed an integrated inversion model for the coseismic ascending and descending interferograms (Table 2) and displacements of the ARKL GNSS station (Figure 6). The inversions of near-field InSAR data from more than one direction are capable of better resolving shallow slip distribution. The applied inversion scheme [48,49,50] is based on a least-squares minimization of the displacement misfit with iterations for the fault geometry. The initial geometry trace is chosen using the linear phase discontinuities in the wrapped InSAR data. We refined the fault plane geometry by allowing changes in the fault dip, strike and location in order to minimize the data misfit in a least-squares inversion. To account for spatial variations in the slip distribution, the fault plane was subdivided into 1.5 × 1.5 km, constant-slip patches, allowing only normal slip mode throughout the inversion. The slip model was smoothed in order to minimize the slip gradients produced by fault. Finally, we choose the slip model that was the smoothest and with the lowest RMS value.

Our best-fit inversion model (Figure 8) accounts for 97% of the observed deformation from both InSAR interferograms and GNSS data with an RMS value of 1.3 cm (Figure 7). The model indicates coseismic dip-slip displacements during the 27 September 2021, an earthquake on a 13 km length of a fault plane that dips 55° to the west and strike of N195° E. The geometrical characteristics of the modelled fault plane (strike, dip-angle, dip-direction) better fit the moment tensor solutions of NOA/UOA/KOERI (see Table 1, ranging N185° E to N207° E) than those of global networks (GFZ, UGSG etc.). We attribute this difference to the use of regional network data (distances up to 500 km; periods less than 125 s) from the local networks than the use of teleseismic data used by global networks. The strike of the modelled fault N195° E differs by about 20° from the modelled fault of Triantafyllou et al. (2022; N216° E); the reason probably found in the method of geodetic inversion used.

The modelled fault parameters are presented in

Table 3. Parameters of the modelled fault based on the combined inversion of both InSAR and GNSS data.

Centre of the Top-Fault			Size of The Fault		Fault Orientation		Fault Slip
Long.	Lat.	Top	Length	Width	Azimuth	Dip	Rake	Max.Slip
°	°	km	km	km	°	°	°	m
25.309	35.138	1.2	13.4	10.5	195	55	-90	1.2

. The fault tip is located 1.2 km beneath the surface. Most of the slip is distributed between 1.2 and 10 km depth, with maximum slip of about 1.2 m at 4.9–6.2 km depth. The geodetic inversion results are compatible to the hypocentre of the relocated mainshock (5 km; Figure 4) and the aftershock distribution which is also concentrated in depth up to 10 km (Figure 4). The major slip at a depth of ~5–6 km coincides with a less dense concentration of aftershocks (Figure 4), implying a higher energy release by the mainshock at this depth. We note that the top edge of the modelled fault plane (line AB in Figure 8) aligns with the eastern margin of the Kastelli basin, which is fault-controlled [9,16] (Figure 9). Therefore, our inversion modelling suggests the activation of one normal fault, which may be identified as the west-dipping Avli (or Lagouta) normal fault (fault GR2465 in the NOAFAULTS v4.0 database; [47]; AF in Figure 1; LF in Figure 2), to the south of the Kastelli Fault (KF in Figure 1).

Moreover, the modelled slip above 1.2 km (Figure 8; patches from AC to BD) is negligibly small which indicates that the rupture did not reach the ground surface. The total calculated geodetic moment release is Mo = 1.1 × 10¹⁸ N·m (corresponding to a magnitude Mw = 5.96), which is comparable to the reported seismic moment release of Mw = 5.9 or 6.0 (see Table 1). With the uniform slip model, the fitting model displacement for both orbits and residuals are shown in Figure 7B,C,E,F, respectively). It shows that the model obtained by the inversion fits well with the InSAR data sets, and the residuals are also relatively small. A postseismic interferogram of the image pair 1–13 October 2021 (descending orbit; 3–16 days after the earthquake; Figure S6) shows no detectable signal for surface deformation, thus we conclude that the postseismic surface deformation during the 1–13 October interval is small (less than one fringe) compared to the coseismic signal. This result was confirmed by the small postseismic motion detected in the GNSS time series of station ARKL (see Figure S5).

5. Discussion

5.1. The Arkalochori Rupture Characteristics

The obtained slip distribution of the mainshock from the inversion of the geodetic data can be visualised in relation to the distribution of relocated aftershocks and the location of the hypocentre of the mainshock (Figure 9). The aftershock data are plotted in the perimeter of the mainshock slip patches carrying the largest amount of slip, at depths larger than 4–5 km (the NE cluster is deeper than 6 km, see cross-section C1-C2 in Figure 4). Thus, the aftershock pattern is arranged around the main rupture area (see red and brown patches in Figure 9) while most aftershocks occurred mostly down-dip of the slip patches (at depths larger than 8 km). In particular, we observe that at depths of 3–6 km there is a seismicity gap that extends about 5 km in the N-S direction, i.e., along strike (Figure 9). This gap may be due to the complete rupture of the main asperity along the fault plane and subsequent aftershock activity triggering around the co-seismic slip patches (e.g., [51,52,53,54,55]). In addition, a small post-seismic deformation (~2 cm on the vertical component) was detected on the GNSS data (Figure S5) but not on InSAR images (Figure S6). The small post-seismic signal may be due to afterslip along the fault plane and not to a gravitational instability of the ground hosting the GNSS antenna. Unfortunately, we do not have additional GNSS data to model the post-seismic deformation and related afterslip.

Given the shallow depth of the hypocentre (~5 km; Figure 4) in relation to (a) the earthquake magnitude (M5.9) (b) the high dip-angle of the fault, well-constraint by geodetic data (55°; 53° obtained by [2]) and (c) the shallow depth of the main slip patch (~5 km; Figure 8) it was reasonable to expect that surface ruptures could be found in the area towards the east of Arkalochori, near the villages Avli and Nipiditos (Figure 9 and Figure 10; allowing for small uncertainties in the intersection of the seismic fault with the Earth’s surface). However, our field surveys around Arkalochori Kastelli, and Nipiditos (Figure 10) did not show any evidence of tectonic ruptures, i.e., formation of fresh scarps with normal step in colluvium along the mountain front or reactivation of limestone scarps such as those found in several shallow events in Greece, e.g., in Thessaloniki 1978 [56], Pissia—Kaparelli in 1981 [57], in Kalamata (1986; [58]), in Atalanti (1894; [59]) etc. Instead, we found a lot of secondary cracks (both on natural and artificial surfaces, i.e., tarmac roads) and rock-falls due to the strong ground motion. This lack of shallow coseismic slip on the fault plane (uppermost 1–2 km) may be due to velocity-strengthening fault frictional characteristics [60,61] of the fault plane as it emerges towards the Earth’s surface. A decrease in the fault offsets toward the Earth’s surface is likely caused by an increased frictional resistance of the shallow layers to rapid coseismic slip that allows it to arrest the rupture. This is not always the case, off course, as there are numerous examples of tectonic surface breaks along normal faults [56,57,58,59,62,63,64] which indicates the influence of other factors, among them most important the size (magnitude) of the event related to the amount of elastic strain accumulation over the geological time.

In relation to previous studies on this event we note that [2] used geodetic data inversion (only InSAR) and catalogue events for their analysis as opposed to our data (InSAR, GNSS and relocated events). As a conclusion, they obtained a double-fault structure with listric geometry (their Figure 12) which is not supported by our findings (one planar fault; Figure 9). The study by [3] used geological data and seismicity relocation to model the seismic fault. They processed InSAR data, but they did not invert the geodetic data set. In this way, they obtained a NW-dipping fault plane imaged at depths of 5–12 km (their Figure 10), i.e., deeper than our geodetic inversion fault model (1–7 km).

Ground and road cracks/fractures were observed in the field in two clusters mainly: one located near Nipiditos village (Figure 10) and one along the highway south of Arkalochori. Note that the Nipiditos cluster is found close to the projected intersection of the seismic fault with the ground surface (see blue line in Figure 10). At one locality we observed an open fracture, oriented N-S with a length of 27 m, but no shear displacement along it (Figure 11a). It is likely that it is related to strong ground motion that caused “spreading” of the uppermost soil layers. This is also observed in the descending co-seismic interferogram (Figure S7a) where it is shown to occur also in the footwall of the Avli-Lagouta Fault. The second cluster (near Kalo Chorio in Figure 10 and Figure 11c) is probably related to triggered secondary slip [65] along a series of small faults occurring on the hanging wall of the seismic fault (Figure 11d), clearly observable on both interferograms (Figure S7b). Triggered shallow slip (TSS) traces (black dashed lines in Figure S7a are visible as discontinuities in both descending and ascending track phase interferograms (Figure S7a–c). Moreover, landslide locations were mapped in the field (west of Arkalochori) and by comparing pre- and post-event Sentinel-2 optical imagery (image dates 19 September 2021 and 29 September 2021), by tracing the changes in land cover between the two dates. Moreover, at both co-seismic interferograms (ascending & descending), we mapped several distinct surficial deformation areas (see red polygons in Figure 10) with small or incoherent displacement signals. These areas are interpreted as minor surficial displacement of triggered pre-existing large landslides, and in some locations may be linked to deep-seated gravitational displacements (DSGD) [53,66].

5.2. Normal Fault Segmentation and the 27 September 2021 Rupture

We identified the seismic fault as mostly rupturing the southern segment of the Kastelli Fault Zone (KFZ; Figure 9 and Figure 10), so it is of interest to look at the geometric relation of the earthquake rupture to the long-term (neotectonic) geological structures. To map relief due to cumulative offsets along the Kastelli Fault Zone, we used the Copernicus GLO-30 Digital Elevation Model (DEM) with 30 m resolution https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198 (accessed on 1 December 2021). Footwall topography can be used as a proxy for fault throw distribution [67,68,69] and the distribution of the fault throw along the strike may be indicative of the maturity of a fault [59,69]. The Kastelli (“Agnos-Avli”) fault Zone shows geomorphic evidence for recent activity, such as prominent range-front escarpments, V-shaped valleys, triangular facets, and axial drainage (Figure 10 and Figure 11e). In terms of footwall topography shape, which is commonly used as a criterion for fault segmentation [59,67,69], the KFZ comprises two segments, the northern (Kastelli) segment and the southern Avli-Lagouta segment (Figure 10). The Kastelli (also known as Agnos) Fault is an 11 km-long normal fault, west-dipping, with most relief (~200 m) developed towards its terminations (Figure 11b,e and Figure 12). The hanging wall topography increases systematically from North to South. The northern termination of Agnos F. may be associated with an NNS-SSW cluster of microseismicity recorded by a local network in Delibasis et al. ([39]; south of Gouves, see Figure 1). The southern Avli (Lagouta) segment is also of similar length (~12 km), but with less relief (~150 m; Figure 12). Both footwall and hanging wall topography diminish toward its southern termination. The structural boundary between the two segments is marked by the change in the direction of axial drainage on either side of the village of Agia Paraskevi (Figure 10). Given that the depth to bedrock exceeds 350 m in the Kastelli basin (I. Bouloukakis, pers. commun.) we could assume a throw of ~500 ± 100 m for the main part of the Kastelli Fault zone (Figure 12) accumulated over Late Pliocene—Quaternary times (~3 Ma). These displacement data (assuming a 50-degree dip angle of the fault plane) yield an average slip rate between 0.17–0.26 mm/yr for the Kasteli fault, thus similar to the active faults of central Greece ([59,70,71]). The estimated slip rate is close to the post-glacial rate estimated by Nicol et al. ([16]; 0.21 mm/yr) for KFZ which implies that slip rates did not change much over the Quaternary. Incidentally, this result matches the observation by Robertson et al. [15] that uplift rates in couth-central Crete were constant throughout Quaternary. Furthermore, the KFZ 0.17–0.26 mm/yr rate is in the same order with the slip rate estimate by Caputo et al. [13,29] (0.5 mm/yr) of the normal fault in the footwall of the KFZ (Geraki fault; see GFZ in Figure 2). It is reasonable to suggest that KFZ is the main NNE-SSW structure now in central Crete and accommodates most of the accumulated elastic strain. It seems that the Geraki Fault zone becomes progressively inactive with late Quaternary slip rates less than KFZ; GFZ is now located on the footwall of KFZ (Figure 1 and Figure 2) without almost any syn-rift deposits occurring along its trace [9].

The 27 September. 2021 rupture extent (Figure 9) indicates that the structural boundary along KFZ (hanging wall high; Figure 12) was breached. The co-seismic slip attained its highest value west of Nipiditos where both the footwall and hanging wall have their highest elevations. It follows that the KFZ and the NW-SE-striking Nipiditos fault (Figure 2 and Figure 9) do not intersect as suggested by [2]. The coseismic slip continued to the north of the village of Agia Paraskevi (Figure 9 and Figure 10), where long-term geomorphology shows drainage towards the north and west. The rupture caused permanent offsets of the NOA campaign GNSS point AS15 (Figure S8; see Figure 1 for location) that is located near the projected intersection of the modelled rupture with the free surface. Therefore, we suggest that the 23 km-long KFZ comprises a single fault segment that hosts overlapping ruptures. To find a pattern in the seismic behaviour of KFZ, we need more data (trenching or 36CL cosmogenic dating of fault scarps, e.g., [64,72,73]). The size of the 2021 rupture indicates the existence of a locked 10 km patch of the KFZ, towards the north (Agnos; Figure 1) which was loaded with tectonic stress following the Arkalochori event. It is well known that earthquakes with a strong or moderate moment magnitude size (Mw > 5.8) transfer static stresses upon other fault segments which are located near the main shock [52,53,74,75] thus changing recurrence intervals by modifying times to failure (advancing or delaying).

The strike-orthogonal distribution of the 2021 co-seismic displacement as seen in the cross-section (Figure 13) does not indicate a relation with long-term topography, as well. This happened because all coseismic deformation occurred on the hanging wall of Avli-Lagouta Fault (Figure 14). The displacement profile is skewed towards east, as expected since the fault dip-direction is towards the west. The surface projection of the modelled fault fits well with the Avli fault segment of the Kastelli Fault Zone, outcropping west of the village of Nipiditos (Figure 11e). The modelled ground deformation also fits well with the Plio-Pleistocene basin located east of Arkalochori (Figure 2). To match the top of the modelled fault (1.2 km beneath the surface) with the scarp of the Avli fault, in the vicinity of the village of Nipiditos, requires a dip angle of ~55° throughout the upper crust of central Crete, in agreement with the modelled fault plane (Figure 8 and Figure 9) and cancelling the need for listric geometries updip the hypocentre as implied by [2,3]. It also requires a small (~1 km) shift in location of the relocated mainshock towards east, which is an offset feature that we have seen in several earthquake sequences in Greece (e.g., [53,76]) likely caused by heterogeneities in the crustal structure not captured by the 1-D velocity model.

The Arkalochori earthquake along the Kastelli Fault Zone highlights the nearly arc-parallel (E-W to WNW-ESE) extension of the Heraklion basin in central Crete as opposed to the neighbouring E-W Messara basin, to the south (Figure 2). The Messara basin was formed due to N-S extension [10] which is active today further east [15,16,18]. However, the extension direction also switches NW-SE eastwards, on-shore Crete [3,7,8,9,13,14,16,24,26,27,29,58,64] highly oblique to the subduction direction (N5° W; [77,78]). Therefore, despite its proximity to the subduction front, Crete is deforming by extension along shallow normal faults of almost orthogonal orientations (N195° E in Arkalochori; East-west in Messara and along the south coast; [9,15,64,78]). The cause of the block-like onshore extension of Crete is the mechanical response of the brittle Aegean crust, which is quickly overriding a nearly stalled Nubian (African) subducting plate with a funnel shape [77].

5.3. The Relative Value of InSAR vs. GNSS Observations

Nowadays, InSAR technology is mature and provides reliable, spatially dense information on coseismic & postseismic displacements after shallow earthquakes [49,50,53,55,79,80,81,82,83,84,85,86,87,88]. What is evident with InSAR, when the quality is very high, as is the case over Arkalochori, is that the InSAR-based displacement of proximal campaign GPS stations (such as AS15; Figure 1) is nearly as accurate as GNSS measurements in the field. The added value of campaign GPS after events is almost zero when high-quality interferograms exists (which is not always the case, especially when the events are offshore such as the 2017 Kos earthquake; [52]). Moreover, fieldwork costs time and money. So, why do we need to make campaign GPS observations there while InSAR allows us to map ground displacements with an accuracy comparable to what we can achieve with GPS in the field?

If we want to make new contributions on earthquake rupture mechanics, what is important is to deploy GPS stations quickly after earthquakes (within the first few hours) in the near field of the epicentre area, and when possible, place instruments even in advance of peculiar seismic crises starting, as it was the case around Arkalochori since June 2021 and even earlier ([1,2,3]; see earthquake time-series in Figure 5). This will allow us to monitor the dynamics of the ground deformation immediately after earthquakes (and sometimes before), capture the postseismic deformation earlier and better understand how the stress/strain migrates during aftershocks sequences, especially when large aftershocks occur in the days/weeks following the first main shock.

Certainly, GPS observations are useful or needed before and after earthquakes in the near field of earthquakes, for example to assess the secular deformations (see [23,52,76,79]) but when they are outperformed by InSAR for the assessment of the co-seismic displacement, the motivation for performing those measurements must be established. We believe there is a real value in the assessment of the secular deformations; therefore, the selection of the GPS points, their direct connection to the bedrock, their long-term stability and preservation over decades, if not centuries, should be central in the design of the benchmarks and their maintenance. This need makes it necessary to apply very strict rules for the design, observation and maintenance of geodetic networks tailored for the long and very long term. This approach could be central to the way of thinking inside Geophysical Observatories in particular because both short-term and long-term observation is their mission.

6. Conclusions

• The 27 September 2021 M5.9 earthquake produced significant ground deformation in the vicinity of the village of Arkalochori, but no surface ruptures of primary (tectonic) origin.
• The distribution of the relocated aftershocks revealed a west-dipping fault with activity limited to depth up to 12 km. The clustering of aftershocks indicates that another shallow structure or a neighbouring segment of the same fault zone was activated right after the occurrence of the mainshock, to the NE of the rupture.
• We used geodetic data (InSAR & GNSS) to derive a finite fault model. The along-strike averaged coseismic slip has a maximum at depth of 4–6 km, with a maximum amplitude of 1.2 m.
• The geodetic moment of this earthquake was found in the order of 1.1 × 10¹⁸ N·m, which corresponds to an earthquake moment magnitude Mw = 5.96.
• Our data (seismological and geodetic) support the activation of a 13 km-long planar fault striking N195° E and dipping 55° towards WNW which intersects the free surface along the mountain front of the Kastelli basin.
• The seismic fault may be identified with the neotectonic Avli Fault segment (AF) of the NNE-SSW striking, west-dipping, 23 km-long Kastelli Fault zone (KFZ). Part of the rupture propagated along the Kastelli Fault segment (KF) indicating a possible fault segment linkage and, therefore, a history of overlapping ruptures along the KFZ.
• Based on geological data and footwall topography, we estimate an average Quaternary slip rate between 0.17–0.26 mm/yr for the KFZ.
• We found no tectonic surface breaks. The decrease in the fault offsets toward the Earth’s surface is likely caused by an increased frictional resistance of the shallow layer to rapid coseismic slip.
• Satellite observations made in the first month after the earthquake detected no post-seismic deformation (i.e., below one fringe or 2.8 cm) but the GNSS station ARKL recorded a 2 cm postseismic motion with a relaxation period of 75 days.