The Upper Crustal Deformation Field of Greece Inferred From GPS Data and Its Correlation With Earthquake Occurrence

We present a new geodetic strain rate and rotation rate model for Greece that has been derived using a highly dense GPS velocity field. The spatial distribution and the resolved rates of the various velocity gradient tensor quantities provided updated constraints on the present‐day upper crustal deformation in the region and revealed new details not reported previously. The spatial distribution of the second invariant demonstrated that the overall magnitude of strain rates is highest across two well‐defined provinces. The first follows the North Anatolian Fault and its two branches within the north Aegean, crosses central Greece and through the Gulf of Corinth it terminates in western Greece, while the second encompasses the extensional province of western Turkey and the eastern Aegean Sea islands. Our estimates revealed that shearing affects some of the fault‐bounded grabens of central Greece that lie to the SW of the North Aegean Basin implying considerable oblique extension. We identified a narrow region of counterclockwise rotation whose location and kinematics have been induced by the net effect across the intersection of the clockwise rotating domains of western and central Greece. The Aegean microplate and the Anatolian plate are separated by a wide transition zone which accommodates the curved stretching of the entire plate system. In both edges of the Hellenic forearc the dominant mode of crustal strain is E‐W extension. We found that earthquakes of M ≥ 5.6 are spatially well‐correlated with high‐strain areas, indicating that strain rate mapping could be used to inform future probabilistic seismic hazard analyses.


Introduction
Greece comprises a complex natural laboratory within the Eastern Mediterranean whereby various modes of crustal deformation, namely plate convergence, strike-slip and Subduction-Transform-Edge-Propagator (STEP) faulting, continental collision and widespread extension, operate simultaneously within a relatively small geographical area (e.g., Goldsworthy et al., 2002;Jolivet et al., 2013;Özbakir et al., 2020;Sboras et al., 2017;Shaw & Jackson, 2010) (Figure 1).This tectonic activity results in high deformation and seismicity rates throughout the country, while the frequent occurrence of large-magnitude earthquakes classifies Greece as the most earthquake-prone country in Europe (Grünthal & Wahlström, 2012).The two main contributors to the complexity of deformation in the Eastern Mediterranean are (a) the westward tectonic escape of the Anatolian plate into the Aegean along the North Anatolian Fault (NAF) and (b) the convergence of the African and Eurasian plates along the Hellenic subduction system (Barbot & Weiss, 2021;Jolivet et al., 2021 and references therein).These simultaneously operating first-order geodynamic processes are primarily accommodated along major active faults, such as the upper plate splay-thrust faults of the Hellenic subduction system (Mouslopoulou et al., 2015;Saltogianni et al., 2020;Shaw et al., 2008;Özbakır et al., 2013) and the large right-lateral strike-slip faults which are developed at the termination of the Hellenic subduction system to the west (Mouslopoulou et al., 2020;Perouse et al., 2017;Sachpazi et al., 2000) as well as across the north Aegean Sea (Ferentinos et al., 2018;Sboras et al., 2017 and references therein).The aforementioned pattern of compressional and transcurrent deformation is superimposed by widespread extension due to slab rollback of the subducting Nubian plate in the south and the propagation of the NAF into the Aegean in the north (Rodriguez et al., 2023;Schlidgen et al., 2014).As a result, multidirectional normal faulting prevails across continental Greece, Peloponnese, Crete and the entire northern and eastern Aegean (Armijo et al., 1992;Chatzipetros et al., 2013;Nicol et al., 2020;Papanikolaou et al., 2002).Complex faulting patterns within the Aegean upper plate indicate that strain is not only currently accommodated by the well-defined larger structures but also by secondary structures which are traversing the Aegean crust onshore and offshore Greece.Therefore, the integration of multidisciplinary data sets that record crustal deformation over a range of temporal and spatial scales alongside the traditional seismologydriven information is key to acquiring more reliable seismic hazard estimates.
Tectonic strain is commonly accumulated on faults during interseismic periods to be released during earthquakes when fault strength is exceeded (Scholz, 1998).This is why interseismic strain rates have been associated with seismic activity and their detailed mapping through geodetic observations has been used to derive seismicity rate forecasts and quantitative hazard estimates (Bird & Kreemer, 2015;Stevens & Avouac, 2021;Zheng et al., 2018).Rates of crustal deformation are commonly measured using Global Positioning System (GPS) velocity solutions.In the past decades, a number of studies have used GPS velocities to constrain the crustal deformation in Greece.In this context, Kahle et al. (2000) carried out a large-scale monitoring of the deformation field in Eastern Mediterranean and Anatolia using velocities of campaign and permanent GPS sites, among which the former dominated the Greek data set.Kreemer et al. (2004) focused on central Greece, northern Aegean and western Anatolia and presented the principal axes and the second invariant of the strain rate field using a data set which included the velocities of Kahle et al. (2000) for Greece, along with additional campaign-mode velocities from central Greece.Hollenstein et al. (2008) calculated the strain rate field in Greece using a velocity field with roughly the same density as that of Kreemer et al. (2004), comprising mainly of reoccupied campaign GPS sites.Floyd et al. (2010) provided the principal axes of the strain rate tensor along with rotation rates utilizing a denser GPS velocity field around the Corinth Gulf and the Peloponnese but, similarly to the previous studies, the entire northern Greece was under-sampled.Perouse et al. (2012) presented, in a large-scale study, a strain rate model for the central and eastern Mediterranean based on a GPS velocity field that combined all previously published velocities in Greece into one consistent solution.Chousianitis et al. (2015) using exclusively GPS velocities from continuously operating stations, the majority of which were not used in the previous studies, evaluated strain rate and rotation rate tensors in central and western Greece. D'Agostino et al. (2020) also used GPS velocities from permanent sites and analyzed strain rate patterns as well as the relation between long-term and short-term rotations across the SW Balkans.
Although these studies greatly refined our understanding on the key tectonic structures that control the deformation in Greece, they clearly suffer from low data coverage throughout northern Greece and the Aegean Islands.In other words, previous studies provided adequate resolution only in central Greece and part of the Ionian Islands.Herein, we overcome this limitation and provide an improved picture of the present-day crustal deformational pattern of the entire Greek region by resolving the rates and spatial distribution of various velocity gradient tensor quantities at higher resolution than ever before.To this end, we align the latest published GPS velocity fields into a common reference frame to produce one consistent solution which, so far, offers the best possible coverage of the entire Greek region.Using this velocity data set, we derive a new strain rate and rotation rate model to better constrain the locus and quantify the present-day upper crustal deformation in Greece.The improved spatial resolution of the input data in regions that, until now, have been insufficiently sampled, allowed the identification of previously unresolved deformational patterns and the refinement of existing findings.Subsequently, we examine the spatial correlation of the current strain rates with seismicity and explore their potential to serve as an indicator of earthquake occurrence for events above a certain magnitude.In summary, our results offer new insights over two directions.First, we provide a new geodetic strain rate model for Greece, allowing a more comprehensive understanding of the distribution of elastic deformation in the upper crust.
Second, we demonstrate that the calculated strain rates can serve as a proxy to indicate regions where future earthquakes are likely to occur, justifying the usability of strain data in the creation of geodetically constrained seismicity models which could inform future probabilistic seismic hazard analyses.

GPS Velocity Field
To obtain a well-constrained strain rate field, we needed to use a velocity field that provides a dense coverage of the broader Greek region.To this end, we combined the GPS-derived velocities of Chousianitis et al. (2015), Briole et al. (2021), Ergintav et al. (2023), and Kurt et al. (2023) (Figure S1 in Supporting Information S1).The published velocities of these solutions account for coseismic as well as postseismic deformation where necessary, and their Figure 1.Regional tectonic setting of Greece and surrounding areas.Main tectonic structures and active faults compiled from the GreDaSS database (Caputo et al., 2012;Sboras, 2011), the DISS Working Group (2018), Goldsworthy et al. (2002), andPerouse et al. (2017).The main features of the Hellenic subduction system (i.e., the thrust front of the Mediterranean ridge, the backstop front and the major upper plate splay faults) along with the Apulian collision front are depicted with gray lines.AFS: Aedipsos fault system; AFZ: Arkitsa Fault Zone; AG: Amvrakikos Gulf uncertainties were estimated taking into account the presence of time-correlated noise in the GPS time series.However, the velocities of these four solutions had to be aligned into a common reference frame due to the fact that they were obtained by means of different analyses and were expressed at different reference frame realizations.To do so, we used the VELROT 1.01 code (Herring et al., 2018) and merged these solutions into a common Eurasiafixed frame on the basis of the rotation pole for Eurasia as defined by the ITRF2014 (Altamimi et al., 2017).Specifically, we estimated a six parameter transformation (i.e., six components of the rate of change of translation and rotation) between each velocity field through the minimization of the horizontal velocity residuals at common stations.As reference velocity field we used that from Briole et al. (2021), because it includes the highest number of common stations, and aligned the remaining velocity fields on it.The reference frames of the four velocity solutions, the number of common stations used in the alignment and the obtained root mean square (RMS) of the fits for the horizontal components of the common stations during the transformations, are all tabulated in Table 1.The final combined GPS velocity field with respect to stable Eurasia is listed in Table S1 and illustrated in Figure 2.

