On the Patterns and Scaling Properties of the 2021–2022 Arkalochori Earthquake Sequence (Central Crete, Greece) Based on Seismological, Geophysical and Satellite Observations

: The 27 September 2021 damaging mainshock (Mw6.0) close to Arkalochori village is the strongest earthquake that was recorded during the instrumental period of seismicity in Central Crete (Greece). The mainshock was preceded by a signiﬁcant number of foreshocks that lasted nearly four months. Maximum ground subsidence of about 18 cm was estimated from InSAR processing. The aftershock sequence is located in an almost NE-SW direction and divided into two main clusters, the southern and the northern ones. The foreshock activity, the deformation area, and the strongest aftershocks are located within the southern cluster. Based on body-wave travel times, a 3-D velocity model was developed, while using combined space and ground-based geodetic techniques, the co-seismic ground deformation is presented. Moreover, we examined the co-seismic static stress changes with respect to the aftershocks’ spatial distribution during the major events of the foreshocks, the Mw = 6.0 main event as well as the largest aftershock. Both the foreshock and the aftershock sequences obey the scaling law for the frequency-magnitude distribution as derived from the framework of non-extensive statistical physics (NESP). The aftershock production rate decays according to the modiﬁed Omori scaling law, exhibiting various Omori regimes due to the generation of secondary aftershock sequences. The analysis of the inter-event time distribution, based on NESP, further indicates asymptotic power-law scaling and long-range correlations among the events. The spatiotemporal evolution of the aftershock sequence indicates triggering by co-seismic stress transfer, while its slow migration towards the outer edges of the area of the aftershocks, related to the logarithm of time, further indicates a possible afterslip.


Introduction
Greece is located at the southeastern tip of Europe where a variety of geological processes take place, such as the formation of the Alpine Mountain chain from the Western French Alps to the Dinarides in the Balkan Peninsula due to the collision between European and Nubian plates.Furthermore, observed intense deformation in the Aegean and the surrounding regions produces a significant portion of SE Europe's seismicity, concentrated in thrust structures along the Hellenic Arc and smaller extensive ruptures in the area of the Aegean back-arc [1][2][3][4].
The epicentral area is located in the front of the Hellenic arc, near the region where the European and Nubia plates converge, giving rise to large earthquakes [1][2][3][4][5].The most characteristic tectonic features in the vicinity of the rapture zone are the two large tectonic basins located in the northern Heraklion and the southwest Messara.These were formed by extensional forces from an arc-normal pull since 11 Ma which formed the average NNE-SSW direction faults [1,2].The basins are filled with Miocene to Pliocene sediments overlaid by Quaternary deposits and in the north-eastern part, there are exposed nappes tectonic features consisting of Carbonates of the Tripolis and Trypalion units.
Destructive earthquakes occurred in the broader Crete region during the pre-instrumental era [16][17][18].In Central Crete, two major historical earthquakes have been reported, both in the vicinity of Heraklion.The first one was on 1 July 1494 (M = 5.4) [16] while on 26 November 1595, a second event (M ~6.4) took place, both causing severe damage on the island of Crete [16,18].
On 27 September 2021 an Mw = 6.0 event, linked to an approximately N-S trending normal fault at the central part of the island of Crete, Greece, occurred ~20 km to the south of the prefecture's capital.The occurrence of the mainshock took place after a long time of foreshocks since the first half of June 2021 with criticality characteristics [19] and a long aftershock sequence with the strongest event occurring on 28 September with local magnitude M L 5.3.Following the recent report by ITSAK [20], the recorded PGA at the epicentral area (Arkalochori) was 0.62 g in the horizontal component (N-S) and 0.82 g in the vertical one, with a duration of strong ground motion (>0.1 g) almost 6 s.Its focal mechanism is characterized by an SSW-NNE to SW-NE-trending, nearly dip-slip normal faulting.Its strike generally ranges from N200 • E-N230 • E and its dip angle varies between 40 • and 60 • .The active fault associated with the main event is the Kastelli Fault, which has a progressive change in the strike from 225 • to 265 • northeastwards and dips between 60 • -80 • northwestwards [9,[21][22][23].
In this work, we present the consistency between seismological, geodetic, satellite, and geophysical data for the 2021-2022 Arkalochori earthquake sequence, highlighting the complementarity of multi-disciplinary approaches.First, we relocate the 2021-2022 Central Crete earthquake sequence to have a more accurate insight on its distribution and scaling properties.Based on body-wave travel times, we have developed a 3D velocity model, while using combined space and ground based geodetic techniques the co-seismic ground deformation is obtained.Moreover, we examined the co-seismic static stress changes with respect to the aftershocks' spatial distribution during the major events of the foreshocks, the Mw = 6.0 main event as well as the largest aftershock.In addition, the scaling properties of the 2021-2022 Central Crete earthquake sequence are investigated thoroughly and we analyzed the aftershocks decay rate based on the modified-Omori formula.Furthermore, the frequency-magnitude distribution parameters of the earthquake sequence along with that of the distribution of inter-event times between the successive events are viewed in terms of non-extensive statistical physics in order to provide more detailed insights into the complex nature of the long-term correlation of the earthquake generation process.

Seismological Data and Earthquake Sequence Analysis
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 since the occurrence of the Mw = 6.0 main event on 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 Research Center (HC network at formerly Technological Appl.Sci.2022, 12, 7716 3 of 32 Educational Institute of Crete [24]), at a distance of about 12-22 km.These data were manually revised using a regional 1D model and were put into a relocation procedure in order to obtain more accurate results, as will be discussed in the next sections.

Data Analysis
The Arkalochori earthquake sequence is divided into two main temporal groups, one that preceded the 27 September Mw = 6.0 mainshock, consisting of 620 events with a significant rise in numbers during July and August 2021, and the aftershock sequence, divided into three sub-groups, composed by 4130 seismic events (Figure 1).A major part of the sequence was recorded by local stations of the regional Hellenic Unified Seismological Network (HUSN), with the nearest stations being KNSS, PFKS (IFEGG), located ~20 km to the SW and NE of the epicenter, respectively.On October 1, 2021, the Geodynamics Institute of the National Observatory of Athens (GI-NOA) installed 4 temporary stations (CRE1-4) around the aftershock zone, contributing to the depth accuracy and providing a preliminary hypocentral solution for this time period.earthquake swarm in early June 2021.The situation, however, changed drama since the occurrence of the Mw = 6.0 main event on 27 September 2021.Although permanent stations of the Hellenic Unified Seismological Network operate on the of Crete at the time of the mainshock, the closest station to the epicentral area was of the Hellenic Mediterranean University Research Center (HC network at fo Technological Educational Institute of Crete [24]), at a distance of about 12-22 km data were manually revised using a regional 1D model and were put into a relo procedure in order to obtain more accurate results, as will be discussed in the ne tions.

