3-D Body-Wave Velocity Structure of the Southern Aegean, Greece

1. Introduction

The Southern Aegean, situated in Greece, boasts the southernmost part of the Hellenic arc, which commences in the area W-SW of Peloponnese and runs southward, dividing into Pliny, Ptolemy, and Strabo trenches. The arc culminates to the east of Rhodes island, branching into Anatolia and South of Cyprus [1,2]. This region is prone to frequent seismic activity, owing to the convergence of the Eurasian and African tectonic plates, as substantiated by significant historical seismic events [3]. This movement is caused by the subduction of the oceanic Mediterranean lithosphere under the Aegean continental one, with a speed of roughly 35 mm/year [4-10].

This study delves into the Southern Aegean region (Figure 1), which has a significant portion of the overall recorded intermediate-depth seismic activity. This particular area is considered one of the most active ones in the Eastern Mediterranean Sea due to intense tectonic movements during the Late Quaternary period, resulting in significant uplifts in local areas and the creation of extensional faults that run in roughly two directions: WNW-ESE and NNE-SSW [4,11-13]. In the Cretan Sea, just north of Crete, there are three primary tectonic depressions between the Southern Cyclades and Crete, Karpathos, and Rhodes. The main tectonic structures in the region are oriented in a NW-SE direction towards the west and a NE-SW direction towards the east, showcasing the curvature of the Hellenic Arc. The broader area of Crete has a complex geodynamic structure where transpressional tectonics dominate the southern part of Crete, in contrast to normal and strike-slip forms on the shallow part of the crust, resulting in a heterogeneous stress field [14,15]. In the northern area, between Santorini and Crete, there is an observation of extension in an E-W direction, while in the western part of Crete, N-forward shortening was seen [16] (Figure 2).

It's interesting to note that the seismic activity in Crete and the neighbouring areas is mainly due to the geodynamics of the Aegean and Mediterranean plate convergence, with events of shallow depth (<25 km) characterizing it. However, the lack of uniform and optimal coverage with seismic stations on the island and the southern Aegean Sea and absence of local velocity models are the main disadvantages in highlighting the earthquake distribution in areas near or along the Hellenic Trench (HT). Fortunately, recent seismic sequences have led to an improvement in the local station coverage, enriching the available catalogues. With improved 1-D regional models and well-localized aftershocks recorded by dense seismic monitoring networks, the region can now provide a basis and baseline for more constrained tomographic studies. The recorded data improved the knowledge of the seismotectonic structures on the island and helped to provide a better 3-D body-wave velocity model through seismic tomography, examining the correlation between seismicity and the activation of the local fault pattern.