Quantification of Geodetic Deformation Rates
Crustal deformation can be quantified by decomposing the velocity gradient tensor into its symmetric and antisymmetric parts, which correspond to the strain rate tensor and the rotation rate tensor respectively, and analyzing the spatial distribution of its derivatives.The dense horizontal velocity field discussed in Section 2 comprises the basic input toward a more detailed description of how deformation is distributed over Greece.To achieve this, we implemented the weighted least squares method of Shen et al. (2015) to interpolate the discrete horizontal velocities into a continuous strain and rotation rate field over a uniform grid of 0.1°spacing.The chosen grid cell size was selected to be comparable to the mean intersite distance of the GPS velocity data, which is 13 km.The weighting function of this approach considers not only the distance of the GPS velocities to the interpolation sites, but also their spatial distribution.The distance-dependent weighting is implemented either through a Gaussian or a quadratic decay function, while the spatial weighting, that improves the interpolation where GPS velocities are heterogeneously distributed, is achieved either through an azimuthal weighting or a Voronoi cell weighting function.The weight reduction with distance in both functions of the distance-dependent weighting is controlled by a spatial smoothing parameter which has units of distance.The values of this parameter reflect the smoothing range of the strain rate field and are optimally calculated in relation to the in situ data strength.This is achieved by specifying an a priori weighting threshold value (W t ) so that the sum of the weighting functions at each interpolation site remains constant.In this way, the algorithm determines the optimum spatial smoothing parameter for each grid point according to the spatial density of the geodetic velocities.
We tested different combinations of the distance-dependent and spatially dependent weighting functions, as well as we used a range of weighting threshold values in an effort to determine the optimal model according to the scale of the problem that the present study addresses; that is, to capture all regional deformation patterns without emphasizing on specific and highly localized strain anomalies.In this context, we produced one set of solutions combining the Gaussian function for distance weighting with azimuthal span and Voronoi cell for spatial weighting, and a second set combining the quadratic decay function for distance weighting with the azimuthal span and Voronoi cell weighting functions.With both distance-dependent weighting functions we tested different weighting threshold values (W t = 1, 3, 6, 12, 18, and 24) so as to explore the impact of spatial smoothing on the interpolation result.In each solution, the optimum value of the spatial smoothing parameter that makes the sum of Note.The minimum and median observational time span of the GPS stations provided in each solution, the processing software, the initial and final reference frames, the number of common stations used to derive the transformation parameters and the resulted RMS fit between these stations are tabulated.
the weighting functions at every interpolation site equal to the W t value, was determined after searching all distances up to 150 km with an increment of 1 km.As expected, lower values of the weighting threshold yielded solutions which emphasized localized strain rate variations and generated more small-wavelength strain signals and anomalies, while larger weighting threshold values produced smoother strain rate fields.This is demonstrated in Figure S2 in Supporting Information S1 where we present the dilatation rates superimposed by the principal strain rates along with the corresponding distribution of the spatial smoothing parameter for progressively increasing W t thresholds after employing the Gaussian and azimuthal weighting functions.This set of solutions suggests that the increase in W t results in smoother and less sharp strain rate patterns where the very localized strain rate variations gradually smear out.This is because the spatial smoothing parameter increases with increase of the weighting threshold (right column of Figure S2 in Supporting Information S1).Accordingly, more GPS data are included in the calculations, leading to smoother solutions especially in areas of dense coverage.However, the less smoothed solutions that were obtained using small W t values (i.e. 3 and 6), suffer by an increase in noise level given that the low values of the spatial smoothing parameter (<25 km across Greek mainland) can easily give rise to artifacts, while at the same time an overestimation of strain rates can arise at certain areas with uneven and/or closely spaced data.Such low weighting threshold values are more suitable for inversions aiming at resolving strain values associated with specific structures and determining very localized strain gradients, objectives which are not the focus of the present study.On the contrary, large W t values (i.e., 18 and 24), tend to produce oversmoothed solutions with less sharp spatial gradients that, as is particularly evident in the W t = 24 solution, yield unrealistic low strain rate values in the highly straining areas due to the contribution of distant data in the calculations.
To further examine the impact of the different weighting threshold values in the strain rate results and determine the optimal model compensating the inherent trade-off between solution stability and resolution that characterizes every strain rate estimation (Maurer & Materna, 2023), we follow Shen et al. (2015) and inspect the differential strain rate fields of the second invariant which are produced as W t decreases from 24 to 3 (Figure S3 in Supporting Information S1).It can be observed that the strain rate pattern in Figure S3a in Supporting Information S1, which corresponds to the decrease of W t from 24 to 18, is smooth and the majority of the differential strain rates occurs across well-know fast deforming regions (e.g., western Greece, Gulf of Corinth, NAF and north Aegean Sea).A similar pattern is obtained when W t decreases from 18 to 12 (Figure S3b in Supporting Information S1) where it is illustrated that the majority of the differential strain rates are resolved across the same areas as in Figure S3a in Supporting Information S1.However, as W t decreases from 12 to 6 (Figure S3c in Supporting Information S1), the differential strain rate pattern starts to deteriorate given that a large number of additional signals emerges throughout the continental part of the study area.These localized spikes of strain rates, which in many cases exceed 30 nstrain/yr, occur not only at the fast straining regions, but also at many slow straining regions where no significant sources of deformation are present, indicating a significant increase in the noise level of the model with W t = 6.Further decrease in W t from 6 to 3 produces an even noisier model which is dominated by localized concentrations of strain rates (Figure S3d in Supporting Information S1).In view of the above and considering the regional-scale of the present study, we select the solution that was obtained using a W t value of 12 as the optimal representation of the present-day crustal deformation pattern in Greece, since it manages to achieve considerable spatial strain rate variations and avoid widespread strong smoothing.The corresponding velocity gradient tensor quantities of this solution are shown in Figure 3.The presented interpolation result was achieved based on the distribution of the spatial smoothing parameter which is illustrated in Figure S4 in Supporting Information S1.For most of mainland Greece and its adjacent regions to the north as well as at western Turkey, this parameter remains well below 45 km, reflecting the consistency of the velocity field interpolation and the relatively low degree of smoothing throughout a wide region.
Following the approach that is described above for the combination of Gaussian and azimuthal weighting functions, we determined the optimal solutions for the other three combinations of data weighting functions which are incorporated in the Shen et al. (2015) code.In this context, we came up to model with W t = 12 when employed Gaussian and Voronoi cell weighting, model with W t = 3 for quadratic and azimuthal data weighting, and model with W t = 3 when we used quadratic and Voronoi cell weighting (Figure S5 in Supporting Information S1).As expected, W t got lower values in the two models with quadratic distance weighting because the quadratic function is much more conservative in weight reduction with distance compared to the Gaussian function (see Shen et al., 2015 for technical details).It is evident from Figure S5 in Supporting Information S1 that the aforementioned models gave very similar results compared to those of Figure 3 in terms of spatial patterns and amplitudes.This consistency underlines the robustness of the deformation field presented in Figure 3 given that even when using different weighting schemes, which differ in the implementation of weight reduction and weight scaling with distance and data distribution respectively, the details of active deformation in Greece remain essentially the same.Additionally, Figure 4 shows the formal uncertainties of the strain and rotation rate products of Figure 3, where it can be observed that throughout the study area the errors are fairly consistent across the study area with median values equal to 8 nstrain/yr, 11 nstrain/yr, 6 nstrain/yr and 0.3°/Myr for the second invariant of horizontal strain rates, dilatation rates, maximum shear strain rates and rotation rates, respectively.Furthermore, as Shen et al. (2015) pointed out, the distribution of the spatial smoothing parameter highlights the spatial resolution of the derived strain model and reflects the quality of the calculations.In this context, Figure S4 in Supporting Information S1 implies that the resolution is consistently high across the entire continental part of the study area where the distribution of the smoothing parameter is characterized by low values.As compared to the most recent strain rate estimations which include Greece and its surroundings, our new model has been derived using GPS velocities that have mean intersite distance equal to 13 km, while across the extent of our study area, this metric was found equal to 44 km in D' Agostino et al. (2020), 17 km in Sparacino et al. (2022), 28 km in Serpelloni et al. (2022) and 25 km in Piña-Valdés et al. (2022).Furthermore, the median value of the smoothing parameter of Figure S4 in Supporting Information S1 is 44 km while in the models presented by D' Agostino et al. (2020) and Piña-Valdés et al. (2022), which also used the Shen et al. (2015) approach, this value was 65 and 76 km, respectively.In what follows, we take advantage of the information incorporated into the velocity gradient tensor quantities presented in Figure 3, to offer insights into the current distribution and patterns of active deformation in the upper crust of the Aegean plate and discuss their tectonic implications.