Data Analysis
The Arkalochori earthquake sequence is divided into two main temporal g one that preceded the 27 September Mw = 6.0 mainshock, consisting of 620 events significant rise in numbers during July and August 2021, and the aftershock seq divided into three sub-groups, composed by 4130 seismic events (Figure 1).A maj of the sequence was recorded by local stations of the regional Hellenic Unified S logical Network (HUSN), with the nearest stations being KNSS, PFKS (IFEGG), l ~20 km to the SW and NE of the epicenter, respectively.On October 1, 2021, the G namics Institute of the National Observatory of Athens (GI-NOA) installed 4 tem stations (CRE1-4) around the aftershock zone, contributing to the depth accura providing a preliminary hypocentral solution for this time period.A total of 4750 events of the 2021-2022 Arkalochori sequence that occurred the period between 13 January and 31 January 2022 (Figure 2) were detected and ally analyzed using real-time waveform data from the Hellenic Unified Seismo Network (HUSN) in the SeisComP3 graphical user interface [25].During the first s the sequence analysis, hypocenters were located in near real-time, by employi Hypo71 single-event algorithm [26] and a custom regional 1D velocity model Hellenic peninsula [27].A total of 4750 events of the 2021-2022 Arkalochori sequence that occurred during the period between 13 January and 31 January 2022 (Figure 2) were detected and manually analyzed using real-time waveform data from the Hellenic Unified Seismological Network (HUSN) in the SeisComP3 graphical user interface [25].During the first stage of the sequence analysis, hypocenters were located in near real-time, by employing the Hypo71 single-event algorithm [26] and a custom regional 1D velocity model for the Hellenic peninsula [27].
In this study, two local 1D velocity models [28,29] have been used through the second stage of the data processing, running the HypoInverse code [30].Residual values from these models were compared (Table 1) with no significant differences, while the epicentral differences were less than 0.5 km, whereas the depths, which are more sensitive to the velocity model, differed by about 1 km on average.The [28] velocity model (Model 1) provided much shallower events than the respective ones of [29] (Model 2) and especially for the stronger events of the sequence (Mw = 6.0 and Mw = 5.3) that were located at a depth shallower than 5 km (2.7 and 0.7 km, respectively), which seemed unrealistic in terms of earthquake physics and the past of the area.In this study, two local 1D velocity models [28,29] have been used through the second stage of the data processing, running the HypoInverse code [30].Residual values from these models were compared (Table 1) with no significant differences, while the epicentral differences were less than 0.5 km, whereas the depths, which are more sensitive to the velocity model, differed by about 1 km on average.The [28] velocity model (Model 1) provided much shallower events than the respective ones of [29] (Model 2) and especially for the stronger events of the sequence (Mw = 6.0 and Mw = 5.3) that were located at a depth shallower than 5 km (2.7 and 0.7 km, respectively), which seemed unrealistic in terms of earthquake physics and the past of the area.The final hypocentral locations were obtained using the velocity model of [29].The range of depth distribution is mainly between 5 and 15 km for the aftershocks of the 27 September Mw = 6.0 event.Figure 2 presents the distribution of 4750 events that were manually revised and relocated with HypoInverse code, along with the faults in the area as extracted from [31].
The seismic sequence was divided into five sub-groups according to its spatiotemporal occurrence (Figure 3):  The final hypocentral locations were obtained using the velocity model of [29].The range of depth distribution is mainly between 5 and 15 km for the aftershocks of the 27 September Mw = 6.0 event.Figure 2 presents the distribution of 4750 events that were manually revised and relocated with HypoInverse code, along with the faults in the area as extracted from [31].
The seismic sequence was divided into five sub-groups according to its spatiotemporal occurrence (Figure 3): 1.
27 September-28 September 2021 (period B), the first day of the aftershock sequence and just a few hours before the greatest aftershock (M5.3), composed of 90 events; 3.
12 October-31 October 2021 (period D), where the M = 4.0 event took place after a significant decay in the aftershocks number in Arkalochori; 5.
1 November-31 January 2022 (period E), which comprises small to moderate magnitude events in a deeper part of the crust (H > 20 km) located near Herakleion.
3. 28 September-12 October 2021 (period C), just after the occurrence of the M = 5.3 aftershock at 04:48 UTC, consisting of 803 events; 4. 12 October -31 October 2021 (period D), where the M = 4.0 event took place after a significant decay in the aftershocks number in Arkalochori; 5. 1 November 1-31 January 2022 (period E), which comprises small to moderate magnitude events in a deeper part of the crust (H > 20 km) located near Herakleion.

Hypocentral Relocation of the Earthquake Sequence
Accurate earthquake hypocenter parameters are required in order to obtain a detailed image of the structural properties and processes that trigger seismic activity.The precision of hypocenter locations and their uncertainties depend on several factors, including the number and quality of available seismic phases, the accuracy with which arrival times are measured, the network geometry, the knowledge of the velocity structure, and the linear approximation to a set of non-linear equations, which is assumed in the inversion.hypoDD [32] is an algorithm that reduces residuals between observed and theoretical differences of travel times (or double differences) for pairs of neighboring events at each station that recorded both events, as can be seen from Equation (1).This way, errors due to unmodeled velocity structures are minimized without station corrections.A minimum 1D layered velocity model is used to predict the travel time differences and partial derivatives (Equation (2)).
where m = (x, y, z) Inter-event distance and misfit weighting are applied to catalog data after the end of each iteration, in order to optimize their quality during the relocation procedure.Hori-