Seismicity studies since 1990 have revealed active Quaternary faults, major compressional structures (i.e. thrusts), and other geologic features in the back-arc basin north of Crete Island. The planes for the available focal mechanism solutions of offshore events near the Hellenic Trench south of Crete (Figure 3) show a direction of compression to be nearly NNE-SSW, while the onshore ones appear to have an almost ESE-WNW to N-S extension, which is associated with fractures along the main mapped faults [3,12,17-19]. Historical records indicate that three large shallow earthquakes occurred in this region, with the first one occurring southwest of Crete in 365 AD, the second one occurring near Rhodes in 1303, and the third one taking place on February 16, 1810, in the area between Crete and Rhodes [19,20]. The existence of permanent stations of the regional Hellenic Unified Seismological Network [21], close to the most significant events, with the nearest stations being Knossos (KNSS), Pefkos (PFKS) (Institute of Physics of the Earth's Interior and Geohazards of the Hellenic Mediterranean University Research Center; IFEGG-HMURC), and temporary ones [22,23] installed after large earthquakes, makes a significant contribution to the depth accuracy of the first hypocentral solution during the routine analysis in near real-time (http://www.geophysics.geol.uoa.gr/stations/gmaps3/leaf_stations.php?map=2&lng=en).

In the area of Aegean Sea (Greece), there have been many studies conducted on regional scale over the past forty years. Seismic velocity models have been obtained for various areas including Western Greece, and Aegean Sea [24-27]. However, it's worth mentioning that the lack of permanent seismic stations in the islands of Southern Aegean has made imaging of the crustal structure more challenging. In this work, we provide new 3-D body-wave velocity models of the crust and uppermost mantle (down to 80 km) utilizing a new big collection of earthquake travel time data. These findings may aid in determining the relationships between surface geology and subsurface structures driven by complex geological processes.

2. Data and Method

The present study focuses on Crete and its neighboring regions in the Southern Aegean Sea. Several permanent stations of the Hellenic Unified Seismological Network [22,23,28-30] operate in this area, complemented by 4 temporary stations (CRE1–4) of the Geodynamics Institute of the National Observatory of Athens (GI-NOA), between September 27, 2021 and June 15, 2022, mainly contributing to the depth accuracy of the Arkalochori aftershock sequence for this period of time (Figure 4) [21,31-33]. The initial dataset comprises 1813 events from June 2018 to February 2023, with M_L≥3.2, recorded by at least 6 stations. Since 2007 continuous waveform recordings have been obtained from HUSN in real-time, comprising of both 1-component (1 Hz) (THR2-9; Santorini Island Seismic Network operated by the Aristotle University of Thessaloniki) and 3-component (either broad-band or short-period) stations, while manual arrival-time picking is applied to obtain hypocentral locations using the SeisComP3 graphical user interface [34] and a custom regional 1-D velocity model for the region of Southern Aegean sea [35]. In the first stage of the sequence analysis, hypocenters were located in near real-time, employing the Hypo71 single-event algorithm [36]. Then they were relocated running the HypoInverse code [37].

In this study, the inversion of body-wave travel-time data was based on LOTOS scheme as described by [38]. This code was previously used for studies of many different subduction zones with similar geometry of observations and area dimensions [39-43]. For the inversion procedure, 1265 manually revised events located on the island of Crete and its neighbouring regions were included during the first step. The dataset comprised 20587 P and 15370 S arrival-times, with at least 15 phases per earthquake (at least 10 for P- and 5 and for S-picks respectively) and ratio of S to P residual smaller than 1.5 (Figure 4). The inversion algorithm provides two different choices: a) inversion for both V_P and V_S (V_P–V_S scheme) using the respective travel-time residuals (dt_P and dt_S), and b) 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 the current study, tomographic inversion was applied for both schemes in order to obtain more information on the body-wave anomalies [38,44].

Sensitivity analysis for the available dataset was performed by applying the checkerboard synthetic test. This method uses alternating anomalies of positive and negative velocity perturbations on an initial 1-D gradient model evenly distributed in a checkerboard pattern throughout the model. The stability of the inversion is controlled by amplitude damping and flattening regularization. The inversion of the large sparse matrix was performed using Least Squares with QR decomposition (LSQR) while the optimal values of P- and S-wave amplitude damping, and smoothing as well as the station, source coordinates and origin times corrections were assessed based on the results of synthetic modeling as it is described in [38, 45-46].

3. Results

The available tomographic software necessitates as input data the longitude and latitude of the available network stations and the arrival times from the recorded seismicity. The coordinates of the hypocenter and the time of origin time are optional, as they are determined during the implementation of calculations. However, if preliminary hypocentral locations are available, they can be utilized to reduce the processing time of the operations. In addition, a generic 1-D velocity model and a set of input parameters for performing the convergence iteration steps, comprising parameterization, grid dimensions, and damping parameters, were considered [38]. As 1-D starting velocity model was set the one derived by [35]. This model was then optimized by running the complete tomographic procedure several times by calculating the average velocities at several depths and using them as a new starting model in the next phase of complete tomographic inversion (Figure 5).

The model parameterization of the velocity field should be able to delineate, according to the local characteristics, the shape and perspective of heterogeneities. A nodal representation was used, since the velocity field reconstructed by a three-dimensional grid does not adopt any specific geometry of heterogeneities [47]. The interval of the grid nodes was kept significantly smaller than the expected size of the anomalies (<20 km) to weaken the distortion of the resulting models due to the grid configuration (Figure 6).

The optimal grid mesh has been determined considering the stations/events geometry. In addition, to further reduce the impact of the model parameterization on the solutions, the inversion was repeated utilizing several grid orientations (e.g. 0°, 22°, 45°, and 67°). The results were aggregated via a 3-D model of the absolute P- and S-wave velocities by simple averaging. Examples of node distributions and ray densities for two of the grid orientations in different depths are presented in Figure 6.

In the vertical direction, the grid spacing was inversely proportional to the ray density, but could not be smaller than a predefined value (10 km in our case; Figure 7).