Principal Strain Rates
The eigenvectors and eigenvalues of the strain rate tensor correspond to the principal axes of the strain rate field and define the coordinate system where shear strains are absent.Strain along these directions equals to the principal strains whose amount and direction provide a quantitative description of the greatest contractional and extensional strain rates.The map in Figure 3a highlights that extensional strain axes are dominant in large areas of Greece and western Turkey.In contrast, the magnitude of the principal compressional axes is considerable only along the large strike-slip structures of north Aegean and western Greece, where they coexist with the extensional axes and are oriented at ∼45°to the strike of these fault systems.The largest extension within the study area occurs along the western Gulf of Corinth where it reaches 270 nstrain/yr.A 50% rate decrease is observed at the eastern section of the Gulf of Corinth, where extension gradually drops to 130 nstrain/yr.Directly to the east of the Corinth Rift, extension rates decrease rapidly (by up to 70%-80% over a distance of ∼100 km) to mark the eastern termination of the rift onshore, within the Metropolitan Athens (Konstantinou et al., 2020).The Cyclades experience low to negligible rates on the order of ∼0-20 nstrain/yr, while the extensional principal strains become again significant along the eastern Aegean islands and the western coastline of Turkey.At this region, the extensional axes experience a counterclockwise rotation of ∼60°(from NE-SW to NW-SE) southward.More precisely, the principal strains at the island of Lesvos and further eastwards, toward the Edremit Gulf and northwestern Turkey, indicate NE-SW extension with significant NW-SE shortening, possibly due to their proximity to the NE-SW oriented strike-slip faulting related to the southern branch of the NAF.South of Chios/ Karaburun peninsula, shortening rate axes become non-significant implying the dominance of extensional tectonics.Near the island of Samos and the Menderes Graben (MeG in Figure 1), the principal extensional axes rotate counterclockwise to become N-S while further south, near the island of Kos and the Gulf of Gökova, they swing toward NW-SE.Throughout the eastern Aegean and western Turkey extension rates vary between 40 and 130 nstrain/yr, reaching maximum values at the Menderes Graben, where the data suggest significant localized extension.
North of the Corinth Rift and the Attica peninsula we document widespread active extension across a distance of ∼500 km, from central Greece to the eastern part of the Republic of North Macedonia and western Bulgaria.However, this extension appears to vary spatially within continental Greece (from south to north): extensional axes first decrease from 130 nstrain/yr within the eastern Corinth Rift to 60 nstrain/yr in Viotia, to increase again to 80 nstrain/yr along the north Evia Gulf rift and 100 nstrain/yr at Pagasitikos Gulf.All these extensional provinces are associated with historic large-magnitude (M > 6) earthquakes on E-W trending normal faults (Ambraseys & Jackson, 1990;Jackson et al., 1982;Papazachos et al., 1983).North of Pagasitikos Gulf, the magnitude of the extensional axes reduces to 40 nstrain/yr within the Mygdonia basin in the central Macedonia, where E-W trending active normal faults have produced a number of strong earthquakes (Caputo et al., 2012;Papanikolaou et al., 2022 and references therein), while further north across the eastern part of the Republic of North Macedonia and western Bulgaria, extension rates decrease by 50%, to about ∼20 nstrain/yr.Thus, the general direction of the extensional axes throughout the aforementioned province is N-S to NNE-SSW.Interestingly, in the vicinity of the broader epicentral area of the 2021 Tyrnavos seismic sequence the extensional axes that we calculated follow a NE-SW direction which deviates from the general N-S orientation of extension in Thessaly that has been estimated from geological data (Caputo & Pavlides, 1993).However, this NE-SW direction is compatible with the NW-SE strike of the coseismic ground ruptures that were revealed by postearthquake field observations and InSAR data (Chatzipetros et al., 2021;Kontoes et al., 2022;Koukouvelas et al., 2021;Mouslopoulou et al., 2022;Sboras et al., 2022) and is in agreement with the NE-SW extensional stress field obtained from the focal mechanisms of the 2021 Tyrnavos seismic sequence (Kassaras et al., 2022) and with geomorphic indicators which suggest repeated past ruptures on the 2021 causative faults (Mouslopoulou et al., 2022).In this context, it turns out that the orientation of the faults that ruptured during the 2021 sequence was not incompatible with the present-day crustal strain orientation as it was initially believed (Lazos et al., 2021;Sboras et al., 2022) and it appears that the occurrence of this earthquake sequence is well-justified by the active strain field.
Across the three southern Balkan countries that border Greece in the north, we document a clockwise rotation of the extension axis from NW-SE in Albania to NE-SW in Bulgaria.Contraction across eastern Albania, the Republic of North Macedonia and Bulgaria is not significant, and the magnitude of the principal compressional axes become considerable only within central and western Albania, where focal mechanisms indicate thrust and strike-slip faulting (Copley et al., 2009).In particular, along the coastal Albania principal strains show a combination of NW-SE extension and NE-SW compression, indicative of a transpressive regime.We note that the compressional axis becomes notably larger compared to the corresponding extensional axis in the epicentral area of the 2019 M w 6.3 Durrës earthquake.At this region we resolve a NE-SW oriented compression of 45 nstrain/yr, which is consistent with the NW-SE trending reverse fault that ruptured during the 2019 mainshock (e.g., Govorčin et al., 2020;Pezzo et al., 2022).The other region where the strain rate field is dominated by compression lies in northwestern Greece and, more specifically, in western Epirus and Paxoi Islands.There, we calculated an ENE-WSW directed compression of ∼90 nstrain/yr, consistent with shortening across the Apulian-Eurasian collision zone.In line with the geodetically derived compression in this area, Louvari et al. (2001) reports focal mechanism solutions with NE-SW-to ENE-WSW-trending σ 1 axes, as well as on 21 March 2020 a blind rupture on a NW-SE-striking east-dipping thrust fault was responsible for the M w 5.6 west Epirus earthquake.
To the south of the Corinth Rift, we resolve extension throughout Peloponnese which, if we overlook its eastern part that is characterized by very low strain rates, ranges from 40 to 70 nstrain/yr.The most striking feature is, however, the clockwise rotation of the extension direction.Indeed, the ∼N-S oriented principal extensional axes in the Gulf of Corinth change to NE-SW in central Peloponnese and to roughly E-W in southern Peloponnese.At the latter region which extends from the Gulf of Laconia to SE of Kythira island, this E-W extension is on the order of 40-50 nstrain/yr and is primarily accommodated by large ∼NNW-SSE oriented active normal faults which have been mapped, both onshore and offshore, using geological, geophysical, bathymetric and seismological data sets (Armijo et al., 1992;de Gelder et al., 2022;Lyberis et al., 1982;Veliz-Borel et al., 2022).
Although normal faults of various orientations (i.e., E-W, NW-SE, etc) traverse southern Peloponnese, Kythira and the seabed across them, the persistence of E-W extension implies that it is the N-S striking faults which play a dominant role in accommodating the present-day upper crustal deformation along the western forearc of the Hellenic subduction system.A similar switch by ∼90°in the extension direction (i.e., from N-S to ∼E-W) is also observed along the eastern section of the Hellenic forearc, near the island of Karpathos.E-W extension proximal to Karpathos is consistent with fault plane solutions, striations on active normal fault planes and the presence of steep NNE-SSW oriented bathymetric escarpments (some exceeding 2 km in size) around the island.
Finally, as regards Crete, the most striking feature of the strain rate field is the prevalence of active ∼N-S extension at the southern part of the island.Nicol et al. (2020) have highlighted a spatial separation of the orientation of normal faulting on Crete, with E-W trending faults extending mostly along the southern coastline, in agreement with our strain rate field.The peak values are resolved at the Messara Basin (MeB in Figure 1), a graben which is controlled by E-W and ENE-WSW trending faults that are still active today (Caputo et al., 2010;Nicol et al., 2020)