Hypocentral Relocation of the Earthquake Sequence
Accurate earthquake hypocenter parameters are required in order to obtain a detailed image of the structural properties and processes that trigger seismic activity.The precision of hypocenter locations and their uncertainties depend on several factors, including the number and quality of available seismic phases, the accuracy with which arrival times are measured, the network geometry, the knowledge of the velocity structure, and the linear approximation to a set of non-linear equations, which is assumed in the inversion.hypoDD [32] is an algorithm that reduces residuals between observed and theoretical differences of travel times (or double differences) for pairs of neighboring events at each station that recorded both events, as can be seen from Equation (1).This way, errors due to unmodeled velocity structures are minimized without station corrections.A minimum 1D layered velocity model is used to predict the travel time differences and partial derivatives (Equation ( 2)).
where m = (x, y, z).Inter-event distance and misfit weighting are applied to catalog data after the end of each iteration, in order to optimize their quality during the relocation procedure.Horizontal and vertical relative spatial errors can be minimized by approximately one order of magnitude under certain conditions.In this study, more than 4500 events of M L ≥ 0.6, comprising 58,777 phases, were relocated with hypoDD.Among the main factors that had been taken into consideration during the relocation procedure were the following: (a) network coverage of the area, (b) the size of the dominant clusters, and (c) their maximum separation distance.This led to the formation of 98,070 P-and 43,812 S-phase pairs, respectively, for the whole volume, while in the central cluster (cluster #1), where the mainshock (M = 6.0) largest aftershock (Mw = 5.3) is situated, 98,070 P-and 18,673 S-phase pairs were formed.
In the area of Central Crete, 4728 out of 4750 events of the initial catalog (M ≥ 0.6), were relocated with hypoDD, giving a first result that could be rated as satisfactory.The mean temporal errors (rms) were reduced from 0.14 s to 0.11 s, while the spatial errors (erx, ery, erz) were decreased from 1.0, 0.8, and 2.0 km to 0.3, 0.3, and 0.4 km, respectively (Table 1).The hypoDD-estimated errors in the final locations were calculated using the LSQR method, which may not be representative of the real ones [32].
The epicenter of the mainshock was located less than 3 km to the SE of Arkalochori (lat: 35.1416 • N, long: 25.2736 • E) at a depth of ~9.6 km, obtained by the double difference algorithm procedure.The optimization of the final results leads to the clustering of the earthquake sequence into four main clusters.A dense cluster of events has occurred west of the mainshock, in an approximately 15 km long area associated with the foreshocks (cluster 1).The epicenters of cluster 2 were mainly distributed in the area between Amourgeles and Parthenio N-S oriented normal faults, in the region to the west of Arkalochori.Further to the NE, another significant cluster of events (cluster 3) was also observed, in the footwall of Agnos NE-SW striking normal fault, near Kastelli (Figure 4).Most seismicity is in a range of focal depths between 7 and 18 km.maximum separation distance.This led to the formation of 98,070 P-and 43,812 S-phase pairs, respectively, for the whole volume, while in the central cluster (cluster #1), where the mainshock (M = 6.0) largest aftershock (Mw = 5.3) is situated, 98,070 P-and 18,673 S-phase pairs were formed.
In the area of Central Crete, 4728 out of 4750 events of the initial catalog (M ≥ 0.6), were relocated with hypoDD, giving a first result that could be rated as satisfactory.The mean temporal errors (rms) were reduced from 0.14 s to 0.11 s, while the spatial errors (erx, ery, erz) were decreased from 1.0, 0.8, and 2.0 km to 0.3, 0.3, and 0.4 km, respectively (Table 1).The hypoDD-estimated errors in the final locations were calculated using the LSQR method, which may not be representative of the real ones [32].
The epicenter of the mainshock was located less than 3 km to the SE of Arkalochori (lat: 35.1416°N, long: 25.2736° E) at a depth of ~9.6 km, obtained by the double difference algorithm procedure.The optimization of the final results leads to the clustering of the earthquake sequence into four main clusters.A dense cluster of events has occurred west of the mainshock, in an approximately 15 km long area associated with the foreshocks (cluster 1).The epicenters of cluster 2 were mainly distributed in the area between Amourgeles and Parthenio N-S oriented normal faults, in the region to the west of Arkalochori.Further to the NE, another significant cluster of events (cluster 3) was also observed, in the footwall of Agnos NE-SW striking normal fault, near Kastelli (Figure 4).Most seismicity is in a range of focal depths between 7 and 18 km.Furthermore, five (5) cross-sections were performed in order to see the impact of the relocation procedure on the sequence hypocentral depths and the discrimination of the local activated structures.Cross-sections 1-2 have an NNE-SSW orientation, and 3-5 WNW-ESE direction, perpendicular to the NE trending faults.The geometry of the hypocenters as they appear in the performed cross-sections, reveal the activation of a fault, dipping ~60 • to the WNW, and a smaller antithetic structure, possibly connected to Galatas N-S striking normal fault (Figure 5).
An almost sub-vertical structure makes an appearance in the cross-sections north of the epicenter of the mainshock (sections B-B', D-D'; Figure 5).The earthquake activity that belongs to cluster 4, started on 16 January 2022, with an event of M L = 3.3 and continued throughout the first two months of 2022 as the latest contribution to the seismic sequence.
The hypocenters appear to be located in deeper parts of the crust and they have an apparent dip towards the WNW.
Furthermore, five (5) cross-sections were performed in order to see the impact o relocation procedure on the sequence hypocentral depths and the discrimination o local activated structures.Cross-sections 1-2 have an NNE-SSW orientation, and WNW-ESE direction, perpendicular to the NE trending faults.The geometry of th pocenters as they appear in the performed cross-sections, reveal the activation of a dipping ~60° to the WNW, and a smaller antithetic structure, possibly connected t latas N-S striking normal fault (Figure 5).A-A', B-B', C-C', D-D  E-E') on the left side of the figure and the results of the SSW-NNE oriented cross-sections (  A-A', B-B', C-C', D-D', and  E-E

Travel Time Tomography
Local Earthquake Tomography (LET) techniques have been successfully applied to reveal the velocity structure in such cases of aftershock sequences.In this paper, the bodywave inversion was based on the LOcal TOmographic Software (LOTOS) by [33].P-and S-phases of more than 800 earthquakes, recorded during the 2021-2022 time period by local and regional stations of the Hellenic Unified Seismological Network (HUSN) and the Hellenic Strong Motion Network (HSMN), located in Southern Greece, were used for the tomographic inversion.Checkerboard tests were performed to set the input parameter values that produced better resolution and increased the fidelity area.Regarding the 3D tomographic inversion, a dataset consisting of 12,236 P-and 9820 S arrival times was chosen, with at least 12 phases per event (Supplementary Material Figure S2).LOTOS code provides two alternative options: inversion for V P and Vs (V P -V S scheme) using P and S travel-time residuals (dt P and dt S ) and inversion for V P and V P /V S ratio (V P -V P /V S scheme) using dt P and differential residuals, dt S -dt P .In this work, inversion was performed for V P -V S and V P -V P /V S schemes, in order to obtain additional constraints concerning the V P and Vs anomalies [33,34].
The stations' coordinates, their elevation, and the body-wave arrival times are essential as input data to the algorithm.The hypocenter locations and the origin times are not necessarily needed, given that their determination is performed during the execution of the calculations.However, if preliminary hypocentral locations are available, they are used to decrease the processing time.Moreover, the available initial 1D velocity model [29] and a set of input parameters, i.e., parameterization, grid dimensions, and damping parameters, are defined by the user [33].A nodal representation was employed, given that the velocity field, reconstructed by a three-dimensional grid, does not assume a specific geometry of heterogeneities [35].The grid spacing (~2 km) was kept considerably smaller than the expected resolution length, to reduce the bias of the resulting models due to the grid configuration.The optimal grid mesh has been determined considering the stations/events geometry.In addition, to further decrease the influence of the model parameterization on the solutions, the inversion was repeated using four grid orientations (0 • , 22 • , 45 • , and 67 • ).The inversion results, obtained for the previously mentioned grids, were stacked into one summary model, reducing the artifacts related to grid orientation, as described by [33].
The values of the P-and S-wave residuals during different iteration steps of the inversion procedure are presented in Table 2.For the P-and S-data, the reduction of the residual is ~10% and 13%, respectively.The resulting P-and S-wave velocity anomalies with respect to the starting 1D velocity model are shown in horizontal depth slices (Figures 5 and 6).The interpretation of the obtained results is limited to the unmasked confidence regions, given that they are characterized by reasonable reconstruction of the checkerboard model.The mean computed P and S anomalies for the study area do not exceed ±13%.Strong NW-SE and E-W oriented negative velocity anomalies predominate at both the upper and the lower crust of Central Crete.These are observed down to 15 km depth in the tomograms of Figures 6 and 7. Deeper than 20 km the model lacks resolution and can only be considered indicative, hence it is not discussed.At the depth slice of 5 km, a NE-SW-trending zone of negative body-wave velocity perturbations appears near the epicentral region of the Mw = 6.0 Arkalochori earthquake (Figures 6 and 7).This anomaly follows the mean distribution of Alluvial deposits and post-alpine sediments which are bounded by positive (~13%) body-wave velocity perturbations, possibly connected to the older post-alpine sediments of Viannos formation and the Mesozoic carbonate rocks to the east and south of Arkalochori basin, respectively [13,22,36,37].In the area north of Figure 6.Lateral V P (%) perturbations at 5, 10, 15, and 20 km depths.Areas with lower resolution are masked (darkened).Fault traces derived by [31].
Strong NW-SE and E-W oriented negative velocity anomalies predominate at both the upper and the lower crust of Central Crete.These are observed down to 15 km depth in the tomograms of Figures 6 and 7. Deeper than 20 km the model lacks resolution and can only be considered indicative, hence it is not discussed.At the depth slice of 5 km, a NE-SWtrending zone of negative body-wave velocity perturbations appears near the epicentral region of the Mw = 6.0 Arkalochori earthquake (Figures 6 and 7).This anomaly follows the mean distribution of Alluvial deposits and post-alpine sediments which are bounded by positive (~13%) body-wave velocity perturbations, possibly connected to the older postalpine sediments of Viannos formation and the Mesozoic carbonate rocks to the east and south of Arkalochori basin, respectively [13,22,36,37].In the area north of Arkalochori, an E-W-trending anticorrelated pattern of negative P-and positive S-wave velocity anomalies are observed at the depth range of In the depth slice of 15 and 20 km, and almost NNE-SSW discontinuity of positive to the west and negative to the east V P anomalies (Figure 6) are identified along the west-dipping Agnos normal fault.Furthermore, cross-sections B-B' and C-C' in both primary (P) and secondary (S) wave velocity anomalies (Figures 8 and 9), reveal this west-dipping structure that may be related to Agnos high-angle (~60 • ) normal fault [13,[21][22][23]36,37].
fault direction, from the area east of Arkalochori towards the town of Malia (35.20°N-35.27°N,25.34°E-25.45°E).In the depth slice of 15 and 20 km, and almost NNE-SSW discontinuity of positive to the west and negative to the east VP anomalies (Figure 6) are identified along the west-dipping Agnos normal fault.Furthermore, cross-sections B-B' and C-C' in both primary (P) and secondary (S) wave velocity anomalies (Figures 8 and  9), reveal this west-dipping structure that may be related to Agnos high-angle (~60°) normal fault [13,[21][22][23]36,37].

Co-Seismic Ground Deformation
To study the co-seismic ground deformation on the epicentral area of the 27 September 2021 event, we combined space and ground-based geodetic techniques.Radar satellite data were used to produce differential interferograms to spatially study the co-seismic deformation.Moreover, geodetic data from continuous Global Navigation

Co-Seismic Ground Deformation
To study the co-seismic ground deformation on the epicentral area of the 27 September 2021 event, we combined space and ground-based geodetic techniques.Radar satellite data were used to produce differential interferograms to spatially study the co-seismic deformation.Moreover, geodetic data from continuous Global Navigation Satellite System (GNSS) stations operating on the broad affected area provided point-wise accurate 3D displacement vectors.

Interferometric Data and Results
The term SAR stands for Synthetic Aperture Radar [38].Differential Interferometry (Differential InSAR technique-DInSAR) is an advanced technique [39] aimed at detecting surface movements due to geophysical phenomena or human interventions.Since the 1990s, the DInSAR technique has proven to be an interesting tool for measuring and observing ground deformation suitable for analyzing geodynamic processes, e.g., [40][41][42][43][44].
To map the co-seismic ground deformation due to the 27 September 2021 mainshock, we used two SAR image pairs, one on ascending and one on descending orbital geometries (Table 3), acquired from ESA's Sentinel-1A and Sentinel-1B satellites (https://scihub.copernicus.eu/,accessed on 31 May 2022).In both cases, the reference image was the one before the earthquake occurrence, while the repeat image was the one that refers to the date after the event.Each reference-repeat pair was processed using the ESA's SNAP software and two individual interferograms were generated.The topographic phase was subtracted using the SRTM 1 arc-second Digital Terrain Model, a 30-m resolution Shuttle Radar Topography Mission Digital Elevation Model (USGS 1 ARC-second SRTM DEM, https://doi.org/10.5066/F7DF6PQS,accessed on 31 May 2022), while the signal to noise ratio was enhanced by applying the adaptive power spectrum filter of [45] with a coherence threshold of 0.4.The two-phase wrapped interferograms were then used to estimate co-seismic ground deformation.To calculate the terrain displacement, an unwrapping process was performed, and the phase unit was transformed into distance units in the satellite line of sight (LoS) for each interferometric pair.Finally, decomposition of the ascending and descending LoS displacement vectors was performed to extract the vertical (up-down) and horizontal (east-west) ground deformation components.
The two wrapped interferograms are of good quality and contain the phase difference, between reference and repeat images, produced by the main seismic event and its aftershocks until the 29th of September.Due to the low temporal geometric baselines (six days), there are no areas of low coherence in the interferograms.Six fringes, forming a lobe, are evident both in ascending and descending wrapped phase interferograms (Figure 10, upper part).This asymmetrical displacement pattern is characteristic of normal-faulting earthquakes indicating that the subsidence is larger than the uplift.Each interference fringe is a phase change that corresponds to a motion of 2.8 cm in the satellite line of sight.
The LoS displacement map in the ascending geometry shows negative LoS displacement values up to 18 cm, after the Mw = 6.0 earthquake including the ground deformation caused by all the seismic events that occurred in the time interval from 23 September to 29 September, (Figure 10, lower part), while the LoS displacement map in the descending geometry shows negative LoS values up to 20 cm after the Mw = 6.0 earthquake, including all the aftershocks, occurred until the 1st of October.On the descending orbital geometry, the maximum value of ground deformation has shifted east of the epicenter of the main earthquake (Figure 10) with respect to the ascending displacement map.After the decomposition of ascending and descending LoS displacement maps, the ground deformation in vertical (up-down) and east-west directions were extracted (Figure 11a).Subsidence up to 18 cm has been calculated from the displacement decomposition in the vertical (up-down) direction while no uplift was detected.The horizontal (east-west) displacement map reveals an eastward motion up to 6 cm in the area west of the epicenter and a westward movement up to 7.6 cm for the area east of Arckalochori (Figure 11b).The latter is in satisfactory agreement with [21] where the displacement is associated with the horizontal motion of the strike-slip component of the mean event focal mechanism.The LoS displacement map in the ascending geometry shows negative LoS displacement values up to 18 cm, after the Mw = 6.0 earthquake including the ground deformation caused by all the seismic events that occurred in the time interval from 23 September to 29 September, (Figure 10, lower part), while the LoS displacement map in the descending geometry shows negative LoS values up to 20 cm after the Mw = 6.0 earthquake, including all the aftershocks, occurred until the 1st of October.On the descending orbital geometry, the maximum value of ground deformation has shifted east of the epicenter of the main earthquake (Figure 10) with respect to the ascending displacement map.After the decomposition of ascending and descending LoS displacement maps, the ground deformation in vertical (up-down) and east-west directions were extracted (Figure 11a).Subsidence up to 18 cm has been calculated from the displacement decomposition in the vertical (up-down) direction while no uplift was detected.The horizontal (east-west) displacement map reveals an eastward motion up to 6 cm in the area west of the epicenter and a westward movement up to 7.6 cm for the area east of Arckalochori (Figure 11b).The latter is in satisfactory agreement with [21] where the displacement is associated with the horizontal motion of the strike-slip component of the mean event focal mechanism.

GNSS Data and Results
GNSS data from four continuous stations (belonging to the private network of MET-RICA SA) on the central-eastern part of Crete were used in the current study to measure the co-seismic displacement.One of the stations (ARKL) is located almost above the hypocenter of the main event, in Arkalochori.The other sites are located in the city of Heraclion (HERA station; ~24 km NNW of Arkalochori), in the village of Moires (MOI1 station; ~38 km WSW from epicenter), and in the Ierapetra region (IERA station; ~49 km ESE from the epicenter).Daily raw GNSS data from these four stations were processed for a period of several years before the 2021 seismic sequence up to 30 April 2022, using the Bernese v5.2 GNSS software [46].
In the processing procedure of the local GNSS data, several stations of the EUREF and IGS were included together with other auxiliary files.The absolute antenna phase center corrections were used, together with precise orbital solutions from the Center for Orbit Determination in Europe (CODE) and the Vienna Mapping Functions for the tropospheric modeling.For the coordinate estimation on static mode solutions, the precise double difference method was used.For ambiguity resolution, numerous strategies were applied, based on the length of the formed baselines between the GNSS stations.The processing resulted in the estimation of high-precision station coordinates.Time series were formed, annual velocities prior to the 2021 seismic event were calculated and coseismic displacements were determined.The daily coordinates of the GNSS stations were estimated on the global ITRF2014 reference frame.
Based on the formed time series of the stations' coordinates (Figures S5 and S6), for the pre-seismic period, all the stations show SE horizontal linear type of motion (with respect to ITRF 2014) and subsiding vertical motion, consistent with the velocity field of the area [47].Calculating the change in the baseline distance between the stations prior to the seismic sequence it is observed that a small extension pattern occurred in the area, since lengthening on the baselines was recorded of small amplitude (~1.5 mm/year).
The main seismic event on Arkarochori caused, as is expected, strong ground displacement in the epicentral station ARKL as well as in the HERA station (Table 4).The vector of the co-seismic displacement was calculated as the static shift of the station coordinates nine days prior to and three days after the 27 September 2021 event, to minimize the effect of possible post-seismic motions.The higher displacement was recorded on the ARKL station, while noticeable displacement occurred in the HERA site (Figure S5).The two other sites (MOI1 and IERA) have not shown any considerable co-seismic motion (Figure S6).The displacement vector on ARKL site shows a strong subsiding component of ~16 cm and significant eastward and northward horizontal motion.Both the E-W and the vertical components deduced by the GNSS analysis agree with the interferometric results.Small discrepancies in the amplitude of the components deduced by the two techniques may be attributed to the incompetence of the DInSAR to define the north motional component that is inherent in the LoS vector, and in our case is quite noticeable (~8 cm).The displacement vector on HERA station shows NW horizontal co-seismic motion and a noticeable upward vertical component.
The overall image of the ground deformation in the epicentral area, based on both interferometric and GNSS results, defines a strong subsiding pattern with substantial horizontal motional component compatible with normal faulting activated structure.
In the period following the mainshock and up to 30 April 2022, the two GNSS sites in Arkalochori and Heraklion city show increased velocity values, compared to the previous period, indicating that the post-seismic relaxation is continuous up to this date.

Spatial Footprint of Coulomb Stress Changes
Numerous studies of strong earthquakes show a correlation between the positive lobe of Coulomb stress changes and the locations of the majority of the most significant aftershocks, (e.g., [48,49]).A moderate earthquake of magnitude Mw = 4.9 occurred on 24 July 2021, 02:07:37 (UTC) accompanied by many aftershocks, before the main and catastrophic earthquake of magnitude Mw = 6.0 and its major aftershock of magnitude Mw = 5.3, which took place on 27 September 2021, 06:17:21 (UTC) and 28 September 2021, 04:48:08 (UTC), respectively.Herein, we examine the co-seismic static stress changes with respect to the aftershock spatial distribution during the event of Mw = 4.9, the Mw = 6.0 main event as well as the Mw = 5.3 major aftershock.The ∆CFS changes were determined via Coulomb3.3software [50] in an elastic half-space and a uniform slip on the rupture planar surfaces.
The Coulomb Failure Stress changes (∆CFS) are given by the Equation ( 3): where ∆τ and ∆σ are the in-shear stress and normal stress, respectively, while the µ f is the effective friction coefficient [51][52][53].For the shear modulus and Poisson's ratio, we used the values of 3.3 MPa and 0.25, respectively, and a mean value for the coefficient of friction equal to µ f = 0.4 [54].
To calculate the subsurface fault's width and length, we used the empirical relations of [55] for each modeled earthquake.In addition, we used the available focal mechanism solutions from various agencies, and we present our preferred models for this study in Table 5. Figure 12 shows the determined co-seismic ∆CFS changes caused by the event of magnitude Mw = 4.9 at a centroid depth of 8.0 km and the vertical cross-sections AB and CD.The spatial distribution of the ∆CFS reveals a stress decrease towards NW and SE and a stress increase towards NE and SW of the ruptured fault.The aftershocks during this period and before the main event of Mw = 6.0 are distributed mainly along and under the fault up to 20 km depth.The ∆CFS values at the hypocenters of the Mw = 6.0 and Mw = 5.4 events were calculated, and the results provide −0.104 MPa and +0.0169 Mpa, respectively.The co-seismic Coulomb stress variations caused by the strong earthquake of Mw = 6.0 and its major aftershock (Mw = 5.4) at centroid depths of 10.0 km and 9.0 km, respectively, as well as the spatial distribution of aftershocks, are presented in Figures 13 and 14.For both seismic events, the same pattern on the spatial distribution of the ∆CFS is observed, which indicates a stress decrease towards NW and SE and a stress increase towards NE and SW of the ruptured faults.3).(Down) Coulomb stress changes along the vertical cross-section AB.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 4.9 and before the major earthquake Mw = 6.0.The green lines show the surface projections of the two fault models.3).(Down) Coulomb stress changes along the vertical cross-section AB.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 4.9 and before the major earthquake Mw = 6.0.The green lines show the surface projections of the two fault models.3).(Right) Coulomb stress changes along the vertical cross-sections A-B, C-D, E-F, and the parallel cross-section G-H (from up to down).The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 6.0 main shock.The green lines show the surface projections of the two fault models.
The focal mechanism of the mainshock was related to normal faulting and only a small portion of the strike-slip component was involved.A similar solution is evident for the largest foreshock (Mw4.9) and aftershock (Mw5.3)(Table 5) although, in the latter, the strike-slip component is increased.
From the obtained co-seismic ∆CFS changes we thus observe that most aftershocks, including those of greater magnitude, occurred within positive static stress changes produced by the major earthquake, and by the strongest aftershock.This suggests that the spatial distribution of aftershocks, including the significant ones, are controlled by the co-seismic Coulomb stress changes produced during the Mw = 6.0 mainshock and the major events of the sequence.The focal mechanism of the mainshock was related to normal faulting and only a small portion of the strike-slip component was involved.A similar solution is evident for the largest foreshock (Mw4.9) and aftershock (Mw5.3)(Table 5) although, in the latter, the strike-slip component is increased.
From the obtained co-seismic ΔCFS changes we thus observe that most aftershocks, including those of greater magnitude, occurred within positive static stress changes produced by the major earthquake, and by the strongest aftershock.This suggests that the spatial distribution of aftershocks, including the significant ones, are controlled by the co-seismic Coulomb stress changes produced during the Mw = 6.0 mainshock and the major events of the sequence.

Frequency-Magnitude Scaling Properties of the Foreshock and Aftershock Sequences in Terms of Non-Extensive Statistical Physics
The frequency-magnitude distribution (FMD), known as G-R law [56], is of vital importance for the characterization of a seismic sequence [56][57][58][59][60][61][62][63][64][65][66][67][68][69] and is expressed by Equation ( 4): Where N is either the discrete frequency of magnitudes M in each magnitude bin, or the cumulative frequency of magnitudes ≥M; α, b are parameters determined by the earthquake catalog data, and the slope b of the G-R law expresses the relative number of the small magnitude earthquakes to the large magnitude ones and α is a measure of the seismicity level.In terms of geophysics, the b-value is considered as a possible stress meter, e.g., [58].In this context, low b indicates high material heterogeneity and concen-  3).(Right) Coulomb stress changes along the cross-sections AB, CD.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 6.0 main shock.The green lines show the surface projections of the two fault models.

Frequency-Magnitude Scaling Properties of the Foreshock and Aftershock Sequences in Terms of Non-Extensive Statistical Physics
The frequency-magnitude distribution (FMD), known as G-R law [56], is of vital importance for the characterization of a seismic sequence [56][57][58][59][60][61][62][63][64][65][66][67][68][69] and is expressed by Equation ( 4): where N is either the discrete frequency of magnitudes M in each magnitude bin, or the cumulative frequency of magnitudes ≥M; α, b are parameters determined by the earthquake catalog data, and the slope b of the G-R law expresses the relative number of the small magnitude earthquakes to the large magnitude ones and α is a measure of the seismicity level.In terms of geophysics, the b-value is considered as a possible stress meter, e.g., [58].In this context, low b indicates high material heterogeneity and concentrated stress while high b implies asymmetrically distributed stress.For details on the b-values for the Arcalochori seismic sequence see [21].
An alternative model that describes the frequency-size distribution of earthquakes from the first principles of non-extensive statistical physics (NESP) was introduced in [70].This model is based on Tsallis Entropy, a generalization of the Boltzmann-Gibbs one [71] that offers a coherent theoretical framework for analyzing complex dynamical systems with fractal features and long-range correlations [72].The approach introduced in [70] considers that the seismic energy E is related to the size of the fragments that fill the space between the activated fault blocks.Then, by considering that the earthquake magnitude is related to the radiated seismic energy as M = 2 3 logE, the cumulative distribution function N(>M) of earthquake magnitudes M can be deduced as [73]: where q m is the entropic index and α m a model parameter that expresses the proportionality between the seismic energy and the size of the fragments.In [74] updated the derived equation to include the minimum earthquake magnitude M 0 in a seismic catalog, which now reads as: The derived model, in the form of Equation ( 6), has extensively been used to describe the earthquake dynamics of local and regional seismicity, (e.g., [75][76][77][78][79][80][81][82][83][84][85][86]).In comparison to the Gutenberg-Richter (G-R) scaling law, the FA model provides a reasonable explanation of recorded earthquake magnitudes over a wider range of scales, while the b-value may be obtained as a special case for values over a certain threshold magnitude [75.86] The result of the application of the NESP model (Equation ( 6)) to the observed cumulative distributions N(>M) of earthquake magnitudes for the foreshock and aftershock sequences of the Mw = 6.0 mainshock, as well as for the NE and SW spatial clusters, are shown in Figure 15.The model provides a good fit to the observed distributions over the entire range of magnitudes, for the model parameters referred to in Table 6 and Figure 15.The greater q m value for the SW cluster indicates greater tectonic instability in this region where the mainshock and the major aftershocks occurred.Table 6.Parameter values for the foreshock and aftershock sequences in Arkalochori, as well as for the NE and SW aftershock clusters.N is the number of events (with M ≥ M c ), Mc the magnitude of completeness, a m , q m the parameters of the NESP model (Equation ( 6)) and τ 0 , q τ the parameters of the q-exponential function for the inter-event time distribution (Equation ( 10)).6) is shown with the solid line, for the parameter values shown in the down left corner and Table 6.The dotted lines represent the 95% confidence intervals.

Aftershock Production Rate and Modelling
It has long been recognized that the number of aftershocks following a major event resembles a power law decay with the time that expresses the relaxation process after the mainshock.This mathematical relationship takes the form of the so-called modified Omori's scaling law [87,88]: where n(t) is the production rate of aftershocks () = ()  ⁄ , N(t) is the number of aftershocks in time t after the mainshock, K and c are constants and p the power law exponent that usually takes values in the range 0.9 < p < 1.6 [89].Moreover, the cumulative frequency of aftershocks N(t) is estimated from n(t) as:  6) is shown with the solid line, for the parameter values shown in the down left corner and Table 6.The dotted lines represent the 95% confidence intervals.

Temporal Properties of the Aftershock Sequence 7.1. Aftershock Production Rate and Modelling
It has long been recognized that the number of aftershocks following a major event resembles a power law decay with the time that expresses the relaxation process after the mainshock.This mathematical relationship takes the form of the so-called modified Omori's scaling law [87,88]: where n(t) is the production rate of aftershocks n(t) = dN(t)/dt, N(t) is the number of aftershocks in time t after the mainshock, K and c are constants and p the power law exponent that usually takes values in the range 0.9 < p < 1.6 [89].Moreover, the cumulative frequency of aftershocks N(t) is estimated from n(t) as: Herein, we use the previous equation to model the evolution of the aftershocks in terms of the cumulative frequency of aftershock activity that followed the Mw6.0 mainshock.We focus on the two major aftershocks clusters, the NE and SW spatial clusters that were discussed previously, and apply the maximum likelihood method to estimate the model parameters of Equation ( 8) [90].In Figure 16 the cumulative number of aftershocks (for M ≥ M c ) with time is shown for the two spatial clusters along with the modified Omori's law (Equation ( 8)), which generally provides a fair fit for the parameter values given in Table 7.However, large aftershocks may trigger secondary aftershock sequences embedded in the aftershock sequence of the mainshock.In this case, several Omori regimes may be used to model the aftershocks production rate n(t) [89][90][91]: (9)  where H(•) denotes a unit step function and t 2 , t 3 indicates the occurrence times of secondary aftershock sequences.In Figure 16, breaks are observed in the cumulative number of aftershocks for both spatial clusters that are associated with strong aftershocks and the generation of secondary aftershock sequences.Hence, we investigate if the composite model of Equation ( 9) fits better the observed distribution.By setting t 2 = 22.8 days and t 3 = 67.4days that designate the occurrence times of strong aftershocks following the mainshock for the NE cluster (Table 7), we find that the composite model provides a better fit to the observed distribution (Figure 16), which is further confirmed by the smaller Akaike Information Criterion (AIC) value in comparison to the single Omori regime (Table 7).7.
A similar result is obtained for the SW cluster.In this case, we used t2 = 24.1 days and t3 = 32.8days which mark the occurrence times of strong aftershocks in the SW cluster (Table 7).The composite model provides a better fit to the observed distribution, in comparison to the single Omori regime, for the parameter values given in Table 7.

The Interevent Times Distributions for the Foreshock and Aftershock Sequences
Furthermore, we study the temporal scaling properties of the foreshock and after-  7.
Table 7.The considered mainshock, the duration (in days), the number of events (N), and the MLE of the modified Omori formula parameters for the NE and SW aftershock clusters, along with their associated uncertainties.AIC is the estimated Akaike Information Criterion for each model.A similar result is obtained for the SW cluster.In this case, we used t 2 = 24.1 days and t 3 = 32.8days which mark the occurrence times of strong aftershocks in the SW cluster (Table 7).The composite model provides a better fit to the observed distribution, in comparison to the single Omori regime, for the parameter values given in Table 7.

The Interevent Times Distributions for the Foreshock and Aftershock Sequences
Furthermore, we study the temporal scaling properties of the foreshock and aftershock sequences by investigating the interevent times (or waiting times) distributions between the successive events.In this analysis, earthquakes are considered as a point process in time, marked by the magnitude of the event, with interevent times τ between the successive events defined as τ i = t i+1 − t i , where t i is the time of occurrence of the ith event {i = 1, 2, . . ., N − 1} and N the total number of events.First, we construct the cumulative distribution of the interevent times (M ≥ M c ) for the foreshock and aftershock sequences, as well as for the NE and SW spatial clusters.We then model the observed distributions with the q-exponential function, derived in the framework of NESP [72,82,84,86,[92][93][94][95][96].It has been shown in various studies that the q-exponential function appropriately describes the distribution of interevent times in global, regional, and volcanic earthquake activity, as well as in aftershock sequences [73][74][75][76][77][78][79][80][81][82][83][84][85][86][91][92][93][94][95][96].
If P(>τ) = N(>τ)/N 0 is the cumulative distribution function (CDF) of the interevent times, with N(>τ) the number of interevent times with a value greater than τ and N 0 their total number, then the q-exponential cumulative distribution is given by [92]: where τ 0 is a constant in time units and exp q (x) is the q-exponential function defined as: when 1 + (1 − q)x ≥ 0 and exp q (x) = 0 in all the other cases.Its inverse is the q-logarithmic function: ln q (x) = 1 1−q x 1−q − 1 .In the limit of q → 1, the q-exponential and q-logarithmic functions lead to the ordinary exponential and logarithmic functions, respectively.
In Figure 17, we present the cumulative distributions P(>τ) for the foreshock and aftershock sequences and for the two spatial clusters in log-log plots.In all four cases, the q-exponential cumulative distribution (Equation ( 10)) provides a good fit to the observed distributions for the parameter values given in Figure 17 and Table 6.This is further confirmed by the expected linear dependence of corresponding q-logarithmic distributions ln q P(>τ) with τ [77], which are shown in the right panels of Figure 17.In all cases, the qlogarithmic function describes the observed distributions with high correlation coefficients, shown in the corresponding panels.The high values of q τ (Table 6) indicate long-range temporal correlations in the evolution of the earthquake activity and further confirm the high q τ values observed in aftershock sequences [84,[96][97][98].

Scaling of the Aftershocks Focal Zone with Time
The growth of the aftershock focal zone with time can provide valuable information regarding the triggering mechanisms of earthquake migration.This migration pattern

Scaling of the Aftershocks Focal Zone with Time
The growth of the aftershock focal zone with time can provide valuable information regarding the triggering mechanisms of earthquake migration.This migration pattern observed in many cases frequently scales as the logarithm of time [97][98][99][100][101][102].Various studies, based on numerical simulations [103,104], as well as on real cases [97,105], suggest that this logarithmic migration pattern signifies that aftershocks migration is driven by afterslip.In this case, aftershocks are generated as the outcome of afterslip propagation along the activated fault.To anticipate the growth of the aftershock focal zone with time as the outcome of afterslip, ref. [106] have recently introduced a numerical model.In this model, asperities on a fault are stressed initially by regional creep occurring at a steady deformation rate during the inter-seismic period.As the mainshock occurs, some of the asperities slip co-seismically, transferring large positive Coulomb stresses to the surrounding creeping regions.During the post-seismic phase, the stress-loaded regions can accommodate large amounts of afterslip and when a critical level of afterslip is reached, aftershocks are triggered.Static stress changes in the model are thought to trigger aftershocks only during the early post-seismic phase so that most aftershocks are triggered by afterslip.In this case, the seismicity rate R(t) can then be proportional to the afterslip rate V(t) [103,106,107]: where V + and V L are the sliding velocity just after the end of co-seismic rupture and the long-term loading velocity after the mainshock and t r the duration of the post-seismic phase.Considering the previous equation, the seismicity rate R(t) can then be given by: where R + and R L are the seismicity rates just after the end of co-seismic rupture and the long-term one after the mainshock, respectively.If .τ is the stressing rate and ∆CFS the co-seismic Coulomb stress changes induced by the mainshock, then the parameters t r and R + are given by t r = A / .τ and R + = R L exp(∆CFS/A ), where A = (a − b)σ, with a and b the rate and state frictional parameters and σ the effective normal stress.For t/t r 1, Equation ( 13) yields a decay rate for R(t) proportional to 1/t, which is consistent with a modified Omori decay rate with p = 1 [87].
With the previous assumptions, the distribution of afterslip velocities can be deduced.Initially, a fault with only depth varying normal stress, stressing rate, and rheological parameter A is considered.If the initial Coulomb stress field varies with the strike direction x, then aftershocks migrate along x, forming the initial distribution of afterslip velocities.Then, the propagation velocity V p of the aftershocks focal zone, in the early stage of the post-seismic phase, which typically lasts several weeks or months after the mainshock, is given by [103]: The expansion of the aftershocks zone L a between time t i and t (t > t i ) is now given by: The latter equation manifests the logarithmic expansion of the aftershocks zone with time and for a smooth co-seismic Coulomb stress field, implies its slow migration.
Since the estimated co-seismic Coulomb stress field ∆CFS (Equation ( 3)) can be significantly different from the "real" one, ref. [103] suggested a mean Coulomb stress gradient to be used.In this case, Equation (15) becomes: In the latter equation, l c is the radius of the co-seismic rupture, ∆σ the mean value of the mean co-seismic stress drops and ζ a constant.For an idealized Coulomb stress field, ζ takes the value of 2.77 [103].
In the light of the previously described model, we investigated the scaling properties of the aftershock focal zone with time for the Arkalochori aftershock sequence.In the analysis, we used the relocated catalog for the events with M ≥ M c = 2.5 to estimate the mean distance of aftershocks from the mainshock ∆L a (t) with time t, along the horizontal dimensions.The result is shown in Figure 18, as a function of the logarithm with time.The expansion of the aftershocks zone becomes apparent, as ∆L a (t) grows systematically with time after the surpass of one day from the mainshock.This growth can well be described by the afterslip front (Equation ( 16)) for over a period of one hundred days (R 2 = 0.97).We note that the logarithmic time dependence starts almost after the first day from the main event, possibly suggesting that after that time the system starts to be driven by an afterslip process.
In the latter equation, lc is the radius of the co-seismic rupture, Δσ the mean value of the mean co-seismic stress drops and ζ a constant.For an idealized Coulomb stress field, ζ takes the value of 2.77 [103].
In the light of the previously described model, we investigated the scaling properties of the aftershock focal zone with time for the Arkalochori aftershock sequence.In the analysis, we used the relocated catalog for the events with M ≥ Mc = 2.5 to estimate the mean distance of aftershocks from the mainshock 〈 ()〉 with time t, along the horizontal dimensions.The result is shown in Figure 18, as a function of the logarithm with time.The expansion of the aftershocks zone becomes apparent, as 〈 ()〉 grows systematically with time after the surpass of one day from the mainshock.This growth can well be described by the afterslip front (Equation ( 16)) for over a period of one hundred days (R 2 = 0.97).We note that the logarithmic time dependence starts almost after the first day from the main event, possibly suggesting that after that time the system starts to be driven by an afterslip process.
Furthermore, from Equation ( 16) the rheological parameter Α′ can be determined once the slope sa of the afterslip front is known.From Figure 18, we obtain sa = 0.320 ± 0.003.Then, from Equation ( 16  Furthermore, from Equation ( 16) the rheological parameter A can be determined once the slope s a of the afterslip front is known.From Figure 18, we obtain s a = 0.320 ± 0.003.Then, from Equation ( 16), s a = d L a (t) ∆σ , where l c is the radius of the co-seismic rupture and ∆σ the co-seismic stress drop.For a simple model of circular rupture, l c can approximately be determined as l c = 7 2 M o ∆σ 1/3 [51].The average stress drop for normal fault earthquakes in Greece is ∆σ = 5.5 ± 1.5 MPa [108,109], while the mainshock's seismic moment is M o = 1.1 × 10 18 Nm (Table 4) [21][22][23].Then, we estimate the value of l c ≈ 8.9 km for the co-seismic rupture.For ζ = 2.77, the rheological parameter A takes the value of A ≈ 0.71 MPa, which is within the range 0.1-1 MPa of A values that are usually found [97,98,103].This value is considerably higher than other A values that were estimated for recent normal fault mainshocks in Greece, as [97] estimated A ≈ 0.041 MPa for the 2020 Mw7.0 Samos earthquake, while [98] estimated the value of A ≈ 0.29 MPa for the 2021 Mw6.3 Northern Thessaly earthquake.

Concluding Remarks
The 2021 Arkalochori earthquake is a characteristic event in the time history of Central Crete.The main event of 27 September 2021 Mw6.0 was a reminder that strong earthquakes do occur onshore Crete.In the present work the patterns and the scaling properties of the 2021-2022 earthquake sequence that occurred at Central Crete, are presented.
A relocation procedure has managed to improve the relative locations of the foreshock epicenters, which are concentrated in the vicinity of the 27 September mainshock.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 [110], the 2017 Kos [111], and the 2020 Samos [112][113][114] earthquakes.
Strong NW-SE and E-W oriented negative velocity anomalies predominate at both the upper and the lower crust of Central Crete.These are observed down to 15 km depth at the tomograms presented.At the depth slice of 5 km, a NE-SW-trending zone of negative body-wave velocity perturbations appears near the epicentral region of the Mw = 6.0 Arkalochori earthquake.This anomaly follows the mean distribution of Alluvial deposits and post-alpine sediments which are bounded by positive (~13%) body-wave velocity perturbations, possibly connected to the older post-alpine sediments of Viannos formation and the Mesozoic carbonate rocks to the east and south of Arkalochori basin, respectively.In the area north of Arkalochori, an E-W-trending anticorrelated pattern of negative Pand positive S-wave velocity anomalies are observed at the depth range of 5-10 km.This anomaly coincides with the eastward bending of Kastelli's normal fault direction, from the area east of Arkalochori towards the town of Malia.In the depth slice of 15 and 20 km, an almost NNE-SSW discontinuity of positive to the west and negative to the east V P anomalies is identified along the west-dipping Agnos normal fault.Furthermore, crosssections created in both primary (P) and secondary (S) wave velocity anomalies, reveal this west-dipping structure that may be related to Agnos high-angle (~60 • ) normal fault.
The co-seismic Coulomb stress variations caused by the strong earthquake of Mw = 6.0 and its major foreshock and aftershock, as well as the spatial distribution of foreshocks and aftershocks are presented.For all the major seismic events analyzed, the same pattern on the spatial distribution of the ∆CFS is observed, which indicates a stress decrease towards NW and SE and a stress increase towards NE and SW of the ruptured faults.The spatiotemporal evolution of the sequence indicated triggering of seismicity throughout most of the aftershock zone soon after the mainshock, attributed to co-seismic stress transfer, followed by slower migration towards its outer edges, indicating possible afterslip.
A straightforward interpretation of horizontal motion patterns derived from the DIn-SAR analysis is mainly in agreement with the aftershock distribution and clustering.The displacement vector on ARKL site show strong subsiding component of ~16 cm and significant eastward and northward horizontal motion.Both the E-W and the vertical components deduced by the GNSS analysis agree with the interferometric results.The overall image of the ground deformation in the epicentral area, based on both interferometric and GNSS results, define a strong subsiding pattern with substantial horizontal motional component compatible with normal faulting activated structure.In the period followed the mainshock and up to 30 April 2022, the two GNSS sites in Arkalochori and Heraklion city show increased velocity values, compared to the previous period, indicating that the post-seismic relaxation is continuous up to this date.
The mainshock was followed by numerous aftershocks, with the largest being an Mw5.3 event on 28 September 2021.Both the foreshock and the aftershock sequences

Figure 1 .
Figure 1.Seismicity rate in terms of events per day (blue vertical bars) and cumulative nu events (solid black line) during June 2021-January 2022 in the area of Arkalohori.The occur events with ML ≥ 4 is marked by red stars (ML magnitude in the red axis to the right).

Figure 1 .
Figure 1.Seismicity rate in terms of events per day (blue vertical bars) and cumulative number of events (solid black line) during June 2021-January 2022 in the area of Arkalohori.The occurrence of events with M L ≥ 4 is marked by red stars (M L magnitude in the red axis to the right).

Figure 2 .
Figure 2. Spatial distribution of the 2021-2022 Arkalochori sequence, for 4750 events that occurred during the period between 13 January 2021 to 31 January 2022.The locations of the permanent (red triangles) and the temporary (blue triangles) stations are presented.The M ≥ 4.0 earthquakes are depicted by yellow stars.Faults are marked as red lines (see text for details).On the top left corner, the location of Greece is indicated in the red triangle while on the bottom left one the study area is included in the red rectangle.

Figure 2 .
Figure 2. Spatial distribution of the 2021-2022 Arkalochori sequence, for 4750 events that occurred during the period between 13 January 2021 to 31 January 2022.The locations of the permanent (red triangles) and the temporary (blue triangles) stations are presented.The M ≥ 4.0 earthquakes are depicted by yellow stars.Faults are marked as red lines (see text for details).On the top left corner, the location of Greece is indicated in the red triangle while on the bottom left one the study area is included in the red rectangle.

Figure 3 .
Figure 3. Temporal distribution of the 2021-2022 Arkalochori sequence for 4750 events.Periods A, B, C, D, and E are marked with orange, yellow, green, purple, and violet circles, respectively.The M ≥ 4.0 earthquakes are depicted by yellow stars.Faults are marked as red lines (see text for details).

Figure 3 .
Figure 3. Temporal distribution of the 2021-2022 Arkalochori sequence for 4750 events.Periods A, B, C, D, and E are marked with orange, yellow, green, purple, and violet circles, respectively.The M ≥ 4.0 earthquakes are depicted by yellow stars.Faults are marked as red lines (see text for details).

Figure 4 .
Figure 4. (a) Location of the 2021-2022 Arkalochori sequence for 4750 events (b) relocated events of the aftershock sequence using hypoDD.The locations of the permanent (red triangles) and the temporary (blue triangles) are presented (see text for details).The M ≥ 4.0 earthquakes are depicted by yellow stars.Faults are marked as red lines (see text and [30] for details).

Figure 4 .
Figure 4. (a) Location of the 2021-2022 Arkalochori sequence for 4750 events (b) relocated events of the aftershock sequence using hypoDD.The locations of the permanent (red triangles) and the temporary (blue triangles) are presented (see text for details).The M ≥ 4.0 earthquakes are depicted by yellow stars.Faults are marked as red lines (see text and [30] for details).

Figure 5 .
Figure 5. Presentation of the map of the 5 cross-sections, 5 km wide (A-A', B-B', C-C', D-D E-E') on the left side of the figure and the results of the SSW-NNE oriented cross-sections (

Figure 5 .
Figure 5. Presentation of the map of the 5 cross-sections, 5 km wide (A-A', B-B', C-C', D-D', and E-E') on the left side of the figure and the results of the SSW-NNE oriented cross-sections (upper panel) and the results of WNW-ESE oriented cross-sections (lower panel) on the right side.The projection of the 27 September mainshock is depicted by yellow star on sections (B-B') and (C-C').
Figure 5. Presentation of the map of the 5 cross-sections, 5 km wide (A-A', B-B', C-C', D-D', and E-E') on the left side of the figure and the results of the SSW-NNE oriented cross-sections (upper panel) and the results of WNW-ESE oriented cross-sections (lower panel) on the right side.The projection of the 27 September mainshock is depicted by yellow star on sections (B-B') and (C-C').

Figure 8 .
Figure 8. Cross-sections of VP (%) perturbations.The cross-section traces are the same with the first four ones of Figure 5.

Figure 8 . 34 Figure 9 .
Figure 8. Cross-sections of V P (%) perturbations.The cross-section traces are the same with the first four ones of Figure 5.Appl.Sci.2022, 12, x FOR PEER REVIEW 12 of 34

Figure 9 .
Figure 9. Cross-sections of vs. (%) perturbations.The cross-section traces are the same with the first four ones of Figure 5.

Figure 9 .
Figure 9. Cross-sections of vs. (%) perturbations.The cross-section traces are the same with the first four ones of Figure 5.

Figure 10 .
Figure 10.Upper maps: wrapped ascending (left) and descending (right) co-seismic interferograms over the Arkalochori area.The interferograms are draped over shaded relief.Lower maps: co-seismic displacement maps generated using the ascending and the descending image pairs and draped over shaded relief.

Figure 10 . 34 Figure 11 .
Figure 10.Upper maps: wrapped ascending (left) and descending (right) co-seismic interferograms over the Arkalochori area.The interferograms are draped over shaded relief.Lower maps: co-seismic displacement maps generated using the ascending and the descending image pairs and draped over shaded relief.Appl.Sci.2022, 12, x FOR PEER REVIEW 15 of 34

Figure 11 .
Figure 11.Displacement maps for (a) the vertical (up-down) and (b) the E-W displacement components for the Arkalochori earthquake overlain by the earthquakes with magnitude greater than 3 and the active faults of the broader area.Positive values on the E-W component indicate eastward motion, while the negative ones describe westward motion.

Figure 12 .
Figure 12. (Up) Coulomb stress changes distribution due to Mw = 4.9 event (yellow star) at centroid depth of 8.0 km.The red rectangle indicates the fault model for the kinematics of Mw = 4.9, while the blue one is the projection of the fault model of Mw = 6.0 main shock (listed in Table3).(Down) Coulomb stress changes along the vertical cross-section AB.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 4.9 and before the major earthquake Mw = 6.0.The green lines show the surface projections of the two fault models.

Figure 12 .
Figure 12. (Up) Coulomb stress changes distribution due to Mw = 4.9 event (yellow star) at centroid depth of 8.0 km.The red rectangle indicates the fault model for the kinematics of Mw = 4.9, while the blue one is the projection of the fault model of Mw = 6.0 main shock (listed in Table3).(Down) Coulomb stress changes along the vertical cross-section AB.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 4.9 and before the major earthquake Mw = 6.0.The green lines show the surface projections of the two fault models.

Figure 13 .
Figure 13.(Left) Coulomb stress changes distribution due to Mw = 6.0 event (yellow star) at centroid depth of 10.0 km.The red rectangle indicates the fault model for the kinematics of Mw = 6.0, while the blue one is the projection of the fault model of Mw = 4.9 events (listed in Table3).(Right) Coulomb stress changes along the vertical cross-sections A-B, C-D, E-F, and the parallel cross-section G-H (from up to down).The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 6.0 main shock.The green lines show the surface projections of the two fault models.

Figure 13 .
Figure 13.(Left) Coulomb stress changes distribution due to Mw = 6.0 event (yellow star) at centroid depth of 10.0 km.The red rectangle indicates the fault model for the kinematics of Mw = 6.0, while the blue one is the projection of the fault model of Mw = 4.9 events (listed in Table3).(Right) Coulomb stress changes along the vertical cross-sections A-B, C-D, E-F, and the parallel cross-section G-H (from up to down).The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 6.0 main shock.The green lines show the surface projections of the two fault models.

Figure 14 .
Figure 14.(Left) Coulomb stress changes distribution due to Mw = 5.4 major aftershock (yellow star) at centroid depth of 9.0 km.The red rectangle indicates the fault model for the kinematics of Mw = 5.4, while the blue ones are the projections of the fault models of Mw = 4.9 and Mw = 6.0 events (listed in Table3).(Right) Coulomb stress changes along the cross-sections AB, CD.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 6.0 main shock.The green lines show the surface projections of the two fault models.

Figure 14 .
Figure 14.(Left) Coulomb stress changes distribution due to Mw = 5.4 major aftershock (yellow star) at centroid depth of 9.0 km.The red rectangle indicates the fault model for the kinematics of Mw = 5.4, while the blue ones are the projections of the fault models of Mw = 4.9 and Mw = 6.0 events (listed in Table3).(Right) Coulomb stress changes along the cross-sections AB, CD.The green circles are the relocated hypocenters of the aftershocks which occurred after the Mw = 6.0 main shock.The green lines show the surface projections of the two fault models.

Figure 15 .
Figure 15.The frequency-magnitude distribution of earthquakes (squares) for the (a) foreshock sequence, (b) aftershock sequence, (c) NE aftershocks cluster, (d) SW aftershocks cluster.The corresponding fit according to Equation (6) is shown with the solid line, for the parameter values shown in the down left corner and Table6.The dotted lines represent the 95% confidence intervals.

Figure 15 .
Figure 15.The frequency-magnitude distribution of earthquakes (squares) for the (a) foreshock sequence, (b) aftershock sequence, (c) NE aftershocks cluster, (d) SW aftershocks cluster.The corresponding fit according to Equation (6) is shown with the solid line, for the parameter values shown in the down left corner and Table6.The dotted lines represent the 95% confidence intervals.

Figure 16 .
Figure 16.The cumulative number of events for M ≥ Mc (symbols) with time that followed the Mw6.0 mainshock in the (a) NE cluster and (b) the SW cluster.The solid lines represent the composite model of three modified Omori regimes, while the dashed line the model for a single modified Omori regime, for the parameter values shown in Table7.

Figure 16 .
Figure 16.The cumulative number of events for M ≥ M c (symbols) with time that followed the Mw6.0 mainshock in the (a) NE cluster and (b) the SW cluster.The solid lines represent the composite model of three modified Omori regimes, while the dashed line the model for a single modified Omori regime, for the parameter values shown in Table7.

34 Figure 17 .
Figure 17.The cumulative distribution function P(>τ) of the inter-event times τ (in minutes) (left panels) and the corresponding q-logarithmic function (right panels), represented by circles, for the (a,b) foreshock sequence, (c,d) aftershock sequence, (e,f) NE aftershocks cluster, (g,h) SW aftershocks cluster.Fitting with the q-exponential function (Equation (10)) is shown with the solid lines, for the parameter values and the corresponding correlation coefficients shown in the down left corners.

Figure 17 .
Figure 17.The cumulative distribution function P(>τ) of the inter-event times τ (in minutes) (left panels) and the corresponding q-logarithmic function (right panels), represented by circles, for the (a,b) foreshock sequence, (c,d) aftershock sequence, (e,f) NE aftershocks cluster, (g,h) SW aftershocks cluster.Fitting with the q-exponential function (Equation (10)) is shown with the solid lines, for the parameter values and the corresponding correlation coefficients shown in the down left corners.

Figure 18 .
Figure 18.The average expansion (in km) of the aftershock zone as function of the logarithm of time (symbols) for Central Crete 2022, Mw6.0 aftershock sequence.The solid line represents the logarithmic growth of the aftershocks zone.

Table 2 .
Average absolute values of P-and S-wave residuals and their cumulative reduction percentage during the inversion of experimental data.

Table 3 .
The main characteristics of the SAR images used in this study.

Table 4 .
GNSS stations' coordinates and the respective vector of the co-seismic displacement.