Figure 8 shows the difference between the initial and the final location of the events derived by LOTOS. In Figure 9 (left) the P- and S-wave travel-time residuals histogram before and after earthquake relocation is shown. The initial P-wave travel time residuals were mainly distributed in the range of −2 to 1.8 s, and the average travel time residuals were 0.46 s, while the P-wave travel time residuals were mainly distributed in the range of 0.6 s after relocation, and the average residuals were reduced to 0.27 s. Figure 9 (right) shows that the residuals of the initial S-wave travel time were also significantly reduced from 0.57 s to 0.28 s. Figure 9 (right) shows that the travel time residuals after relocation are mainly distributed in the range of −2.2 to 2 s. As seen in Figure 8, there is a clear difference between the initial and final locations of the events obtained by LOTOS. Moving on to Figure 9, we see a histogram of the P- and S-wave travel-time residuals before and after the earthquake relocation. The initial P-wave travel time residuals (left) were mainly distributed in the range of −2 to 1.8 s, and the average travel time residuals were 0.46 s. After relocation, the P-wave travel time residuals were mainly distributed in the range of ±0.6 s, and the average residuals were reduced to 0.27 s. In the right side of Figure 9, we see that the residuals of the initial S-wave travel time were also significantly reduced from 0.57 s to 0.28 s (Table 1). The overall seismic travel time residuals after relocation are mainly distributed in the range of −2.2 to 2 s. It appears that there has been a redistribution of the depths of the earthquake foci, with a majority of them now situated within a range of 0-25 km depth. However, it is worth noting that there are still significant concentrations of earthquake foci located below 50 km depth near the volcanic centers of the SAVA.

3.1. Resolution tests

To define the fidelity area in this study, we used the checkerboard method [48] in order to indicate the resolution and errors associated with the inversion. This method uses alternating anomalies of positive (fast) and negative (slow) velocity perturbations evenly spaced throughout the model in a checkerboard pattern. Generally, ray-path distribution, model parameterization, and smoothing determine data resolution [49]. Among the limits imposed of such a test is the fact that the extent of near-vertical alteration of structure is difficult to be fully recognized due to a structure (checkerboard) where the diagonal elements in the vertical plane are strongly associated with dominant ray directions [50].

The error in source mislocations for both vertical profiles was reduced after five iterations of simultaneous inversions of the source and velocity parameters (Figure 10). It is worth noting that the uncertainty of the source depth locations was greater than that of the horizontal coordinates.

Following, a rigorous procedure was performed, ensuring that the synthetic models (checkerboard) accurately reflect the real data processing. The size of cells in the initial models corresponded to the expected anomalies, and we defined spiked regions with ±10% variability in velocity structure compared to the reference 1-D velocity model. To add a realistic element in the procedure, we added noise of 0.27 s for P- and 0.28 s for S-wave arrival-time picks, which was derived from real-data inversion. The travel-times for the paths between the source and receiver were then computed, using 3-D ray tracing that followed the bending algorithm principles, mirroring the real observation system. A test of various parameters was performed and we found the ones that gave us the most accurate results in reconstructing the model of checkerboard anomalies. To control the accuracy of the inversion, we considered errors such as mis-picks, mis-locations, and inaccurately of the determined ray-paths.

The adopted procedure included three different sets of dimensions of anomalies for the horizontal tests (20x20 km² , 40x40 km² and 80x80 km²), in order to define the limitations of our model. The variations (%) of body-wave velocity anomalies (±10%) and the V_P/V_S ratio distribution are presented in Figure 11, Figure 12 and Figure 13, at depths of 10, 40 and 80 km. Based on our findings, it appears that the first set of anomalies with a 20x20 km² grid size is either not resolved or poorly resolved at various depth slices, whereas the larger sets with a 40x40 km² and 80x80 km² cell size provide reasonable resolution throughout most of the study area (in between Pliny Trench and the central-eastern part of SAVA). The 40x40 km² size pattern was particularly robustly resolved in most parts of the study area at shallower depths, up to around 40 km. However, at a depth of approximately 60 km, the 80x80 km² anomalies were more easily resolved than the smaller ones. It seems that the models with smaller anomalies were not as clearly distinguished at deeper depth slices.