Dilatation Rates
The dilatation rate is given by the trace of the strain rate tensor and equals to Δ = εxx + εyy where εxx and εyy are the normal strain rates along the x and y directions, respectively.The dilatation rate is the first invariant of the strain rate tensor, so it does not depend on orientation of coordinate axes, and is associated with the rate of volume change.In a two-dimensional space the dilatation rate indicates contractional and extensional deformation and can therefore be associated with areas dominated by reverse and normal faulting, respectively.
Dilatation rates are pronounced along the Gulf of Corinth, where they peak at its western section (∼300 nstrain/yr) (Figure 3b).Extension along this active asymmetric rift is ∼10-15 mm/yr (Avallone et al., 2004;Briole et al., 2000;Chousianitis et al., 2013) and is primarily accommodated by a system of large active north-dipping E-W and WNW-trending normal faults.This strong positive dilatation along the Gulf of Corinth dominates the current pattern of dilatation rates and highlights its major role in the regional geodynamics of mainland Greece.Significant dilatation rates, on the order of 100 nstrain/yr, are also accommodated directly northeast of the Gulf of Corinth, across a zone that includes the regions of Fthiotida, north Evia Gulf and Pagasitikos Gulf.Interestingly, this extension zone coincides with the three major rift basins in central Greece, that is, the Sperchios graben, and the marine basins of north Evia Gulf and Pagasitikos Gulf (e.g., Goldsworthy et al., 2002;Roberts & Jackson, 1991;Walker et al., 2010).This feature provides strong evidence for active rifting in the aforementioned areas, which are located between the transtensional regime of the North Aegean Basin (NAB) and the strong extensional province of the Corinth Rift.
Two further regions that accommodate similar amounts of dilatation rates (∼100 nstrain/yr) with those resolved on the rift basins of central Greece are the Menderes Graben, a well-known E-W trending active rift basin in southwestern Turkey, and the offshore basin between Rhodes and Karpathos.Both areas host significant E-W and ∼N-S trending active normal faults, respectively (Aktug et al., 2009;Sakellariou & Tsampouraki-Kraounaki, 2019).The extension along the edge of the eastern Hellenic forearc (Karpathos, Rhodes) likely reflects arcparallel stretching in response to the African-Eurasian convergence.Furthermore, it is evident from Figure 3b that the broader region of Dodecanese islands, that is, from Samos in the north to Karpathos and Rhodes in the south, as well as from southwestern Turkey to Naxos and Amorgos in the Cyclades islands, comprises a broad area of positive dilatation with values ranging from ∼50 nstrain/yr to 100 nstrain/yr.This extension, that manifests itself through sharp E-W and NE-SW trending bathymetric escarpments within the eastern Aegean (Nomikou et al., 2018), accommodates partly the westward "escape" of the Anatolian plate (Bozkurt, 2000).
On the other hand, a ca.NNE-SSW narrow zone of zero dilatation across central Cyclades, which widens in the Cretan Sea, turns into slow dilatation onshore Crete.However, despite the low-level of upper plate deformation, the island's central part was stricken on 27 September 2021 by the M w 6.0 Arkalochori earthquake that ruptured a NNE-SSE-striking normal fault (Vassilakis et al., 2022).Prior to this, no significant shallow earthquake with M > 5.0 had occurred onshore Crete, at least during the instrumental period, although at least 25 roughly N-S striking active normal faults with postglacial activity have been mapped along the northern coastline of the island (Nicol et al., 2020).This slow deforming area may express a slow stress loading region, where stress is released with strong (i.e., characteristic) events with long recurrence intervals.A similar hypothesis can be made for the slow straining area of northeastern Greece, where a large angular, and a smaller, but yet quite prominent, generally E-W-striking normal fault zone (i.e., the "Thrace" and the "Drama" fault zones, respectively) were responsible for at least two strong (M ≥ 6.0), undeniable historic events that occurred in 1784 and 1829 (see Caputo et al., 2012 and references therein for further details).The very small dilatation rates evident in Figure 3b and the low rate of instrumental seismicity in this region may imply, as in the case of central Crete, long recurrence intervals of the large-magnitude earthquakes.
Dilatation rates in central Ionian islands (i.e., Cephalonia and Lefkada) are negative, indicating an excess of horizontal contraction, which becomes stronger onshore Lefkada in accordance to previous findings (Chousianitis et al., 2015;Ganas, Marinou, et al., 2013).Further north, significant compression, which amounts up to 80 nstrain/yr, is also found across a zone that extends from Corfu and Paxoi islands in the west to Epirus in the east.This ENE-WSW trending compressive zone can be associated with the active collision of the Apulian continental lithosphere against northwestern Greece, given that the western boundary of this zone matches the NNW-SSE oriented collision front and the compression is developed toward the ENE, that is, across the direction of motion of the Apulian continental block.Negative dilatation rates are maintained across Epirus where they progressively decrease and die out approximately 60 km inland.The resolved compressional regime onshore Epirus highlights the dominant role of the NNW-SSE trending thrust faults in this region of tectonic diversity, where the aforementioned thrusts are intersected by large E-W trending strike-slip faults and smaller normal and obliquenormal faults (Hatzfeld et al., 1995;Ntokos, 2018;Tselentis et al., 2006).Immediately to the south of the compressional regime of Epirus, dilatation rates become again positive, indicating the cease of active contraction resulting from the Apulian collision and the establishment of extensional tectonics associated with the E-Wtrending half-graben of Amvrakikos Gulf.

Maximum Shear Strain Rates
The maximum shear strain rate is defined as γmax = √ where εxx and εyy are the normal strain rates along the x and y directions, and εxy and is the shear strain acting on the xy plane.The maximum shear strain rate is associated with localized shear deformation and implies increased rupture likelihood for faults with the same orientation as the direction of the maximum shear strain rates.It also highlights regions with dense fault distribution.
The distribution of the maximum shear strain rates in Figure 3c points to faults/fault systems associated either with rifting or with strike-slip tectonics.The former is primarily related to the Corinth Rift which accommodates high-magnitude shear strain rates that range from 60 nstrain/yr in the east to 280 nstrain/yr in the west.The westward increase of maximum shear strain rates within the Corinth Rift is accompanied by increased microseismicity and more frequent earthquake swarms (e.g., Chouliaras et al., 2015;Ganas, Chousianitis, et al., 2013;Kaviris et al., 2021;Pacchiani & Lyon-Caen, 2010).Although both parts of the Corinth Rift experience some creep (Chousianitis et al., 2015), the different rheological conditions and the presence of a low-angle active layer below the rift (Bernard et al., 2006;Bourouis & Cornet, 2009;Kaviris et al., 2018;Rigo et al., 1996) seem to contribute to the observed high seismicity rates of the western Gulf of Corinth.Smaller amounts of shear strain rates, on the order of 60 nstrain/yr, are found along the north Evia Gulf rift and the Pagasitikos Gulf in central Greece, as well as at the Menderes graben in SW Turkey.These locations comprise active rift systems where extension is accommodated by a combination of localized dilatation (see Figure 3a) and shear deformation.Of particular interest is the former region, where shear strain rates are mainly associated with the normal fault zone of Arkitsa (AFZ in Figure 1) and the Aedipsos fault system (AFS in Figure 1) which currently appears to involve significant strike-slip component (Caroir et al., 2023).In conjunction with other smaller structures onshore NW Evia island as well as offshore within the north Evia Gulf which are also characterized by oblique-to strike-slip kinematics (Ganas et al., 2016), we conclude that the maximum shear strain rates resolved in this region imply a significant amount of oblique extension.
The rest of the locations with significant maximum shear strain rates correlate with major strike-slip faults in Greece and Turkey.The most striking feature is the wide zone of high shear strain rates that encompasses the NAF in Turkey and the North Aegean Trough (NAT) and the NAB in Greece.The NAF accommodates >100 nstrain/yr of shear across the Sea of Marmara and within the Gulf of Saros, in accordance with the results of Weiss et al. (2020), immediately before it enters the north Aegean Sea.The NAT, which comprises a large transtensional shear zone west of the propagating tip of the NAF, accommodates a smaller, but yet, significant amount of shear strain rate (80-100 nstrain/yr).The 24 May 2014 M w 6.9 North Aegean earthquake revealed moment tensor solutions highly compatible with this shear regime (Evangelidis, 2015;Konca et al., 2018).The NAB, which lies west of the NAT, accommodates 60-80 nstrain/yr.It is notable, however, that west of the central Sporades islands, approximately at 23.8°E, shearing within the NAB decreases implying the prevalence of extension over transtension, as also suggested by the distribution of dilatation rates (Figure 3b).The southwestern and southern limits of the wide zone of high maximum shear strain rates within the north Aegean Sea coincide with the sinistral NW-SE trending Skyros fault and the dextral NE-SW trending southern branch of the NAF that is entering the Aegean, respectively.The Skyros Basin forms at the junction of these two conjugate strike-slip faults, in a similar manner to that of the NAB, albeit it is characterized by shallower depths and lower seismic activity (Mascle & Martin, 1990;Papanikolaou et al., 2002Papanikolaou et al., , 2019)).Maximum shear strain rates of 60-80 nstrain/yr and low dilatation rates within the Skyros Basin suggest that, although both basins are characterized by oblique opening and host normal and strike-slip faults, the Skyros Basin accommodates primarily shearing with only small amounts of extension.This differentiates the Skyros Basin from the NAB whose western part accommodates primarily extension and only its eastern part is dominated by shearing.Lastly, at NE Aegean Sea, maximum shear strain rates of roughly 60 nstrain/yr are found between Lesvos island and Karaburun peninsula.This illustrates the accumulation of notable amounts of shearing at this area of transtensional tectonics where normal and strike-slip faulting coexist.The importance of strike-slip tectonics was highlighted during the 2017 M w 6.3 Lesvos earthquake, where strike-slip faults accommodated, in conjunction with normal faulting, the relative motion of Lesvos island with respect to the Karaburun Peninsula (Chousianitis & Konca, 2018;Nomikou et al., 2021).Finally, the central and southern Aegean (from south Evia to Crete and Karpathos), appear to accommodate insignificant shear strain rates.
As we move to western Greece, the high shear strain rates of the western tip of the Gulf of Corinth continue unceasingly to the northwest, forming a very clear zone of high shear strain rates along the NNW-trending Katouna Fault System (KFS in Figure 1; Pérouse et al., 2017;Underhill, 1989).This left-lateral structure intersects at high angles (∼90°) the western part of the zone of influence of the Corinth Rift (evident in Figures 3b  and 3c), while extension rates decrease considerably west of this junction (Figure 3b).This is possibly because slip is transferred from the rift onto the strike-slip system through their mutual intersection, facilitating the former to terminate (Mouslopoulou et al., 2007).The 100-240 nstrain/yr which are accommodated across the KFS, emphasize the dominant role of the coexistent active pure left-lateral shear and left-lateral transtension in the kinematics of this region (Haslinger et al., 1999;Perouse et al., 2017;Vassilakis et al., 2011) and confirm, in terms of strain rates, the behavior of this fault system as a boundary that separates pure extensional tectonics in the east from strike-slip dominated tectonics in the west.A second less profound shear zone is recorded in NW Peloponnese.With maximum shear strain rates exceeding 60 nstrain/yr, the NW part of Peloponnese shears considerably compared to the rest, where rates are either low (south and southeast) or negligible (east).The shearing has a general NE-SW orientation and is thought to align with the Kato Achaia-Amaliada Fault Zone (KAFZ in Figure 1) which was responsible for the 2008 M w 6.4 strike-slip earthquake (Serpetsidaki et al., 2014 and references therein).Furthermore, Figure 3c shows that maximum shear strain rates of >60 nstrain/yr extend from NW Peloponnese to throughout the central Ionian islands, a region where focal mechanisms suggest distributed strike-slip faulting (Shaw & Jackson, 2010).Thus, seismological and geodetic data demonstrate that the area bounded by the Cephalonia Transform Fault Zone (CTFZ in Figure 1), the KAFZ and Zakynthos island, is a region of distributed shear.Finally, as suggested by the negative dilatation rates observed in the central Ionian islands (i.e., Cephalonia and Lefkada) and offshore NW Peloponnese (Figure 3b), it is clear that apart from pure strike-slip deformation, oblique thrusting should also be expected in this area (Chousianitis & Konca, 2019;Mouslopoulou et al., 2020;Sokos et al., 2020).

Second Invariant of Strain Rates
The second invariant of horizontal strain rates equals to İ2 = and as in the previous relationships, εxx , εyy , and εxy represent the strain rate components.The second invariant is a scalar quantity that can be considered a measure of the total strain rate field and correlates well with seismicity and increased seismic hazard (Bird & Kreemer, 2015;Elliott et al., 2016;Kreemer & Young, 2022;Stevens & Avouac, 2021).It can be observed that the resolved pattern of the second invariant of horizontal strain rates (Figure 3a) shares common features with both the dilatation rates and the maximum shear strain rates.In this context, Figure 3a displays the main characteristics of Figure 3c combined with additional features taken from Figure 3b.This underlines that the deformation field of Greece is not entirely dominated by shear strain as it has been observed elsewhere, such as continental China and California (e.g., Wang & Shen, 2020;Zeng, 2022), but dilatational strain also plays an important role.
The most striking feature of the second invariant of strain rates is a high-magnitude strain rate "ridge" (>50 nstrain/yr), which extends from the NAF in the east to the Ionian islands in the west, separating two lowmagnitude strain rate "depressions" (<20 nstrain/yr) that occupy N/NW Greece and central/southern Aegean Sea.In this context, Figure 3a reveals that the bulk of strain rate in Greece occurs as pervasive deformation across a zone that is dominated, from east to west, by shearing (across the NAF, the NAT and the NAB), extension (throughout the fault-bounded grabens of central Greece) and again shearing (across the Katouna fault system and throughout the region of distributed strike-slip faulting from KAFZ to CTFZ).This pattern remarkably reproduces, in terms of strain rates, the westward propagation of the NAF across the north Aegean which separates the central/southern Aegean province from Eurasia and allows its trenchward motion in a quasi-rigid manner.An interesting pattern that arises from a closer examination of Figure 3a relates to a narrow NW-SE trending zone in central Greece, immediately north of Pagasitikos Gulf, that suggests NE-SW extension (of 50-80 nstrain/yr) along the newly identified (Mouslopoulou et al., 2022) ∼100 km long relay fault zone that accommodated the three normal fault earthquakes during the 2021 Tyrnavos seismic sequence.
Furthermore, the NNW-SSE oriented wide zone of high-magnitude strain rates that is developed between the central/southern Aegean and the Anatolia effectively highlights the widespread active continental extension in western Turkey, which is related to onshore "basin & range" structures.We observe also that these rates are not restricted only along the extensional province of western Turkey, but extend offshore, across the eastern and northeastern Aegean where transtensional deformation is partitioned between ∼E-W normal and ∼NE-SW conjugate strike-slip faults.This indicates that the active structures in the eastern and northeastern Aegean accommodate significant strain and are associated with elevated earthquake hazard.This increased strain accumulation was notoriously expressed the previous years by the occurrence of the 2017 M w 6.3 Lesvos and M w 6.6 Kos earthquakes (e.g., Kiratzi, 2018;Konca et al., 2019;Sboras et al., 2020) and the 2020 M w 7.0 Samos earthquake (e.g., Chousianitis & Konca, 2021;Hu et al., 2022;Sboras et al., 2021).
Finally, in Figure 3a it is evident that throughout Albania, high strain rates are seen only in the epicentral area of the 2019 M w 6.3 Durrës earthquake, where negative dilatation rates indicate strong compressional deformation (Figure 3b) in accordance to the thrust-faulting character of this earthquake (e.g., Govorčin et al., 2020;Pezzo et al., 2022).

Rotation Rates
The rotation rate (vorticity) field can be obtained from the antisymmetric part of the velocity gradient tensor.It provides information about the spatial distribution, direction and magnitude of the present-day vertical axis rotation rates which comprise an important component of crustal deformation and active tectonics.Their comparison with paleomagnetic rotations can reveal whether specific rotational patterns persist today and provide insights into the evolution of deformation through the geological timescale.
The results in Figure 3d show that a number of clockwise and counterclockwise vertical axis rotations operate today, implying that the deformation field is largely block-like.A narrow band of fast clockwise rotation with rates between 3°/Myr and 6°/Myr is resolved along the NAF and its extension into the NE Aegean Sea.This pattern is consistent with right-lateral slip along this major strike-slip fault system.Clockwise rotations are preserved throughout northern Aegean, given that this area is primarily governed by right-lateral shear zones with strike-to oblique-slip kinematics, but with lower rates (∼2°/Myr).Moreover, the band of clockwise rotation that encompasses the NAF and the NE Aegean Sea borders the opposite-rotating western Turkey domain, where counterclockwise rates increase toward the south reaching ∼4°/Myr at SW Turkey.The boundary between these two domains is defined by the southern branch of the NAF which appears to comprise the northern limit of the broad counterclockwise rotation that characterizes the entire western Turkey and the eastern Aegean islands.
In central and western Greece, the present-day rotational pattern appears to be more complex than previously thought.Although the first-order patterns are similar to those discussed in earlier studies (Chousianitis et al., 2015;Hollenstein et al., 2008), the dense GPS site spacing allows here the identification of a region of coherent counterclockwise vertical axis rotation NW of the Corinth Rift, which is in sharp contrast to the otherwise clockwise dominated rotational field of central and western Greece.This clockwise rotation primarily occurs over two domains: (a) Across central Greece, North Evia and the broader area of Pagasitikos Gulf and (b) across the central and southern Ionian Sea and the NW tip of Peloponnese.The former domain comprises an extended coherent region of clockwise rotation that occupies 200 km in a NNE-SSW direction, with average rates of rotation ∼4°/Myr.This feature occurs in the rotation rate field because this domain accommodates the NE-SW shear originating from the right-lateral strike-slip faulting systems of the NAT and the southern branch of the NAF (Goldsworthy et al., 2002;Shaw & Jackson, 2010).Paleomagnetic data suggest that central Greece experienced ∼25°of post-Middle Pliocene clockwise rotation with estimated mean rates on the order of 8-10°/Myr (Bradley et al., 2013).By comparison, the vertical axis rotation rates that we obtained are at least half the paleomagnetic rate, indicating that although the rotation of this domain is still active today, it has been significantly reduced.
Immediately to the south, however, paleomagnetic and geodetic data are largely inconsistent with each other.Specifically, Mattei et al. (2004) provide paleomagnetic evidence for vertical axis rotations on the order of 19°for western Attica (Megara basin) and 15°for NE Peloponnese during Plio-Pleistocene (i.e., ∼5 Ma) which do not seem to be continuing today given that our GPS-derived rotation rates for these regions are quite low (<1°/Myr).
The other region which rotates clockwise about the vertical axis is associated with the large dextral strike-slip faults of western Greece, namely the CTFZ and the KAFZ, which fairly delineate the NW and SE boundaries of the rotating domain, respectively.This region exhibits clockwise rates of 5-8°/Myr, values significantly larger compared to those resolved in central Greece.Another feature that differentiates this rotating domain from that of central Greece, underlining the different nature of their driving mechanisms, is that the rotation in western Greece is accompanied by internal contraction and not extension, as in the case of central Greece.Paleomagnetic data suggest that the Ionian Sea together with the western part of the Greek mainland experienced a net clockwise rotation of 50°which was initiated in middle Miocene (i.e., 15-13 Ma) and occurred in two episodes (Duermeijer et al., 2000;Kissel et al., 2003;Kissel & Laj, 1988;van Hinsbergen et al., 2005).The most recent episode took place during Plio-Pleistocene and while it was initially believed that it was associated with a rapid clockwise rotation of 25°, a reevaluation led van Hinsbergen et al. (2005) to conclude that this large rate reflected primarily local phenomena and that the regional rotation during this phase was actually on the order of 10°.Our GPSderived rates for the western Greece rotational domain suggest that although decreased, they are similar to the paleomagnetic rates indicating that this clockwise rotation continues today and contemporary GPS data capture long-term motions.
Between the two aforementioned clockwise rotating domains, there is a narrow region that undergoes counterclockwise rotation with rates up to 4°/Myr which was not identified in previous strain rate models.The location and orientation of this zone is consistent with left-lateral slip on the NNW-trending KFS.The rotation rate field of Figure 3d shows that the formation and the kinematics of this structure are associated with the clockwise rotation of the two aforementioned domains in western Greece and central Greece, whose net effect across their intersection is to cause the opposite sense of rotation and consequently induce faulting with left-lateral sense of slip.In fact, the opposite rotating region resembles a simple three-gear train: one of the bilateral clockwise rotating regions (central/western Greece) could be considered as the driver gear and the other as the driven gear with counterclockwise motion.In view of that, the KFS which accommodates counterclockwise rotation and high shear rates (see Section 3.3), appears to be a notable tectonic boundary in continental western Greece.This was independently shown from block modeling (Vassilakis et al., 2011) and from GPS velocity analysis (Perouse et al., 2017), while with the present study we demonstrate the same outcome on the basis of strain and rotation rates.Whether or not an eastern boundary of this narrow counterclockwise rotating region with the central Greece domain exists is unclear, given that no active mapped faults are known from the literature.However, we believe that the possibility for the existence of blind and/or immature lateral slipping fault(s) cannot be ruled out.
In northern Greece and southern Bulgaria, a region of slow counterclockwise rotation is noticeable.The peak of this rotation, which amounts to about 1°/Myr, coincides with the NW-SE trending Strymon fault system (SFS in Figure 1) which has given evidence for left-lateral strike-slip kinematics (Lazos et al., 2022;Mouslopoulou et al., 2014).The counterclockwise rotation resolved here confirms its sinistral movement and taking into account the coexistence of positive dilatation rates (Figure 3b) which imply coupling with extension, it can be argued that oblique left-lateral kinematics characterize the deformation pattern across this region.