It's worth noting that the anomalies of 20 km board size were not well-resolved for either the P- or S-wave velocity models, while the respective ones of 40 and 80 km were reasonably restored in all depth intervals (10-80 km), in the major part of the study area (34.75°-35.9°N, 23.0°-26.1°E). The results of the checkerboard tests suggest that there is a horizontal smearing towards the NE part of the study area, likely due to a lack of seismological stations and sparse epicentral distribution.

To estimate the vertical resolution, we performed a synthetic test with anomalies defined along vertical profiles (Figure 14). Based on the patterns that we needed to investigate, we began setting the size of the checkerboard anomalies. We began with a board of 20 km in the horizontal direction while in the vertical one, the signs of the anomalies changed at depths of 20, 40, 60 and 80 km. It seems that for the shallower part of the model (~60 km), the vertical resolution was sufficient to detect the change in the anomaly sign at a depth of 20, 40 and 60 km. However, the model did not have sufficient resolution to recognize anomaly sign changes below 80 km. It seems that we have some horizontal smearing towards the edges of the study area (SW and NE), potentially due to a lack of seismological stations and reliable hypocentral solutions. It's interesting to note that the transition between negative and positive velocity anomalies is reconstructed sufficiently for the main part of the study area in profile 1 (west), with some horizontal smearing along the ray paths towards the NE. Profile 2 (east) appears to be less resolved than profile 1, particularly to the northern part of the section and for depths below 60 km. It's possible that the lower vertical resolution than the respective horizontal one may be due to source and velocity parameters trade off, as suggested by previous studies [51]. The synthetic tests we conducted showed that the absolute amplitudes of the body-wave anomalies were up to 4-5% smaller than the respective ones of the starting checkerboard grid. Overall, the checkerboard model was recovered to a satisfactory level, despite the presence of errors and outliers.

3.2. Tomographic inversion results

The generated P-, S-velocity anomalies (%), and the V_P/V_S ratio distribution are depicted in four horizontal slices at 10, 25, 40 and 80 km depth (Figure 15), and two cross-sections (profile 1 and 2) perpendicular to the western and eastern part of the HT (Figure 16). The evaluation of the acquired values in both horizontal slices and vertical profiles is restricted to model portions with a sufficient reconstruction of the checkerboard (Figure 14). For P- and S-waves, velocity perturbations in the shallow depth layers range between -10 and +10%, while the V_P/V_S ratio varies between 1.66 and 1.98. According to synthetic modeling tests, the velocity structure is resolvable for depths ranging from 25 to 60 km and fairly reconstructed at 80 km depth for the central part of the area of study (34.75º-36.15ºN, 23.2º-26.25ºE). This fact, together with intense ray coverage down to a depth of 40 km (Figure 6), lends confidence to the interpretation of the 3-D inversion results.

It's interesting to note the prominent features identified by the present tomographic model, particularly the continuous fast velocity anomalies located north of Crete. This pattern was also identified in previous studies [25,52,53,54], lending credibility to the model's accuracy. These high velocities are mainly attributed to be the subduction of the Tethys oceanic lithosphere beneath the Aegean, a conclusion that is consistent with previous research. Additionally, it's worth noting that the front of the Hellenic Subduction Zone (HSZ) slab appears to be parallel to the HT in the region between Chania and Antikythera. Finally, it's intriguing that the high velocity anomalies are distributed in an amphitheatre shape that is almost parallel to the HSZ, with the area of the positive anomalies increasing with depth.

Based on the results of this study, there are some noteworthy discontinuities, highlighted by the slow and fast body-waves (P,S) velocity anomalies towards the western, central-western region of Crete, specifically within shallow depths (10-25 km). It is plausible that this pattern is correlated with the local neotectonic faults striking in that direction, as illustrated in Figure 15. These faults divide the Mesozoic Pindos Unit basement from the Neogene and Quaternary deposits. A similar pattern of body-wave velocity (P, S) anomalies is observed in the same depth slices (10 and 25 km) near the town of Chania, with slow and fast deviations from the reference 1-D model to the west and east of Kissamos gulf. This shallow, broad, slow V_P anomaly down to 25 km depth characterizes the collapsed zones.

4. Discussion