Implications for Large-Scale Tectonics in the Broader Greek Region
The different aspects and the principal short-and long-wavelength features of crustal deformation throughout Greece were presented and thoroughly discussed in Section 3.However, the mapping of the various velocity gradient tensor quantities provides some worth mentioning implications associated with large-scale tectonics acting across Greece and western Turkey.In this context, our findings suggest that shear deformation is mostly distributed over two regions which correlate with strike-slip structures, specifically (a) along a broad zone that includes the NAF in Turkey and its two branches within the north Aegean Sea and (b) throughout an area in western Greece that encompasses the central Ionian islands and is bounded by the CTFZ, the KAFZ and the KFS.Large shear strain rate values are not maintained across the central Greece mainland, as has already been pointed out by Perouse et al. ( 2012) and Chousianitis et al. (2015), and the two aforementioned regions are not merged via a high magnitude shear strain rate zone occupying the so-called Grecian Shear Zone (GSZ) or Central Hellenic Shear Zone (CHSZ) as was named by Şengör (1979; see their Figure 4) and Papanikolaou and Royden (2007; see their Figure 2), respectively (although in terms of the second invariant of strain rates they do connect, as we'll discuss below).Nevertheless, our dense GPS data set allowed the identification of sizable amounts of accumulated shear strain rates to the SW of the NAB, namely across the north Evia Gulf rift and the Pagasitikos Gulf in central Greece.This feature, which was not detected in earlier geodetic studies, implies that across this part of the extensional province of central Greece, where some important rift basins have been developed in an oblique direction to the strike of the large transcurrent structures of north Aegean, earthquakes with significant strike-slip component can also be expected.
Contrary to the aforementioned distribution of the shearing component, the overall magnitude of strain rates, as it is quantified by the second invariant, remains increased across the aforementioned zone in central Greece that passes through the GSZ (or in its alternate name, the CHSZ), in a way that high strain values continue unceasingly from north Aegean to western Gulf of Corinth and central Ionian islands (Figure 3a).Within this transition zone that lies between the subparallel strike-slip structures of north Aegean (i.e., NAF branches) and western Greece (i.e., CTFZ and KAFZ), our results suggest that the dominant type of deformation which is responsible for the increased magnitude of the strain rates is distributed ∼N-S extension accompanied by notable (∼4°/Myr) clockwise rotation (see Figures 3a and 3d).We note that tectonic and geodetic studies contradict each other on the subject of the dominant type of deformation that is accommodated across this transition zone in central Greece (Briole et al., 2021;Papanikolaou & Royden, 2007;Perouse et al., 2012;Royden & Papanikolaou, 2011).The former advocate the occurrence of significant dextral shear that coexists with extension, while the latter support either solely N-S extension or coupled with negligible dextral deformation, so that the connection of the NAF system with western Greece through a wide shear zone is not justified.However, our findings suggest that the actual situation appears to fall somewhere in the middle of these two end-member theories.The rotating distributed extension that we resolved in this transition zone in conjunction with the locally accumulated shear strain rates at north Evia Gulf rift and Pagasitikos Gulf imply oblique-slip/transtensional tectonics, but not with the typical sense, because the much stronger extensional component masks out the weaker signal of the strike-slip motions.Additionally, contrary to the majority of transtensional regimes, the dominant structures across this transition zone are the extensional (i.e., active normal faults that bound large graben systems) and not the strikeslip structures, given that no major transcurrent faults have been clearly mapped in the field and the recorded strike-slip seismicity in the area is associated only with small faults with limited potential.In conclusion, our strain rate mapping implies that although this region in central Greece accumulates mostly ∼N-S extensional deformation, it also hosts sizable amounts of transtensive deformation and the occurrence of strike-slip earthquakes to the SW of the NAB cannot be ruled out, at least as far as the north Evia Gulf rift.
Apart from the well-defined high-strain zone that follows the NAF and its two branches within the north Aegean, crosses central Greece and through the Gulf of Corinth it ultimately terminates in western Greece, the mapping of the second invariant of strain rates highlighted a second high-strain zone that encompasses the extensional province of western Turkey and the eastern Aegean Sea islands near the Turkish coastline.This NNW-SSE oriented zone differs from the aforementioned in that it is much wider and the deformation here, although less intense, is diffusely distributed without prominent localization along specific major structures and with rates that do not exhibit large variability throughout the extent of the zone.We interpret this pattern as the manifestation of a "soft" boundary between the Aegean microplate which moves rapidly to the SW with rates that reach ∼34 mm/yr (McClusky et al., 2000) and the Anatolian plate that lies to the east of this high-strain zone and moves westwards with rates on the order of 21 mm/yr (Reilinger et al., 2006).We use the characterization "soft" to describe a diffused tectonic boundary, not marked by a discrete fault zone, but consisting of multiple, smaller and scattered faults which individually accommodate most of the region's deformation.These faults are mainly pure ∼E-W normal structures onshore western Turkey, and as we approach the Turkish coasts and move westwards offshore, the typical pure ∼E-W normal faults coexist with conjugate ∼NE-SW strike-slip structures (Aktar et al., 2007;Aktuğ et al., 2009;Uzel et al., 2013;Yilmaz & Gürer, 2023).This diffuse zone of deformation and the faults that receive it could possibly reflect the fact that the Aegean microplate and the Anatolian plate are not separated in a straightforward way, that is, by means of a well-defined tectonic boundary, but with this wide transition zone which accommodates the curved stretching of the entire plate system whose direction and rate of motion changes as we pass from the westward moving Anatolian plate to the faster and SW moving Aegean microplate.
Our findings revealed that to the south of the two areas that currently accommodate the most intense extension, that is, the extensional provinces of central Greece and western Anatolia, we observe an identical pattern in which the principal extension rate axes rotate from ∼N-S to E-W.In the former region, this rotation takes place across the entire Peloponnese in a clockwise sense, while in the latter region the rotation starts somewhere between Samos and Kos islands and is implemented in a counterclockwise sense.Accordingly, it is illustrated that in both edges of the forearc of the Hellenic subduction system the dominant mode of crustal strain is E-W extension and hence the most important faults which accommodate the present-day active deformation are N-S trending normal faults.Across the central part of the forearc that is occupied by the island of Crete, the situation is less clear given that our strain rate model, in agreement with recent earthquake focal mechanisms (see Figure 1; Shaw & Jackson, 2010) and tectonic/geomorphological field observations (Caputo et al., 2010;Gallen et al., 2014;Nicol et al., 2020), highlighted localization of strain at certain regions and showed that both E-W and N-S striking extensional faults are presently active.Thus, although the entire forearc accumulates primarily extensional deformation, it seems that localized phenomena affect the strain rate field in its central part, while deformation across the forearc edges is driven by regional forces.Such forces include slab rollback of the subducting Nubian plate (van Hinsbergen & Schmid, 2012), incipient continental collision (Armijo et al., 1992;Lyon-Caen et al., 1988), NAF propagation (Armijo et al., 1999) and slab tearing at both edges of the Hellenic subduction system (Faccenna et al., 2014;Jolivet et al., 2013).

Correlation of Geodetic Strain Rates With Seismicity Distribution
During the last decade or so, geodetic strain rate mapping has become an integrated component in seismic hazard analyses and associated models especially in regions with dense geodetic coverage (Bird & Kreemer, 2015;Johnson et al., 2022;Shen et al., 2007;Stevens & Avouac, 2021;Ward, 2007).This is because strain rate estimates exhibit positive correlation with seismicity rates (Kreemer & Young, 2022;Shen et al., 2007;Zeng et al., 2018) and are often used to calibrate the location, magnitude and odds of future large-magnitude earthquakes (Riguzzi et al., 2012;Rollins & Avouac, 2019;Xiong et al., 2021;Zheng et al., 2018).Indeed, detailed strain rate mapping provides earthquake potential models which account for unmapped and/or too slow faults overcoming, thus, the limitations of traditional seismological approaches.In addition, agreement between geodetically derived slip rate estimates (i.e., short-term estimates using data sets that span a few years only) and geologically derived displacement rates (i.e., long-term estimates using data sets that span thousands to millions of years) has been observed for major faults and/or fault systems in various regions worldwide, such as in the Eastern Mediterranean (Reilinger et al., 2010), the Arabia-Eurasia collision zone (Allen et al., 2004) and the India-Eurasia collision zone (Zheng et al., 2017).This observed correlation implies that interseismic geodetic strain rates remain roughly stationary in the interim between large earthquakes and are representative of the tectonic loading which triggers them.
In this framework, the combination of earthquake catalogs (as complete as possible) with geodetically constrained rates of elastic strain buildup during the interseismic phase of the earthquake cycle has the potential to elevate strain rate into a spatial predictor for earthquake occurrence.As regards Greece and its surroundings, previous studies have examined either the consistency between geodetic and seismic moment rates (Chousianitis et al., 2015;Jenny et al., 2004;Sparacino et al., 2022) or compared smoothed seismicity rates with deformation models (Piña-Valdés et al., 2022).Nevertheless, the correlation between the spatial pattern of the regional seismicity and strain rates has not been explicitly investigated.Given that the strain rate mapping implemented in this study has high resolution and nation-wide coverage, we explore to which degree seismicity distribution correlates with our geodetic strain rate model.To do so, we use the earthquake catalog of the National Observatory of Athens (NOA) for the period from 1964 to December 2022 (available at https://www.gein.noa.gr/en/services-products/earthquake-catalogs/).This earthquake catalog within 19°-29°E and 34°-42°N contains (as of December 2022) more than 350,000 events and is the most complete catalog for the region of investigation.It's spatial and temporal characteristics regarding homogeneity and completeness have been extensively analyzed and documented in previous studies (Chouliaras, 2009;Chouliaras et al., 2013;Mignan & Chouliaras, 2014).We consider earthquakes with depths less or equal to 40 km and with magnitudes larger than 4.0 given that, since 1964, the magnitude of completeness is well below 3.8 (Mignan & Chouliaras, 2014;Figure S6 in Supporting Information S1).
We declustered the catalog following the method of Reasenberg (1985) so as to exclude clustered earthquakes (i.e., such as those occurring in aftershock sequences and swarms) and retain only independent events (the obtained declustered catalog is listed in Table S2).This decision is driven by the rationale of the study, that is, to examine if strain rates can potentially be used as a proxy for identifying the location and magnitude of future earthquakes; hence, the incorporation of foreshocks and aftershocks might bias the results.Figure S7 in Supporting Information S1 shows an example of Reasenberg's declustering method for the period after 2008, the year when the development of the Hellenic Unified Seismological Network (HUSN) started to take place and the density of seismic stations used to determine earthquake location increased significantly compared to the preceding period.It is evident in Figure S7 in Supporting Information S1 that the curve showing the cumulative number of earthquakes with time, includes several discontinuities, the largest of which are related to strong earthquakes with intensive aftershock series.On the contrary, the respective curve which has been produced using the declustered catalog displays a linear behavior which is indicative of the removal of clustered earthquakes.
At first glance, the distribution of earthquakes with M ≥ 4.0 seems to show a positive spatial correlation with current strain rates given that a considerable amount of earthquakes occur in high strain rate areas (Figure 5).The events presented in Figure 5 occurred after 2008, when, as stated previously, a significant network expansion associated with the formation of the HUSN resulted in improved quality of the NOA seismic catalog, better detectability of small events and shifting toward lower completeness magnitudes (Mignan & Chouliaras, 2014).Yet, to link the recorded seismicity with the strain rate model more efficiently, we calculate for each grid cell the geodetic potency rate as in Zeng et al. (2018) along with the corresponding earthquake count.We define, however, potency rates on the basis of the second invariant of the strain rate tensor instead of the maximum shear strain rates because stick-slip earthquakes in Greece occur at all kinds of tectonic environments (extensional, contractional and transcurrent).Hence, a quantity that describes the total strain rate such as the second invariant appears to be more appropriate than using dilatation rates or maximum shear strain rates which are associated with volumetric changes and shear deformation, respectively.
Next, we sorted the grid cells in order of decreasing second invariant of strain rate values and plotted the cumulative geodetic potency rates and the cumulative earthquake counts per cell, after normalizing both quantities to a maximum value of one (Figure 6).The cumulative earthquake count curves were calculated for various cutoff magnitudes, because we intend to examine not only whether current strain rates are spatially correlated with seismicity, but also to investigate if there is a specific magnitude above which such correlation becomes maximal and strain rates acquire the potential to serve as a proxy of where future earthquakes may occur.In the diagrams of Figure 6 the cumulative geodetic potency rate curve is illustrated with blue color and remains the same as the cutoff magnitude changes, while the earthquake count curve is depicted in red color and its shape depends on the respective cutoff magnitude.The light blue diagonal straight line that passes through the points (0, 0) and (N, 1), where N is the total number of grid cells, would be expected if earthquakes occurred randomly and were completely uncorrelated with the strain rate data.The Poisson curve, which is depicted with the dashed black line, tests the assumption that the events occur randomly in space via a numerical experiment.More precisely, we estimate the rate λ of a Poisson distribution as the ratio of the total number of events divided by the number of equal-area cells.This gives the average number of events per cell.We then generate for each cell a random number of events by drawing random deviates from the Poisson distribution with rate λ.Next, we calculate the cumulative number of events (i.e., as we count more cells, we accumulate the number of simulated events in these cells), and we normalize by dividing with the total number of events generated in all the cells.The normalized cumulative number of events per cell thus obtained is the Poisson benchmark in the case of random distribution of events per cell.Ideally, it should coincide with the theoretical no-correlation curve (light blue straight line).As evidenced in all plots of Figure 6, the random and Poisson lines are quite close.The fluctuations of the Poisson lines for the larger magnitudes are expected since the total number of events in these cases is small, and the fluctuations of the realizations drawn from the Poisson distribution become evident.Overall, the simulation-based Poisson curves confirm that the theoretical straight line in the case of randomly occurring events is valid.
The extent to which current strain rates are correlated with seismicity can be assessed by comparing the earthquake count and the geodetic potency rate curves of Figure 6.If the former appears well below the latter curve, the seismicity occurs in a random and irregular manner with respect to the strain rates, and earthquakes occur indistinguishably within low and high strain rate areas.In contrast, when the earthquake count curve matches or exceeds the geodetic potency rate curve, the seismicity distribution correlates well with strain rates, and earthquakes predominantly occur in fast-straining regions, indicating that strain rates could be potentially used in earthquake likelihood and seismic hazard models.From Figure 6 it becomes evident that as we move from lower to higher cutoff magnitudes, the correlation between strain rates and seismicity becomes stronger.It is also noticeable that above a certain cutoff magnitude the bulk of the earthquakes are spatially well correlated with higher strain rates.This is because the majority of the larger-magnitude events occur in the left part of the horizontal axis, where also the grid cells with the higher strain rate values extend, contributing to the fast steepening of the earthquake count curve.
To quantify the degree of correlation for each cutoff magnitude, we calculate the area that lies between the geodetic potency rate and the earthquake count curves taking into account that at a given point the latter can be Figure 6.Diagrams showing normalized cumulative geodetic potency rates (blue curve) and normalized cumulative earthquake counts (red curve) for various cutoff magnitudes reported on top of each plot.Both curves were obtained after sorting the grid cells in order of decreasing strain rates (i.e., the largest strain rates occur to the left of the x-axis).The light blue diagonal straight line would be expected if earthquakes occur in a completely random manner and the dashed black line checks this assumption through a numerical experiment based on the Poisson distribution.The difference between the earthquake count and geodetic potency rate curves is illustrated with the green shaded region whose area is quantified by the dA M metric (cf.Equation 1) which is reported within each plot.
either below or above the former.Accordingly, we quantify the difference between the two curves for a given cutoff magnitude M using the following relation: where the variable x ∈ [0,N] denotes the number of cells, and C M (x),S(x) ∈ [0,1] denote respectively the normalized cumulative earthquake counts and geodetic potency rates per cell.Since the functions S(x) and C M (x) are only known at a number of discrete points x i ,i = 0, 1, …,N the integrals are evaluated by means of the trapezoidal rule.If S(x) = C M (x) for all x, then dA M = 0. Positive (or negative) values of dA M indicate that the geodetic potency rate curve lies above (or below) the earthquake count curve.It is also possible that the relation between S(x) and C M (x) changes sign for some x i values.The sign of the integral dA M thus depends on the overall (signed) difference between the two curves.
The cumulative plots in Figure 6 reveal a consistent pattern, according to which the earthquake count curve is initially (for low cutoff magnitudes) positioned well below the geodetic potency rate curve and as we move toward larger magnitudes it progressively rises until it starts to locally exceed the geodetic potency rate curve.For lower magnitudes and particularly for earthquakes of M ≤ 5.0, this pattern suggests minimal spatial dependence of seismicity with strain rates, highlighting that earthquakes within this magnitude range exhibit a diffused spatial behavior and occur in both high and low strain rate areas.Interestingly, for larger magnitudes the earthquake count curve increases significantly faster, and it quickly approaches and exceeds the geodetic potency rate curve, implying that the spatial distribution of earthquakes correlates well with regions characterized by higher strain rates.This is appropriately quantified by the dA M metric which decreases monotonically as magnitude increases (Figure S8 in Supporting Information S1).More precisely, dA M starts with positive values which indicate that the largest part of the shaded (green) area corresponding to the inter-curve difference in Figure 6 lies beneath the geodetic potency rate curve; dA M becomes zero at M = 5.6 and then it becomes increasingly negative signifying that the majority of the green area lies above the geodetic potency rate curve.This implies that M = 5.6 marks a threshold magnitude above which correlation becomes significant and progressively stronger, suggesting that strain rates constitute a potential indicator for the location of future earthquakes with magnitudes above this threshold value.It is indicative that about 90% of the earthquakes of M ≥ 5.6, for each cutoff magnitude, occurred within the first 50% of the sorted grid cells, which correspond to areas where the second invariant of strain rates exceeds 40 nstrain/yr (see in Figure 5 the 40 nstrain/yr contour and the corresponding area which is illustrated with yellow and orange tones).Finally, so as to evaluate the impact of a different cell size, we repeated the aforementioned analysis using a finer and a coarser grid.In the first case we decreased the cell size to half and we set it equal to 0.05°× 0.05°, while in the second case we doubled the cell size and set it equal to 0.2°× 0.2°.The variation of the dA M metric with magnitude for both cases is illustrated in Figure S9 in Supporting Information S1 where it can be observed that the results are consistent with those in Figure S8 in Supporting Information S1 given that the highest correlation is achieved after M = 5.6.

Conclusions
In this work we produced a new strain rate and rotation rate model for Greece using a highly dense GPS velocity field.The improved data coverage, especially in formerly under-sampled regions, allowed us to provide a more refined view of the rates and spatial distribution of various velocity gradient tensor quantities and highlight some deformational patterns that, so far, were either poorly resolved or unidentified.We found the magnitude of the strain rate field, as it is quantified by the second invariant, to be highest across two well-defined provinces.The first follows the NAF and its two branches within the north Aegean, crosses central Greece and through the Gulf of Corinth it ultimately terminates in western Greece.The second encompasses the extensional province of western Turkey and the eastern Aegean Sea islands.These two high-strain provinces are in sharp contrast with the low-strain areas of northern Greece and central/southern Aegean Sea, where the rates which are accommodated are well below 20 nstrain/yr.
Across the southern part of the three southern Balkan countries that border Greece, we observed a smooth rotation of the extension axis from NW-SE in Albania to NE-SW in Bulgaria.In northern Greece and southern Bulgaria we document slow counterclockwise rotation whose peak (1°/Myr) correlates with the NW-SE trending Strymon fault system.As we move southwards, the high-strain province that extends all the way from the NAF in the east to the Ionian islands in the west, remarkably reproduces the westward propagation of the NAF into north Aegean Sea and the separation of the Aegean microplate from Eurasia.This zone which encompasses from east to west the NAF, the NAT, the NAB, the extensional province of central Greece, the Gulf of Corinth and the large transcurrent structures of western Greece, accommodates the largest amounts of strain rates and rotations rates within the study area.Strong shearing accompanied by fast clockwise rotation is acting across the NAF and the NAT.
The NAF accommodates >100 nstrain/yr of shear across the Sea of Marmara and within the Gulf of Saros, while the NAT accommodates 80-100 nstrain/yr.Also, these major interconnected dextral shear zones induce a narrow band of clockwise rotation with rates that range from 3°/Myr to 6°/Myr.Further to the SW, increased shear strain rates (60-80 nstrain/yr) are maintained throughout the entire Skyros Basin and within the eastern part of the NAB approximately up to 23.8°E.West of this longitude, however, shearing within the NAB decreases and the increased dilatational component of the strain rate tensor with rates of ∼80 nstrain/yr imply the prevalence of extensional deformation.This pattern highlights the different character of these two basins which both host normal and strike-slip faults, as well as marks the onset of the ∼N-S oriented extension that dominates central Greece.
The extensional province of central Greece is characterized by N-S to NNE-SSW extension accompanied by notable (∼4°/Myr) clockwise rotation and forms a transition zone between the quasi-parallel transcurrent structures of north Aegean and western Greece.The largest extensional rates occur at the western Gulf of Corinth, where they reach 270 nstrain/yr, while across the eastern Gulf of Corinth the principal extensional axes experience a ∼50% rate decrease becoming ∼130 nstrain/yr.To the north of the Corinth Rift, the widespread active extension throughout central Greece appears to be more complex than a straightforward distributed N-S extension, as it has been considered so far in geodetic studies.In this context, across this region which is located between the strongly extensional system of the Corinth Rift and the transtensional regime of the NAB, the accommodated dilatation rates vary spatially from 60 nstrain/yr to 100 nstrain/yr and in some places these rates are coupled with sizable amounts of shear strain rates.At these locations, which include the north Evia Gulf rift and the Pagasitikos Gulf, this strain pattern implies the existence of a significant amount of oblique extension and the occurrence of earthquakes with large strike-slip component to the SW of the NAB cannot be considered unexpected.
To the northwest of the western tip of the Gulf of Corinth, we document a very clear zone along the NNWtrending KFS that accommodates 100-240 nstrain/yr of shear strain rates and undergoes counterclockwise rotation with rates up to 4°/Myr.This configuration is consistent with pure left-lateral shear across this fault zone which is situated at the intersection of the clockwise rotating domains of central and western Greece.Accordingly, we conclude that the KFS has been formed in response to the movement of two same-sense rotating domains whose motion induced the opposite sense of rotation across their intersection.Apart from the area in continental western Greece that hosts the aforementioned tectonic boundary along the KFS, high shear strain rates exceeding 60 nstrain/yr are maintained throughout the central and southern Ionian Sea and across the NW part of Peloponnese where they continue as far as the NE-SW trending KAFZ.Apart from wide-scale shearing, this domain experiences fast clockwise rotation with rates of 5-8°/Myr.In the region of the central Ionian islands, this rotation is accompanied by internal contraction and not extension as in the case of central Greece, something which emphasizes the different nature of their driving mechanisms.Negative dilatation rates were also resolved in northwestern Greece where a ENE-WSW trending compressive zone of ∼90 nstrain/yr can be interpreted as the footprint of the active collision of the Apulian continental lithosphere with northwestern Greece.
The rotation rate field showed that the southern branch of the NAF comprises the northern limit of the broad counterclockwise rotation that characterizes the extensional province of western Turkey and the eastern Aegean Sea islands near the Turkish coastline.Across this region, the mapping of the second invariant of strain rates highlighted that a NNW-SSE oriented wide zone of high-magnitude strain separates the Aegean microplate from the Anatolian plate.This zone which is situated between the westward moving Anatolian plate and the faster and SW moving Aegean microplate accommodates the curved stretching of the entire plate system as it moves in response to the fast Nubia-Aegean convergence.We suggest that this wide high-strain province could represent a "soft" boundary between Aegean and Anatolia under the concept that a number of smaller faults, instead of one distinct and well-defined tectonic boundary, accommodate most of the region's deformation which is induced by the aforementioned kinematic framework.Regarding the forearc of the Hellenic subduction system we showed that in both edges the dominant mode of crustal strain is E-W extension and as a result the faults that accommodate the deformation there are the N-S trending normal faults.
Finally, we explored to which degree the GPS strain rates correlate with seismicity.Using our model for the distribution of the second invariant of strain rates and a declustered earthquake catalog of the regional seismicity, we showed that the vast majority of earthquakes with M ≥ 5.6 are located at high strain rate areas where the second invariant exceeds 40 nstrain/yr.On the contrary, lower magnitude earthquakes and particularly those with M ≤ 5.0 appear uncorrelated with strain rates and occur in both high and low strain rate areas exhibiting a diffuse spatial behavior.The strong positive correlation that we demonstrated for events with M ≥ 5.6 verifies a strong spatial dependence of strain rates with seismicity above this threshold magnitude.This implies that geodetic strain rate data comprise a suitable proxy for the prediction of the location of future moderate and large earthquakes above the aforementioned magnitude, as well as justifies their potential to become a useful complement in seismic hazard assessment models.

Figure 2 .
Figure 2. Horizontal GPS velocity vectors with respect to stable Eurasia.Uncertainty ellipses correspond to the 95 per cent confidence level.

Figure 5 .
Figure 5. Second invariant of strain rates superimposed by earthquakes with M ≥ 4.0 and hypocentral depths up to 40 km for the time period 2008-2022.Seismicity located outside the strain rate grid cells is not shown.

Table 1
Individual GPS Velocity Fields That Were Aligned Into a Common Reference Frame