The derived anomalies of body waves (dV_P, dV_S) and, V_P/V_S ratio provided important information about the Southern Aegean regional tectonics, and secondarily active faults of smaller scale (>20 km). In the study area, the seismic tomography model revealed a predominance of high-velocity anomalies North of Crete, in the region between Kythera and Karpathos, while close to the SAVA, the region is marked by significant low-velocity anomalies in the crust and uppermost mantle, beneath the active arc volcanoes. The seismicity related to the HSZ is connected to high-velocity anomalies in the Sea of Crete and the islands south of SAVA. Its extension reaches up to 120-150 km (Figure 16), while the observed termination of the seismic activity in even greater depths may be due to plate dehydration embrittlement or runaway thermal shear stress, as highlighted by studies in similar environments [57].

More specifically, the main anomalies that have been identified in the region are (a) a large zone of fast velocity anomalies (dV_P, dV_S), and low V_P/V_S ratio in the shallow part of the Hellenic island-arc area (Kythera-Crete-Karpathos, Figure 17) (b) a region of slow (negative) anomalies in the SAVA that may reflect hot upwelling flow related to slab dehydration, causing arc magmatism and volcanism in this region [44,46,55,56] (c) an area of high V_P/V_S ratio (~1.85) north of the Sea of Crete, in the SE part of Peloponnese and its offshore area in the depth slices of 10 and 25 km. [26] have related this part to fluid-saturated sediments, non-consolidated soils, and/or cracked subducted crust. The results in profile 1 (Figure 16) for depths ranging from 10 to 30 km may indicate a possible upward movement of fluids connected to this mechanism. However, it's important to keep in mind that the limited vertical resolution provided by the synthetic tests makes it challenging to give a definite interpretation of these results. In Figure 17 we mark the upper limit of the western segment of the HSZ slab. In the SW-NE direction, the slab’s dip decreases from ~60° (<40 km) to ~30-40° in depths greater than 70 km. The subduction depth is consistent with the results of [58].

In shallow depths (10-25 km), near the island of Crete, some anomalies of smaller size have been noticed. The contrast between fast and slow anomalies is mainly linked to the local fault pattern, which produces a comparable distribution of V_P/V_S ratio (Figure 16; Matala, Ptolemy Trench: 34.75º-35.15ºN, 25.25º-25.75ºE). This anomaly distribution is mainly affected by local surface geology at these depths [44,45]. At shallow depth (10 km), close to the city of Heraklion, a discontinuity of fast to the east and slow to the west of V_P anomalies is discovered, possibly related to the 2021-2022 earthquake sequence of Arkalochori [31-33,59] The formed anticorrelation between slow P- and fast S-wave velocity perturbations form low values of the V_P/V_S ratio (~1.70-1.73). This result supports the high-density fault pattern in this area as tectonic studies have shown [9,18] and the thick Neogene to Quaternary sediments [31,60].

5. Conclusions

Based on the travel times of the manually revised events from 2018 to 2023, we obtained a new high-resolution 3-D V_P model of the crust and upper mantle on regional scale for the region of the Southern Aegean. The resulting 3-D velocity models aimed to the improvement of the hypocentral parameters of the seismic events that occurred in this area. This process served as a springboard to a more detailed image of velocity anomalies that were obtained in the following stage of the tomographic inversion.

Both P- and S-wave velocity anomalies, as well as the distribution V_P/V_S ratio, support the existence of a thin crust west of the island, in contrast to a thicker one under Crete, which constitutes the accretionary prism in this portion of the Eastern Mediterranean Sea.

The regional distribution of seismic activity, and body-wave velocity anomalies (P, S), and V_P/V_S ratio values revealed:

• A complex shallow (<10 km) structure in Crete's central region mainly attributed to the dense pattern of neotectonic faults;
• A region of significant low-velocity anomalies in the crust and uppermost mantle, close to the SAVA, marked by the active arc volcanoes;
• A smoother and more continuous image in deeper slices (>40km) where high V_P/V_S ratio (>1.82) distribution is observed;
• The existence of a low-angle positive anomaly of V_P correlated to the observed intermediate-depth seismicity (H≥40 km) in this part of the study area. This result could be related to the diving HSZ slab.

It seems that additional research is necessary to fully understand the Seismic Hazard of Crete. Additionally, expanding the research areas will be useful in investigating the geometry of the Wadati-Benioff zone in the Southern Aegean, the properties of the seismogenic layer, and comparing this data with past studies and other regions globally. It may also be worth exploring the possibility of a tear in the slab near Karpathos and Rhodes Island in the eastern portion of the HSZ.