TY - JOUR AU - Papadimitriou,, P AB - SUMMARY This research provides new constraints on the intermediate depth upper-mantle structure of the Hellenic lithosphere using a three-step Rayleigh-wave tomography. Broadband waveforms of about 1000 teleseismic events, recorded by ∼200 permanent broad-band stations between 2010 and 2018 were acquired and processed. Through a multichannel cross-correlation technique, the fundamental mode Rayleigh-wave phase-velocity dispersion curves in the period range 30–90 s were derived. The phase-velocities were inverted and a 3-D shear velocity model was obtained down to the depth of 140 km. The applied method has provided 3-D constraints on large-scale characteristics of the lithosphere and the upper mantle of the Hellenic region. Highlighted resolved features include the continental and oceanic subducting slabs in the region, the result of convergence between Adria and Africa plates with the Aegean. The boundary between the oceanic and continental subduction is suggested to exist along a trench-perpendicular line that connects NW Peloponnese with N. Euboea, bridging the Hellenic Trench with the North Aegean Trough. No clear evidence for trench-perpendicular vertical slab tearing was resolved along the western part of Hellenic Subduction Zone; however, subcrustal seismicity observed along the inferred continental–oceanic subduction boundary indicates that such an implication should not be excluded. The 3-D shear velocity model supports an N–S vertical slab tear beneath SW Anatolia that justifies deepening, increase of dip and change of dip direction of the Wadati-Benioff Zone. Low velocities found at depths <50 km beneath the island and the backarc, interrelated with recent/remnant volcanism in the Aegean and W. Anatolia, are explained by convection from a shallow asthenosphere. Structure of the Earth, Tomography, Surface waves and free oscillations, Dynamics of lithosphere and mantle, Subduction zone processes 1 INTRODUCTION In this work, the intermediate depth seismic structure of the Hellenic lithosphere and its relation to seismicity is studied. The Hellenic lithosphere includes the Hellenic Arc and trench system, as well as the backarc Aegean microplate (e.g. Cianetti et al. 2001). This region is dominated by the presence of a complex system of lithospheres, that is Adriatic, Eurasian, Aegean and African plates. Greece is one of the most rapidly deforming parts of the continents globally, capable of generating earthquakes even larger than M∼8.0 (e.g. Papazachos & Papazachou 1997). Complexity of the Hellenic lithospheric system results from slab retreat of the N–NE subducting African Plate along the Hellenic Subduction Zone (HSZ) and subsequent accretion of its surficial crustal sediments, resulting in the closure of the Tethys oceanic domains (Dercourt et al. 1993; Channell & Kozur 1997; Robertson 2004). Extensive research engaging geological, seismological and geodetic approaches has explained the geodynamics of the Hellenic lithosphere as a result of the northward advancement of the African and Ionian crust since the Cretaceous (e.g. Jolivet & Brun 2010). Across this region, Africa and Eurasia converge at a rate of ∼25–40 mm yr−1 (McClusky et al. 2000; Reilinger et al. 2006, 2010) forced by the subduction of the Tethys oceanic slab under the Aegean microcontinent, resulting in roll-back of the subduction zone. Overall SW propagation dominates the regional kinematics, explained as the result of three fundamental processes related to convergence between the African, Aegean, and Anatolian plates (McClusky et al. 2000). The first process is the subduction of the African lithosphere and subsequent accretion of sediments. The second process regards the westward propagation of the Anatolian Plate into the Aegean due to the northward push from the Arabic Plate. The third process is the extension of the overriding plate in the backarc domain. According to geodesy, the average horizontal velocity rate across the Greek region is about 10 mm yr−1, which becomes more significant along the Hellenic Arc (∼40 mm yr−1; McClusky et al. 2000; Reilinger et al. 2006), which is the active margin between Aegean and African lithospheres, marked by a 1300-km-long trench extending from Corfu to Rhodes, separated and offset by 100–150 km by the dextral Cephalonia Transform Fault (CTF, e.g. Louvari et al. 1999, Fig. 1). Figure 1. Open in new tabDownload slide Map of major tectonic features in the region of Greece and placenames discussed in the text. Selected seismicity (M ≥ 3.0) is from routine analysis at NKUA-SL (period 2012–2017) and from the catalogues of Bocchini et al. (2018). Fault lines from Karakonstantis & Papadimitriou (2010). CTF: Cephalonia-Transform-Fault, CLAB: Cephalonia-Lefkada-Akarnania Block, SAAVA: South Aegean Active Volcanic Arc, CHSZ: Central Hellenic Shear Zone, MFZ: Movri Fault Zone, HSZ: Hellenic Subduction Zone, NAT: North Aegean Trough, NAF: North Anatolian Fault, CG: Corinth Gulf. Figure 1. Open in new tabDownload slide Map of major tectonic features in the region of Greece and placenames discussed in the text. Selected seismicity (M ≥ 3.0) is from routine analysis at NKUA-SL (period 2012–2017) and from the catalogues of Bocchini et al. (2018). Fault lines from Karakonstantis & Papadimitriou (2010). CTF: Cephalonia-Transform-Fault, CLAB: Cephalonia-Lefkada-Akarnania Block, SAAVA: South Aegean Active Volcanic Arc, CHSZ: Central Hellenic Shear Zone, MFZ: Movri Fault Zone, HSZ: Hellenic Subduction Zone, NAT: North Aegean Trough, NAF: North Anatolian Fault, CG: Corinth Gulf. HSZ is the most prominent feature within the Hellenic lithosphere, widely recognized in all global and regional tomographic models. A commonly accepted configuration of HSZ concerns a high-velocity anomaly associated with lithosphere of continental and oceanic origin, extending along the concave part of the Hellenic Arc. The footprints of the slab, found much deeper than the lowest intermediate seismicity (∼180 km depth) of the Wadati-Benioff Zone (WBZ), penetrating below the upper-lower mantle discontinuity, has led to the assumption that the slab is continuous and involves alternating subduction of both oceanic and continental crusts along the HSZ (Faccenna et al. 2003; Meier et al. 2004; van Hinsbergen et al. 2005; Jolivet & Brun 2010). Several arguments have been proposed against this generic model, concerning the affinity and the vertical and horizontal continuity of the slab, which has been the subject of scientific debates during the last decades. Although there is abundant intermediate depth seismicity evidence for oceanic subduction at the southern part of HSZ, the lack of subcrustal seismicity in the northern part puts in question the relation of the therein high-velocities with a subduction process. Recent studies (Pearce et al. 2012; Halpaap et al. 2018) have imaged a subducted continental crust as a low-velocity layer beneath NW Greece, which, contrary to the southern oceanic crust (Suckale et al. 2009; Pearce et al. 2012), does not host subcrustal seismicity. Although these findings provided evidences on the affinity of the slab, the geometry of the transition between the two types of lithosphere has not yet been clarified and the question about the slab continuity remains. Moreover, the lower boundary of the lithospheric plates of the region, a very important parameter for understanding wide-scale tectonics, is still a poorly known quantity (Kind et al. 2015). In this work, we measure dispersion properties of long-period Rayleigh surface waves, which are essential for modeling 3-D upper mantle and asthenospheric shear wave velocity structures with particularly good lateral resolution (e.g. Laske & Masters 1996), towards shedding light to the above as well as other issues related to the geodynamic frame of the Hellenic lithosphere. To this purpose, we have exploited new data of teleseismic Rayleigh surface waves acquired from modern broadband deployments and applied a recently developed approach for surface wave tomography (Jin & Gaherty 2015). We take advantage of the capability of the regional Hellenic Unified Seismological Network (HUSN) of broadband receivers to roughly act as an array in order to correlate similar waveform recordings of teleseismic Rayleigh waves at neighbouring stations and perform phase-delay measurements. Thus, a significant increase of input data and consequent improvement of the 3-D shear velocity model is achieved with respect to the classical surface wave two-station method that has been previously implemented in the Hellenic region (Bourova et al. 2005; Kassaras et al. 2009; Salaün et al. 2012). The coherent results from multiple teleseismic events are stacked to produce phase velocity maps while the derived dispersion curves are then inverted to construct 3-D, vertically polarized shear wave velocity (VSV) models throughout the region. Finally, the produced tomographic images are combined with reliable hypocentral locations of subcrustal earthquakes registered in a newly compiled seismic catalogue. 2 SEISMIC STRUCTURE AND GEODYNAMIC SETTING The seismic structure of the Hellenic lithosphere has been constrained using various approaches of passive and active imaging. Passive imaging concerns teleseismic (e.g. Spakman et al. 1988, 1993; Bijwaard et al. 1998; Wortel & Spakman 2000; Faccenna et al. 2003; Piromallo & Morelli 2003; van Hinsbergen et al. 2005; Koulakov et al. 2009; Zhu et al. 2015) and local/regional tomography (Papazachos & Nolet 1997; Tiberi et al. 2000; Karakonstantis 2017; Halpaap et al. 2018; van der Meer et al. 2018; Hansen et al. 2019), receiver function analyses (e.g. Li et al. 2003; Gesret et al. 2011; Sodoudi et al. 2015; Sachpazi et al. 2016), inversion of scattered teleseismic waves (Suckale et al. 2009; Pearce et al. 2012), active reflection (e.g. von Huene et al. 1997; Finetti & Del Ben 2005; Kokinou et al. 2005, 2006; Zelt et al. 2005) and refraction studies (e.g. Bohnhoff et al. 2001). Surface waves have also been exploited to investigate the S-wave structure of the Aegean area (Bourova et al. 2005; Karagianni et al. 2005; Kassaras et al. 2005, 2009), with the most recent comprehensive broad-scale tomography by Salaün et al. (2012) not sufficiently sampling the western part of HSZ. 2.1 The West Hellenic Subduction Zone (WHSZ) WHSZ extends ∼400 km along the western coast of Greece, from Corfu to western Crete, where convergence is documented, associated with the subduction of multiple ocean basins and continental units, from the Jurassic–Cretaceous until today (e.g. van Hinsbergen et al. 2005). The rate of subduction varies along the WHSZ, ranging from ∼8 mm yr−1 beneath NW Greece and ∼40 mm yr−1 beneath SW Greece (e.g. McClusky et al. 2000). A less dense continental lithosphere, associated with the Adriatic microplate, which has subducted north of ∼N38° can explain this variation, while the denser oceanic lithosphere subducts further to the south. However, it is still under question how the continental/oceanic lithospheric boundary manifests itself at depth. WHSZ has been the subject of great interest during the last decades. The deep structure of the WHSZ has been investigated by broad-scale surface and body wave tomographic studies (Spakman et al. 1988, 1993; Bijwaard et al. 1998; Pasyanos & Walter 2002; Piromallo & Morelli 2003; Schmid et al. 2006; Di Luccio & Pasyanos 2007; Endrun et al. 2008; Hosa 2008; Konstantinou & Melis 2008; Li et al. 2008; Koulakov et al. 2009; Salaün et al. 2012; Zhu et al. 2015). These studies show a common feature of systematic fast velocities extending down to ∼1400 km into the lower mantle beneath WHSZ, interpreted as the subducted slab. The most recent high-resolution tomographic constraints of the shallow part of the system (<250 km depth) were derived by Suckale et al. (2009) and Pearce et al. (2012), using an inversion technique based on the Generalized Radon Transform (GRT-RF) on teleseismic recordings obtained from the MEDUSA experiment, Halpaap et al. (2018), who applied a body-wave local tomography method, and Hansen et al. (2019), who used teleseismic and regional P-wave traveltime tomography over an adaptive grid. According to the above studies, WHSZ is composed of two lithospheric segments, continental to the north and oceanic to the south, with the separating boundary defined roughly at the northern tip of CTF (Halpaap et al. 2018). Several authors (e.g. Spakman et al. 1988; Piromallo & Morelli 2003; Govers & Wortel 2005; Suckale et al. 2009; Royden & Papanikolaou 2011) consider this boundary as a Subduction-Transform-Edge-Propagator (STEP) fault due to a trench-perpendicular tearing of the slab or a ramp-type structure (Pearce et al. 2012; Halpaap et al. 2018) beneath CTF, bridging NAT with the Hellenic Trench (Jolivet et al. 2013). Hansen et al. (2019) argue that a vertical tear occurs beneath CTF, but a trench-parallel one exists from northern Greece to north Peloponnese. Interpretations of this matter based on mantle anisotropy vary (Hatzfeld et al. 2001; Endrun et al. 2011; Olive et al. 2014; Paul et al. 2014; Evangelidis 2017; Kaviris et al. 2018). Sachpazi et al. (2016) suggest multiple vertical tearing (segmentation) of the slab occurring on along-dip, transverse faults. Moreover, vertical tearing along-dip of the slab has been proposed to exist beneath the Movri Fault Zone (MFZ), a NE–SW primary dextral strike-slip fault system in NW Peloponnese (Ganas et al. 2009), interpreted as the result of differential shear slip at the base of the two crusts (Durand et al. 2014; Sachpazi et al. 2016). 2.2 The East Hellenic Subduction Zone (EHSZ) The EHSZ extends from Crete to Rhodes. Seismicity, volcanism and tomography imply a different structure between the western and eastern part of the arc, which has been a matter of debate. However, the location of the transition between the two segments is unresolved. Papazachos & Nolet (1997), Meier et al. (2007), Brüstle (2012) and Sodoudi et al. (2015) explain the slab's increasing dipping angle from west to east by a first-order segmentation in SE Aegean. Tomographic studies have suggested the presence of a N–S vertical tear separating the HSZ from the Western Cyprus Subduction Zone (e.g. Piromallo & Morelli 2003; van Hinsbergen et al. 2010; Biryol et al. 2011; Legendre et al. 2012; Govers & Fichtner 2016). On the contrary, Faccenna et al. (2014) have suggested horizontal tearing of the slab in the area of Rhodes to explain the presence of a NW dipping slab. Endrun et al. (2008) in a surface wave study, suggest this boundary to occur at the longitude of central Crete. Bocchini et al. (2018), using well-located hypocentres from global and local seismicity catalogues, suggest a first-order slab segmentation east of Crete, with the eastern segment dipping with a steeper angle compared to the western one. A prominent feature of EHSZ is the South Aegean Active Volcanic Arc (SAAVA), a complex of Plio-Quaternary volcanoes aligned in an arcuate shape, whose activity is considered to be related to the subduction of the Tethys beneath the Aegean. SAAVA extends for 450 km, from Methana to Nisyros (Fig. 1), and is characterized by calc-alkaline volcanic activity (e.g. Fytikas et al. 1985; Mitropoulos et al. 1987; Mitropoulos & Tamey 1992; Francalanci et al. 2005). Variation in magmatism is observed between the western and eastern part of SAAVA. The western part of the arc (Aegina, Methana) is related to Pliocene age andesite-dacite volcanism, while geochemical studies have showed that magmas result from the upwards migration of fluids in the asthenospheric mantle wedge (Pe-Piper & Piper 2002; Tzanis et al. 2018). The central and eastern parts of the arc (Milos, Santorini, Nisyros) consist of lavas such as tholeiitic and calc-alkaline minor basalt, andesite, dacite and minor rhyolite, originated by both hydrated (calc-alkaline composition) and depleted asthenospheric (tholeiitic magmas) mantle. The northern margin of SAAVA correlates well with a reversal of rotation (Kissel & Laj 1988), the transition of directions of the maximum horizontal compressional stress axis at crustal depths (SHmax; Kapetanidis & Kassaras 2019; see Fig. S5a in the Supporting Information) and the respective strain-rate field (εH1, Chousianitis et al. 2015) from roughly E–W compression in the north to almost N–S in the south. Moreover, this transition is related to a sharp contrast of the crustal stress-shape (Kapetanidis & Kassaras 2019) and deflection of SKS (Evangelidis 2017; Kaviris et al. 2018) and Rayleigh wave (Endrun et al. 2011) anisotropy orientation. 2.3 The backarc region The backarc region includes N. Greece, central/north Aegean and N. Anatolia. The westward motion of Anatolia and extension of the Aegean/N. Greece dominates the geodynamics of this region, partly driven by asthenospheric flow and extrusion of a lid of rigid crust (Jolivet et al. 2013). Major normal faults, which accommodate crustal extension in the backarc, control the formation of ‘domino-style’ onshore and offshore Neogene sedimentary basins (e.g. Westaway 1991) in a series of E–W to NW–SE half-grabens. The most important one is the Corinth rift, which opens rapidly in a ∼N–S direction (e.g. Armijo et al. 1996) at rates reaching 15 mm yr−1 at its western part (Briole et al. 2000; Avallone et al. 2004). The North Aegean Trough (NAT) marks the westward propagation of the North Anatolian Fault (NAF), a chief characteristic of the tectonics of the SE Mediterranean (Taymaz et al. 1991; Papazachos & Kiratzi 1996; Papadimitriou & Sykes 2001; Goldsworthy et al. 2002). NAT and CTF systems are bridged by extensional grabens of the Central Hellenic Shear Zone (CHSZ), which facilitate the relay of plate motion (Shaw & Jackson 2010) and result in a general clockwise rotation of the area (Chousianitis et al. 2015). Anatolia is part of the Alpine–Himalayan orogenic belt, consisting of continental lithosphere formed during the Cenozoic (e.g. Robertson & Dixon 1984). It is located between Arabia and Eurasia, the collision of which leads to lateral extrusion that translates Anatolia westwards, towards the Aegean, along major strike-slip faults such as the North Anatolian Fault (NAF, Fig. 1, van Hinsbergen et al. 2009). The geological architecture of Anatolia is different from that of Greece, being composed of a complex distribution of continental and oceanic lithosphere (Jolivet et al. 2004) while underthrusting below Eurasia is associated with high-temperature metamorphism and granitic magmatism (van Hinsbergen et al. 2010). 3 DATA AND METHOD The surface wave analysis aims at estimating the seismic shear wave velocity by solving the inverse problem based on measured dispersion curves (Foti et al. 2018). The analysis is typically a three-step procedure: acquisition/processing of waveform recordings, dispersion curves estimation and inversion of Rayleigh wave dispersion curves for the determination of VSV velocity models with depth. We have used recordings from broadband seismological networks operating in the region of Greece, complemented by networks of neighbouring countries, that is Turkey, Bulgaria and Albania, of earthquakes that have occurred at epicentral distances between 5° and 160°, mostly teleseismic ones (Fig. 2a). The waveform data were acquired from the European Integrated waveform Data Archive (EIDA) node of GI-NOA (http://eida.gein.noa.gr/) via the International Federation of Digital Seismograph Networks web-services (FDSN-WS) data request protocol. The facilities of Bogazici University Kandilli Observatory and Earthquake Research Institute Data Services (KOERI DS) were also used for accessing waveform and metadata. KOERI waveforms were retrieved through the EIDA node located at KOERI, via the web services provided by the FDSN-WS. Additional recordings were acquired from stations of the Corinth Rift Laboratory Network (CRLN), operating in Central Greece, as well as of GEOFON and MEDNET, mainly to cover southern Greece and the northern borders of Greece. Figure 2. Open in new tabDownload slide (a) Locations of the earthquakes used for the dispersion analysis. The red rectangular indicates the study area. Inset rose-diagram shows the distribution of event backazimuths with respect to the centre of the study area. (b) The valid number of CC per station used. Figure 2. Open in new tabDownload slide (a) Locations of the earthquakes used for the dispersion analysis. The red rectangular indicates the study area. Inset rose-diagram shows the distribution of event backazimuths with respect to the centre of the study area. (b) The valid number of CC per station used. The data set extracted from the above databanks consists of 2-hr waveform windows for ∼1000 regional and teleseismic earthquakes with Mw ≥ 5.0 and Mw ≥ 6.0, respectively, that have occurred during 2010–2018 at epicentral distances between 5° and 160° from the centre of the network (Fig. 2a). Given that this study focuses on the seismic structure of the upper mantle of the Greek region, the analysis is performed in the range 30–90 s. Over 150 000 vertical component teleseismic seismograms recorded by ∼200 stations were pre-filtered between 0.01 and 0.33 Hz, decimated down to 1 sample s–1 and corrected for the instrumental response to unify recordings of different types of seismometers and changes in the instrumentation throughout the years. The second step of the surface wave analysis is the estimation of Rayleigh wave interstation phase velocities. Herein, we constrain dispersion curves by applying the Automated Surface Wave Measuring System algorithm (ASWMS; Jin & Gaherty 2015), based on the Generalized Seismological Data Functionals (GSDF) proposed by Gee & Jordan (1992). In this method, phase-delays are measured between all combinations of neighbouring stations of the network to avoid systematic biases due to propagation effects, rather than only for pairs that are aligned with the epicentre, within a few degrees of tolerance along a common great-circle path, an assumption that would greatly reduce the number of quasi-valid measurements. Besides, the simultaneous consideration of paired relative amplitudes allows greater precision for the Rayleigh wavefield, compared to applying only the phase information (e.g. Pollitz & Snoke 2010). The ASWMS method has the advantage that it requires minimal analyst interaction and can be applied to arrays with apertures ranging from continental scale, such as the EarthScope's TA, down to a few hundred metres (Jin & Gaherty 2015), yielding consistent results (http://ds.iris.edu/ds/products/aswms). The procedure requires the input of broad-band displacement recordings at a large multitude of stations and involves phase and amplitude difference measurements for phases of interest between each pair of seismograms (S1, S2) recorded on neighbouring stations. First, a window function WS is applied to crop the signal of interest (on S2), roughly based on the theoretical arrival-time of the primary surface wave and a length broad enough to include most of its coda, which can be highly correlated at interstation distances within one to two wavelengths (Jin & Gaherty 2015). WS is generated by frequency–time analysis to each waveform (Levshin et al. 1992) and group delays are isolated for the selected range of frequencies. Next, the cross-correlation function, C(t), between the two seismograms is calculated, containing the time-delays of coherent signals. Each cross-correlogram C(t) is then cropped with a bilateral ‘Hanning taper’ windowing function Wc centred on the C(t) maximum, which corresponds to a group delay whose central frequency depends on the dominant energy of the input data. The windowed C(t) is then convolved with a sequence of narrow-band Gaussian filters, for different central frequencies or periods of interest, providing a series of filtered cross-correlograms which contain information about the frequency-dependent group and phase-delays for the specific interstation path. Following GSDF (Gee & Jordan 1992), each filtered cross-correlogram is then modeled, using non-linear least-squares inversion, by a five-parameter wavelet. The latter is the product of a Gaussian envelope and a cosine function parametrized by the group and phase-delays between the two stations, a scale factor, the half-bandwidth and the central frequency of the Gaussian filter (Gee & Jordan 1992). Data selection is performed using a quality control criterion based on frequency-dependent coherence. For this purpose, the autocorrelation function (AC) of each signal is calculated and fitted by the five-parameter wavelet. Coherence, C, is estimated for each period by comparing the amplitude XC of the respective narrow-band filtered cross-correlation between the autoselected waveform windows for pairs of nearby stations and the respective amplitudes A1, A2 of the autocorrelations, as C = XC2/(A1 A2). Pairs with C < 0.5 were discarded, thus eliminating measurements with low SNR and from malfunctioning stations. The constructed 2-D network of frequency-dependent interstation phase-delays for each earthquake creates a phase gradient that can be inverted to provide a field of slowness that describes the surface wave propagation, with the vector length on each node of the grid being reciprocal to the apparent phase-velocity, according to the Eikonal equation (Jin & Gaherty 2015). This procedure is appropriate on the assumption that the used network of seismological stations can be considered as a seismic array. The latter is true for interstation distances smaller than three wavelengths for the smallest period analysed, that is ∼360 km for surface waves with a period of 30 s, taking into account an average propagation velocity of 4.0 km s–1. For this purpose, signals from station-pairs with a distance larger than 300 km are discarded and do not pass through the cross-correlation stage. To further avoid cycle-skipping issues, which can be a problem especially on shorter periods, phase-delays are compared to the expected traveltime corresponding to the interstation distance and, if needed, are corrected by ±1 or up to ±2 periods. Phase-velocities are derived from the inverted phase-delays in terms of apparent and structural velocities, approximated by the application of the Eikonal and Helmholtz equation, respectively. The Eikonal equation is valid for smoothly varying velocity structures; however, the resulting apparent phase-velocities can likely be contaminated by multipathing or backscattering effects, which invalidate the latter assumption and create biases (Lin et al. 2009). To avoid this issue, the structural phase-velocities are calculated by means of an approximation of the Helmholtz equation that takes into account amplitude correction terms (Jin & Gaherty 2015). For the calculation of the latter, the amplitude parameter of the wavelets which have been fitted on the narrow-band filtered AC is taken into account as the amplitude measurement at each station. Moreover, site conditions and problematic instrument responses, which can at times introduce non-negligible biases, were taken into account by determining the ratio between single-station amplitude measurements and the median amplitude values of all other stations for each event/period. These ratio values were found to be consistent at each station, except for a few outliers on some events. No significant changes were found throughout time that could be attributed to a serious issue on a station's sensor or to undocumented changes in the instrumentation. The bias of each station was derived by the median of the log10 of the individual ratios for each period, found to be in the range ±0.17, or, in true values, varying from 0.68 to 1.46, further reduced to a range from 0.9 to 1.1 after the application of minimum curvature surface fit smoothing on the grid nodes (for example see Fig. S10 in the Supporting Information). After a first run of the procedure, the calculated offsets were used to remove the systematic biases from the measured amplitudes on individual events (see e.g. Accardo et al. 2017). Amplitude values which lie below ½ or are over two times the median of amplitudes of neighbouring stations at valid cross-correlation pairs, are discarded. The retained amplitudes are then ‘interpolated’ at each grid point by fitting a minimum curvature surface using a fourth order partial derivative equation (PDE, Lin & Ritzwoller 2011; Jin & Gaherty 2015), that reduces sharp effects often caused by second-order nonlinear PDEs (Lŭ et al. 2008). The correction term is calculated and then subjected to strong smoothing (see e.g. Fig. A2 and Figs S13, S15 in the Supporting Information), to reduce short-wavelength variances relative to the resolution of the network (Jin & Gaherty 2015) before it is finally applied, yielding the structural phase-velocity maps for each event/period. The last stage of the dispersion analysis involves stacking, to produce the final structural phase-velocity maps for each period. The phase-velocity maps of individual events are averaged through weighted stacking, with the weight of each node determined from the ray density during the slowness inversion (Jin & Gaherty 2015). The presented stacked maps of either apparent or structural phase-velocity finally pass through a last smoothing filter, with a smoothing length of ¼ the average wavelength. 3.1 Phase-delay measurements and resolution Following the above-mentioned methodology, a total of 4320 765 CC measurements were found to be valid across the network over the period range of interest (30–90 s, Fig. 3a), with an average of ∼617 252 valid CC per examined period, selected with a frequency-dependent coherence above 0.5. The adopted cross-correlograms were passed through a sequence of various narrow-band, zero-phase Gaussian filters, with their shape constrained between 6 per cent and 10 per cent of their central frequency, the latter being defined by each period of interest (30–90 s). Then, phase-delays from all the station-pairs were assembled and their linear trend against interstation distance was calculated for each period. Measurements having a phase-delay misfit higher than 10 s, relative to the least-squares phase-delay prediction, were rejected (Fig. 3b). This is a conservative threshold value, as it rejects phase-delay measurements which are offset by over 20 per cent for station-pairs which are less than 200 km apart (Jin & Gaherty 2015). Having started with a database of ∼1000 regional/teleseismic events with M ≥ 6, almost 20 per cent of the initial events data set (∼200 events) was rejected after the application of the selection criteria (see Section 3). Figure 3. Open in new tabDownload slide (a) Number of valid CC measurements across the region, (b) example of data selection in the period range of interest by rejecting cross-correlation measurements having a phase delay misfit higher than 10 s (grey circles). Figure 3. Open in new tabDownload slide (a) Number of valid CC measurements across the region, (b) example of data selection in the period range of interest by rejecting cross-correlation measurements having a phase delay misfit higher than 10 s (grey circles). Periods <30 s were initially taken into account, but a large proportion of them were discarded by the applied criteria and it was decided to proceed only for T ≥ 30 s. The main issue for discrepancies in the lower periods is the large spacing between the stations due to morphology limitations, that is the Aegean Sea, often larger than the wavelength of the high-frequency surface wave of interest, but also strong attenuation and scattering at short periods and teleseismic distances. Large period maps (>90 s) were also discarded due to the scarcity of paths sampled by long-period instruments, which is why many of the resulting cross-correlations were considered as unreliable. After a first series of cross-correlation measurements was performed, the results were examined to detect stations which could be problematic, that is producing a large amount of invalid measurements compared to the acceptable ones. This is important, as the ratio between ‘good’ and ‘bad’ measurements is a quality control criterion that can reject a whole phase-velocity map for a certain event/period, while problematic data can negatively affect the final maps, biasing average values and increasing the standard deviations and relative errors. A total of 34 stations, whose pairs produced a valid/invalid ratio of less than 0.20, were removed and the whole procedure was repeated without them. The seismological network, herein considered as an array, has an aperture broad enough (1000 × 600 km2, see Fig. 2b) to resolve lateral phase-velocity variations up to ∼100 s period, with a lateral resolution which is mainly governed by the interstation distance (e.g. Pedersen et al. 2003; Bodin & Maupin 2008). It is apparent from Fig. 4 that the region that is better sampled by interstation links is particularly well constrained around the Corinth gulf and W. Anatolia, while N. Greece and S. Peloponnese are partially covered. On the other hand, the regions of Southern Aegean and Dodecanese Islands are poorly sampled, since the network is very sparse there due to the sea cover. It is noticeable that valid pathways appear increased at higher periods, given that dispersion measurements are more stable over the long wavefield (e.g. Shapiro & Ritzwoller 2002). A criterion of the contribution of teleseismic events per grid cell was applied as a first order approximation of the ray density towards the herein resolution analysis. A threshold value of 10 contributing events per grid cell has been considered towards retaining/rejecting measured regions (Fig. 5). Figure 4. Open in new tabDownload slide Valid CCs per interstation pathway for selected periods of interest in the study area. Figure 4. Open in new tabDownload slide Valid CCs per interstation pathway for selected periods of interest in the study area. Figure 5. Open in new tabDownload slide Seismic events used per grid cell for various period (T) intervals. Figure 5. Open in new tabDownload slide Seismic events used per grid cell for various period (T) intervals. An important parameter for the inversion of phase-delays to apparent phase-velocities is the damping parameter. Weak damping can result in overfitting, producing small data misfit but an unrealistic model, while strong damping produces a smoother model-data convergence but also larger errors between the two. It is a common practice during inversions to select a damping parameter which produces a compromise between a smaller ‘model length’ (a less complex model) and small data misfit, that is the ‘knee’ of the, optimally, L-shaped trade-off curve. To define the appropriate smoothing scheme, the inversion was repeated for each period, by varying the damping factor over a wide range of values. Data misfit was calculated by the differences between predicted and measured phase-delays, with the errors being weighted by the amplitude of the respective cross-correlation maxima. Model length was measured as the variance of the final apparent phase-velocity maps for each event and period. Data misfit and model length were then averaged over all events for each period and damping factor. Diagrams of the trade-off between model length and data misfit were then used to determine the optimal damping parameter (Fig. S11). These diagrams normally have the shape of an ‘L’ and the optimal parameter is the one that corresponds to the ‘knee of the curve’. Disregarding discrepancies observed for very small damping parameters, the curves for periods between 30 and 55 s have the appearance of a smooth L-shaped curve. Larger periods exhibit increasing or almost constant model variance for increasing damping parameter, while the data misfit also becomes larger. In this case, values that produced a smaller model length after its initial increase were selected. 3.2 Checkerboard tests Checkerboard tests (Spakman & Nolet 1988) are the most commonly used model for synthetic tests to investigate the robustness of a tomographic solution along each spatial dimension of an a priori velocity model, that is geometry, polarity and amplitude of the reconstructed perturbations. To reconstruct the synthetic model, a procedure proposed by Accardo et al. (2017) is applied. The method examines, for each event, the ray coverage of the real data, in terms of -valid- cross-correlations between pairs of stations for every period. For each of these valid combinations, the real (observed) relative travel-ime differences are taken into account from which the real slowness field is calculated for each cell. The slowness vectors (reciprocals of the phase-delay gradient) are normalized by dividing with the real slowness and are then multiplied by the respective slowness value of the a priori checkerboard model and converted back into values of the phase-delay gradient. The time-lag weights that are used are the same as those measured with the real data. The rest of the procedure is then similar to the one followed for the determination of the apparent phase-velocity from the phase-delay gradient through inversion, with a map produced for each event and period that is then stacked to produce the final reconstruction. During stacking, the results are smoothened by a factor that is proportional to ¼ of the wavelength, as is also performed for the respective apparent phase-velocity maps with the real data. Although this is not a typical model reconstruction (no real synthetic travel-times are calculated or noise applied) it provides a qualitative estimate of how well the available rays of valid cross-correlations between stations, their weights and the applied smoothing, is able to reconstruct a certain pattern of velocity anomalies. We performed several checkerboard tests in order to obtain the optimum cell size of the resolved models, in terms of the degree of the achieved consistency of the inverted synthetic anomalies with the initial checkerboard model for both signed amplitude and shape. An average velocity of 3.9 km s–1 was considered for each examined period, with the synthetic anomalies alternating between 3.5 and 4.3 km s–1, representing a ±10 per cent perturbation relative to the average value of the observed phase-velocity and its variation. Fig. 6 presents the synthetic tests, showing the capability of the method to resolve structures given the available real data distribution, which is over phase-velocity anomalies with dimensions 1.35° × 1.35° for periods 30–50 s and 1.80° × 1.80° for periods 60–90 s. It can be observed that for such a checkerboard test configuration, the shape and polarity of the periodic intervals are well reconstructed in some parts of the study area, including continental Greece and W. Anatolia. In the Aegean area, smearing occurs due to the sea coverage, while cells in N. Greece and S. Peloponnese are not optimally recovered. The limited resolution in these regions is taken into account in the interpretation of the results. Figure 6. Open in new tabDownload slide Results of a checkerboard resolution test using the source-receiver geometry of the valid CC pairs, (top panel) with phase-velocity anomaly size 1.35°×1.35° for periods 30–50 s, (bottom panel) with phase-velocity anomaly size 1.80°×1.80° for periods 60–90 s. The a priori phase-velocity pattern is shown at the panels on the left. Unconstrained regions of the model have been omitted. Figure 6. Open in new tabDownload slide Results of a checkerboard resolution test using the source-receiver geometry of the valid CC pairs, (top panel) with phase-velocity anomaly size 1.35°×1.35° for periods 30–50 s, (bottom panel) with phase-velocity anomaly size 1.80°×1.80° for periods 60–90 s. The a priori phase-velocity pattern is shown at the panels on the left. Unconstrained regions of the model have been omitted. 4 INVERSION OF PHASE-VELOCITIES FOR SHEAR WAVE VELOCITY (VSV) The third step of the analysis concerns the depth inversion of the phase-velocity maps in order to obtain the vertically polarized shear velocity (VSV) with depth, following a scheme implemented in the ASWSM suite (Jin et al. 2015). The surface wave depth inversion problem is highly non-linear and is affected by solution non-uniqueness. Although the error of phase-velocities propagates into the final model, the main sources of ambiguity in the final model are the non-uniqueness of the solution due to the starting model, the smooth nature of the 1-D sensitivity kernels and data errors. Different starting models may yield different final models but with similar goodness of fit with the experimental data (Maupin 2016). Herein, the problem of non-uniqueness of the inversion due to the starting model is treated by a global search method applied to evaluate large ensembles of 1-D earth models within defined parameter ranges, searching for those that produce reasonable misfits with the experimental data. To build starting 1-D models, for each grid cell, a grid-search and forward modeling is applied by trying out combinations of different velocities and thicknesses on a three-layered 1-D VSV model, representing a layer of sediments, crust and uppermost mantle, selecting as the best model the one with the smaller misfit in the derived dispersion curve with respect to the one obtained by the structural phase-velocity maps. To allow smooth velocity variations with depth, we subdivide the crustal and mantle layers of each initial model into sub-layers of even thickness, that is layers with a thickness of 10 km—based, however, on the results of the crude three-layer model—and a more detailed ‘final model’ is constructed by applying inversion using the surf96 algorithm. Information on the Moho topography is taken into account, disabling smoothing on the layer that corresponds to the Moho discontinuity, to obtain a sharp velocity contrast in the ‘final model’. For this purpose, the CRUST1.0 model (Laske et al. 2013; Fig. S3 in the Supporting Information) is used. CRUST1.0 (https://igppweb.ucsd.edu/∼gabi/crust1.html) is a 1° × 1° global crustal model, based on constraints from active seismic sources, receiver functions, and gravity. It is considered as a realistic compromise of first-order crustal thickness variations in the area of study, since it is generally in good agreement with results from pertinent research for the region, implying that the Moho depth in the area ranges from 20 to 45 km (Papazachos & Nolet 1997; Clement et al. 2004; Karagianni et al. 2005; Zelt et al. 2005; Sodoudi et al. 2006; Sachpazi et al. 2007). Disabling smoothing on the Moho discontinuity produces an abrupt change in the VSV model, with the layer below Moho being more influenced by velocities of deeper layers rather than of crustal layers. On the other hand, shallower layers above Moho mainly refer to depths <40 km, which are not well resolved due to lack of short period data. This creates a trade-off between crustal and upper-mantle structure which may cause differences between phase-velocity maps and VSV at shallow depths to be observed. A simple inversion using a single initial model provides some first-order information on the problem, such as the resolving kernel, which shows the sensitivity of the obtained results for each layer to a range of real depths. Ideally, the resolution matrix would be an identity matrix and each layer would be sensitive to its own depth range and not affected by the structure at different depths; in such a case, each row of the resolution matrix would contain a single spike on the resolution matrix diagonal element (RDE). In practice, however, the resolution kernel is represented by a smooth curve, with a maximum ideally near the RDE and non-zero off-diagonal elements (ROEs, e.g. Fig. S17). The width of the curve or the offset of its maximum from the diagonal indicate a degree of influence/biases of the result on a certain layer with respect to the structure in other layers (more information on the resolution matrices in the Supporting Information; Figs S4, S17–S21). However, an inversion based on a single starting model can be greatly biased by the choice of that model. To remove the factor of starting model selection, a Monte Carlo technique is applied. The main parameters (shear velocity and thickness) of the multilayered (as in the previous step) best 1-D initial model are perturbed within a limited range to generate a set of 100 random models. Forward modeling is applied on each of these 100 initial models to examine whether the resulting rms value is smaller than a threshold value set at 0.3 km s–1, otherwise the model is discarded. For starting models with rms < 0.3 km s–1, an iterative damped least-squares inversion is applied on the observed dispersion curves, taking also into account the respective standard deviations obtained during stacking of the Helmholtz maps as the error of the observed velocity measurements, to derive the final shear wave model with depth, together with its statistical error distribution, by applying the surf96 code (Herrmann 2013). After 10 iterations, the inversion produces an 1-D VSV ‘final’ model for each of the 100 starting models per grid-cell. Finally, another error filtering is applied on each node, by rejecting ‘final’ models which have an rms between observed and calculated dispersion curves larger than 20 per cent with respect to the median rms misfits of ‘final’ models with rms < 0.3 km s–1. The retained models are then interpolated into a common, more detailed depth grid and averaged into a final VSV model and its errors in terms of the standard deviation for each layer. Similarly to the phase-velocity maps, the horizontal resolution of the 3-D model is considered to be of the same order as that of the recovered checkerboards (Fig. 6). However, the vertical resolution of the shear velocity from the surface wave analysis is highly variable with depth, being reversely proportional to the wavelength, and, consequently, depth. The above implies that surface wave provides information about bulk properties in a given depth range rather than sharp images of the subsurface (Maupin 2016). We addressed this issue by quantifying the precisely resolved depth range and the resolution obtained at each depth, by examining the resolution matrix per grid node, constructed using the surf96 code and the velocity resolving kernel output. The latter analysis showed that in our model best resolution is achieved at depths between 40 and 100 km, while the structure above 40 km and below 140 km is not represented by our data set (Fig. S4 in the Supporitng Informtion). 5 AVERAGE STRUCTURE The average apparent and structural phase-velocities per period interval over the whole study area (Table 1) appear consistent with Bourova et al. (2005), and Salaün et al. (2012). Fig. 7(a) presents the average apparent and structural phase-velocity curves for the Greek region, showing that the Helmholtz tomography has yielded slightly lower values for periods >50 s, by 0.01–0.03 km s–1. The latter is negligible, taking into account that this difference is much smaller than one standard deviation. The herein Helmholtz velocities (Fig. 7b) appear to be lower, when compared to several global and regional continental dispersion models, consistently with tomographic results that show slow upper-mantle velocities for the most active features in Europe (e.g. Hearn 1999; Piromallo & Morelli 2003; Koulakov et al. 2009; Zhu et al. 2015). Figure. 7. Open in new tabDownload slide (a) The average apparent (Eikonal) and structural (Helmholtz) phase-velocity curves derived in this study for the study region; (b) Comparison of the average structural (Helmholtz) phase-velocity model with global and regional dispersion curves. Figure. 7. Open in new tabDownload slide (a) The average apparent (Eikonal) and structural (Helmholtz) phase-velocity curves derived in this study for the study region; (b) Comparison of the average structural (Helmholtz) phase-velocity model with global and regional dispersion curves. Table 1. The average Rayleigh-wave phase-velocity dispersion model of the study region in 10 s intervals. . . Apparent velocity (Eikonal) . Structural velocity (Helmholtz) . Period (s) . No of valid CCs . Average PhV (km s–1) . St dev. (km s–1) . Average PhV (km s–1) . St dev. (km s–1) . 30 527,726 3.71 0.09 3.71 0.09 40 630,728 3.83 0.06 3.83 0.06 50 653,278 3.91 0.08 3.91 0.08 60 645,078 3.97 0.09 3.96 0.10 70 634,302 4.01 0.10 4.01 0.10 80 621,901 4.05 0.10 4.04 0.10 90 607,752 4.08 0.10 4.07 0.09 . . Apparent velocity (Eikonal) . Structural velocity (Helmholtz) . Period (s) . No of valid CCs . Average PhV (km s–1) . St dev. (km s–1) . Average PhV (km s–1) . St dev. (km s–1) . 30 527,726 3.71 0.09 3.71 0.09 40 630,728 3.83 0.06 3.83 0.06 50 653,278 3.91 0.08 3.91 0.08 60 645,078 3.97 0.09 3.96 0.10 70 634,302 4.01 0.10 4.01 0.10 80 621,901 4.05 0.10 4.04 0.10 90 607,752 4.08 0.10 4.07 0.09 Open in new tab Table 1. The average Rayleigh-wave phase-velocity dispersion model of the study region in 10 s intervals. . . Apparent velocity (Eikonal) . Structural velocity (Helmholtz) . Period (s) . No of valid CCs . Average PhV (km s–1) . St dev. (km s–1) . Average PhV (km s–1) . St dev. (km s–1) . 30 527,726 3.71 0.09 3.71 0.09 40 630,728 3.83 0.06 3.83 0.06 50 653,278 3.91 0.08 3.91 0.08 60 645,078 3.97 0.09 3.96 0.10 70 634,302 4.01 0.10 4.01 0.10 80 621,901 4.05 0.10 4.04 0.10 90 607,752 4.08 0.10 4.07 0.09 . . Apparent velocity (Eikonal) . Structural velocity (Helmholtz) . Period (s) . No of valid CCs . Average PhV (km s–1) . St dev. (km s–1) . Average PhV (km s–1) . St dev. (km s–1) . 30 527,726 3.71 0.09 3.71 0.09 40 630,728 3.83 0.06 3.83 0.06 50 653,278 3.91 0.08 3.91 0.08 60 645,078 3.97 0.09 3.96 0.10 70 634,302 4.01 0.10 4.01 0.10 80 621,901 4.05 0.10 4.04 0.10 90 607,752 4.08 0.10 4.07 0.09 Open in new tab In Fig. 8 we compare the inverted average VSV structure of the study region with two 1-D shear velocity global reference models, namely the AK135 (Kennett et al. 1995) and SPREM-VSV (Ho et al. 2016) models. The average 1-D model exhibits slow velocities for the whole depth range concerning AK135, and at depths ∼40–80 km compared to SPREM-VSV. Low S-wave uppermost mantle values were also found by Karagianni & Papazachos (2007) and Salaün et al. (2012) for the Aegean region and by Mindevalli & Mitchell (1989) for Turkey. Additional average velocity models and dispersion curves for selected subregions are included in the Supporting Information (Figs S5b and S6, Tables S1 and S2). Figure 8. Open in new tabDownload slide Average 1-D model of weighted VSV determined using [1–st.dev.(VSV)/mean(VSV)] as weight, in juxtaposition with AK135 and SPREM reference models. Grey bold dashed lines indicate the depth range where the VSV model is resolved (40–140 km). Figure 8. Open in new tabDownload slide Average 1-D model of weighted VSV determined using [1–st.dev.(VSV)/mean(VSV)] as weight, in juxtaposition with AK135 and SPREM reference models. Grey bold dashed lines indicate the depth range where the VSV model is resolved (40–140 km). 6 STACKED APPARENT AND STRUCTURAL PHASE VELOCITY MAPS We examine apparent and structural phase-velocity variations for periods of 30–90 s, that is for wavelengths between approximately 120 and 450 km. The final tomograms are the result of weighted stacking of the respective single-event Eikonal and Helmholtz maps. Some examples of single-event phase-velocity maps are presented in Appendix A and in the Supporting Information, to demonstrate different steps of the procedure as well as some characteristics of the Rayleigh wavefield as it propagates through the Aegean region. Figs 9 and 10 present maps of apparent and structural phase-velocities, respectively, at various periods of interest, derived from the Helmholtz tomography. These maps were constructed using a colour scale that changes depending on the average phase-velocity at each period. As a consequence of the approach, spatial variations of phase-velocities appear similar to velocity perturbations, while absolute phase-velocity values are displayed simultaneously. Figure 9. Open in new tabDownload slide Rayleigh-wave apparent (Eikonal) phase-velocity maps at different periods. The maps are plotted for event number ≥10. The colour scale varies with the period (T). White colour corresponds to the average phase-velocity shown in the lower right-hand corner of each map. Figure 9. Open in new tabDownload slide Rayleigh-wave apparent (Eikonal) phase-velocity maps at different periods. The maps are plotted for event number ≥10. The colour scale varies with the period (T). White colour corresponds to the average phase-velocity shown in the lower right-hand corner of each map. Figure 10. Open in new tabDownload slide Rayleigh-wave structural (Helmholtz) phase-velocity maps at different periods. The maps are plotted for event number ≥10. The colour scale varies with the period (T). White colour corresponds to the average phase-velocity shown in the lower right corner of each map. Green contours indicate relative errors (per cent) with respect to the average phase-velocity per period of interest. Tick lines of the contours indicate the direction of error reduction. Figure 10. Open in new tabDownload slide Rayleigh-wave structural (Helmholtz) phase-velocity maps at different periods. The maps are plotted for event number ≥10. The colour scale varies with the period (T). White colour corresponds to the average phase-velocity shown in the lower right corner of each map. Green contours indicate relative errors (per cent) with respect to the average phase-velocity per period of interest. Tick lines of the contours indicate the direction of error reduction. Although the average values of phase velocities are similar, the Eikonal maps of Fig. 9 exhibit an amount of smearing for all periods of interest throughout the region. The source of the differentiation between Eikonal and Helmholtz maps is the application of the amplitude correction term on the apparent phase velocities, which produces finer spatial variations on the structural phase-velocities (see Figs S8, S13, S15 in the Supporting Information). Contours of relative error, in terms of the ratio between the standard deviation on each grid cell and the average phase-velocity value for each period of interest, are shown in Fig. 10 to indicate regions where the observed anomalies are significant, that is the relative error is smaller than the value of the anomaly. As shown in Fig. 10, the relative error appears systematically large in S. Aegean, likely due to the poor network coverage; however, a similar pattern has been previously explained by multipathing effects due to strong heterogeneities beneath this area (Bourova et al. 2005). Several features of morphologic, tectonic and geologic significance are indicated in the obtained phase-velocity maps. The geometry of these features is not well-recovered at the edge of our model and in regions with poor station coverage at sea. Thus, for these regions, interpretation is roughly based on their apparent continuity with well resolved neighbouring anomalies. The interpretation is thus limited to phase-velocity patterns that present smaller relative error than the amplitude of the anomaly. Both Eikonal and Helmholtz approaches depicts sufficiently the lithospheric plates that interact in the region, emphasizing the adequacy of surface wave analysis to reveal the lithospheric structure in complex regimes, such as the Hellenic region. Structures such as the low-velocity wedge above the slab that is prominent along the forearc at low periods (30 s), the subducted oceanic plate extending from N. Peloponnese to the coast of Turkey, the continental subduction in N. Greece, the overriding Aegean and the Anatolian Plate are visible in the phase-velocity maps. No high velocities are observed beneath southwestern Anatolia but only south of it, which implies that the vertical slab tear, that has been proposed by several authors (e.g. Piromallo & Morelli 2003; van Hinsbergen et al. 2010; Biryol et al. 2011; Legendre et al. 2012; Govers & Fichtner 2016), is likely limited to regions beneath continental Anatolia. Concerning the Helmholtz correction term, the regions with the largest anomalies, in absolute terms, are CLAB and S. Aegean (Fig. S8). Reversal of the anomaly polarity (towards higher structural velocity) is observed in S. Aegean, an area poorly sampled by the network, which indicates that these corrections may not be reliable. For smaller periods (T < 60 s) the anomalies in continental Greece north of Peloponnese propagation are uniform and discrete from the ones to the south. 7 RESULTS OF THE 3-D VSV MODEL—DISCUSSION The new 3-D VSV model yields interesting large-scale features of the Hellenic lithosphere, which are interrelated with the geodynamics of the region, and have also been the subject of scientific discussions/debates throughout the years. In the following paragraphs we focus on the main structures which are sufficiently resolved by the depth inversion of the Helmholtz phase-velocities: (a) the oceanic–continental lithosphere transition, (b) slab segmentation and (c) lithospheric thickness and the arrangement of the asthenosphere. These issues are discussed in juxtaposition with information on the distribution of subcrustal seismicity (catalogues compiled in this study), volcanism (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015) and the crustal stress-field inferred from the inversion of focal mechanisms (Kapetanidis & Kassaras 2019). The results of our 3-D VSV model are presented as horizontal depth slices between 40 and 140 km depth (Fig. 11). A series of profiles at various directions (Fig. 12) are used for the construction of vertical cross-sections (Fig. 13). Horizontal slices of absolute and relative VSV (Fig. 11) are drawn with a common colour scale at several depth intervals between 40 and 140 km depth. The colour scale of the maps is centred at the average velocity of each layer. Percentage of misfit for each specific depth is overlain on the maps as contours of relative error (standard deviation of the Monte Carlo tests with respect to the average VSV of each depth). Misfit values were found in the range 0.5–3 per cent, reaching 5 per cent in N. Greece only for the depth of 40 km, while the amplitude of the velocity anomalies ranges ± 9 per cent throughout the 3-D model. It should be noted that regions/depths with significant misfit in the phase-velocity model can have small misfit in the VSV model, as the relative errors of the latter only reflect the instability of the Monte Carlo solutions, which can be very low at large depths, despite the lower resolution. It is worth noting the difference between the shear velocity map at 40 km depth (Fig. 11) with phase-velocity maps (Figs 9 and 10), apparently the result of trade-off between the CRUST1.0 model (see Fig. S3 in the Supporting Information) and upper-mantle structure derived from herein Rayleigh wave dispersion data; this issue is expected to be resolved when lower period phase-velocities (10–25 s) become available in a forthcoming ambient noise tomography study. Figure 11. Open in new tabDownload slide Maps of absolute vertically polarized shear wave velocity (VSV) at different depths. The colour scale varies with depth Z. White colour corresponds to the average VSV as given in the lower right corner of each map. The colour bars tick lines are placed at ±3 per cent, ±6 per cent and ±9 per cent intervals relative to the corresponding average VSV. Green contours indicate the relative error (per cent) with respect to the average VSV per depth. Tick lines on the contours indicate the direction of error reduction. VSV velocities are smoothened horizontally within circles of radius 100 km on each grid cell. Figure 11. Open in new tabDownload slide Maps of absolute vertically polarized shear wave velocity (VSV) at different depths. The colour scale varies with depth Z. White colour corresponds to the average VSV as given in the lower right corner of each map. The colour bars tick lines are placed at ±3 per cent, ±6 per cent and ±9 per cent intervals relative to the corresponding average VSV. Green contours indicate the relative error (per cent) with respect to the average VSV per depth. Tick lines on the contours indicate the direction of error reduction. VSV velocities are smoothened horizontally within circles of radius 100 km on each grid cell. Figure 12. Open in new tabDownload slide Map of intermediate depth seismicity (Focal Depth ≥ 40 km) included in the herein used catalogues. Coloured triangles are volcanoes and thermal springs (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015). Black solid lines indicate profiles of cross-sections shown in Fig. 13. For more placenames discussed in the text refer to Fig. 1. Figure 12. Open in new tabDownload slide Map of intermediate depth seismicity (Focal Depth ≥ 40 km) included in the herein used catalogues. Coloured triangles are volcanoes and thermal springs (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015). Black solid lines indicate profiles of cross-sections shown in Fig. 13. For more placenames discussed in the text refer to Fig. 1. Figure 13. Open in new tabDownload slide Open in new tabDownload slide Cross sections on the 3-D VSV model along profiles shown on the map of Fig. 12, presented as velocity perturbations relative to the mean VSV value of the respective depth. Green contour lines and percentage labels represent the relative error of phase-velocities, translated at depth. Seismicity (this study, Karakonstantis 2017; Bocchini et al. 2018) as well as volcanoes and thermal springs (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015) are also plotted. Black dashed lines denote the Moho interface of the subducted lithosphere (after Halpaap et al. 2018). The white dashed line corresponds to the CRUST1.0 model of the upper plate Moho surface. VSV velocities are smoothened horizontally within circles of radius 100 km on each grid cell. See also Fig. S7 in the Supporting Information for cross-sections of absolute VSV values. Figure 13. Open in new tabDownload slide Open in new tabDownload slide Cross sections on the 3-D VSV model along profiles shown on the map of Fig. 12, presented as velocity perturbations relative to the mean VSV value of the respective depth. Green contour lines and percentage labels represent the relative error of phase-velocities, translated at depth. Seismicity (this study, Karakonstantis 2017; Bocchini et al. 2018) as well as volcanoes and thermal springs (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015) are also plotted. Black dashed lines denote the Moho interface of the subducted lithosphere (after Halpaap et al. 2018). The white dashed line corresponds to the CRUST1.0 model of the upper plate Moho surface. VSV velocities are smoothened horizontally within circles of radius 100 km on each grid cell. See also Fig. S7 in the Supporting Information for cross-sections of absolute VSV values. The global pattern of shear velocity anomalies in Fig. 11 is similar to that obtained by surface wave analysis, that is Bourova et al. (2005) and Kassaras et al. (2009), regarding the Aegean region and Salaün et al. (2012) over the Greek region and W. Anatolia. Compared to the latter, our model exhibits better resolution due to the applied methodology, which allows much larger input of dispersion measurements with respect to the typical two-station method. Resolved velocity anomalies are consistent with results from body-wave tomography (e.g. Papazachos & Nolet 1997; Halpaap et al. 2018; Hansen et al. 2019), regarding the imaging of lithospheric margins and slabs, although not in detail, which suggests that our data set and surface wave analysis are robust enough to depict the large-scale structure of the region's upper mantle for the intermediate depth. Our 3-D model clearly images the boundary between the subducting African slab and the Aegean lithospheric mantle, having an arcuate shape that crosses from N. Greece to Rhodes, separating two distinct structures, which are strongly linked with deformation of the crust. The latter is supported by the differences in the directions of crustal stress and strain-rate fields (i.e. the directions of the maximum horizontal stress/strain and the stress-shape values; Fig. S5a) between the northern and southern sides of this boundary (Kassaras et al. 2016; Kapetanidis & Kassaras 2019). To compare the velocity patterns obtained by teleseismic surface wave tomography with mantle dynamics and crustal deformation, we have collected four catalogues of reliable seismicity data (Figs 12 and S9): Hypocentral locations obtained for the broader study area between 2012 and 2017, using the manually picked P- and S-wave phases at the Seismological Laboratory of the National and Kapodistrian University of Athens (NKUA-SL), an optimum 1-D velocity model (Karakonstantis 2017) and HYPOINVERSE software (Klein 1989). Average hypocentral uncertainties of this catalogue are 0.9 and 3.0 km for the horizontal and vertical directions, respectively, and mean rms = 0.25 s; (see Fig. S2 in the Supporting Information for more statistics) the EGELADOS catalogue (Friederich & Meier 2008; Brüstle 2012); the LIBNET catalogue (Becker et al. 2010); the CYCNET catalogue (Bohnhoff et al. 2004, 2006; Brüstle 2012). The latter three catalogues include solutions having hypocentral uncertainties ≤15 km and rms ≤ 0.7 s (Bocchini et al. 2018). Fig. 13 illustrates cross-sections of shear velocity perturbations with respect to the average 1-D model along various profiles shown in Fig. 12. Cross-sections of absolute shear velocity along the same profiles are drawn in Fig. S7 of the Supporting Information. Only velocity anomalies with similar or larger dimensions than the ones retrieved by the synthetic checkerboard tests (Section 3.2) are taken into account. To facilitate the discussion and remove small-scale anomalies, additional smoothing is applied in the 3-D VSV model, averaging velocities on each horizontal layer within circles of radius 100 km centred on each grid point. Moreover, unrealistic anomalies that correspond to an error of phase-velocities larger or equal to the amplitude of perturbation are not taken into consideration. Anomalies with a perturbation percentage exceeding the relative error but which still represent plausible geological structures are retained and interpreted with caution. 7.1 The continental–oceanic subduction boundary The new 3-D shear velocity model provides constraints on the continental–oceanic boundary (COB) of the Hellenic lithosphere, which are discussed below. In the northwestern extremity of HSZ, slow/fast velocities are in agreement with the dip of Moho interface of the continental slab (Halpaap et al. 2018) down to the maximum depth resolved by the 3-D model (140 km, cross-section b1–b2 in Fig. 13). Sparse subcrustal seismicity distributed between 60 and 90 km depth also follows this pattern. Despite the reduced vertical resolution of the structure, inherent in surface wave analysis, this finding is consistent with Pearce et al. (2012) and Halpaap et al. (2018), who imaged the subducted continental crust in NW Greece as a low velocity zone (LVZ), dipping at a similar angle with the inferred velocities configuration. The pattern of consistency of the Moho interface (Halpaap et al. 2018) and intermediate seismicity with shear velocity contrasts continues to the south. Beneath N. Peloponnese, the arrangement of the above indicate a shallow-dipping slab (15°–20°) penetrating upper-mantle deeper than 140 km depth of our model (cross-sections d1–d2, f1–f2, in Fig. 13). In this region, WBZ is relatively shallow, reaching 90 km depth. This observation is consistent with Suckale et al. (2009), who interpret a shallow-dipping LVZ (∼20°) in this area as the hydrated crust of the Hellenic slab. Halpaap et al. (2018), attribute a lack of seismicity beyond 100 km depth for the northern HSZ to the effect of complete eclogitization of the crustal material, which limits dehydration embrittlement and consequently the potential occurrence of intermediate depth earthquakes. Eclogitization could reasonably be presumed, since it has been broadly adopted as a process to explain subduction of continental crust at convergence plate boundaries, when the oceanic crust has been consumed and slab-pull forces continue to exist. The overall arrangement of mantle velocities and the WBZ (Fig. 13, b1–b2, d1–d2, f1–f2) in the northern part of HSZ could be explained by a low-dipping continental lithosphere, whose eastern margin marks the western termination of NAT. In the area of S. Peloponnese, a ∼30° NE-dipping zone (Fig. 13, g1–g2, h1–h2), associated with the oceanic subduction zone, is defined by the Moho interface (Halpaap et al. 2018), the WBZ, that reaches 180 km depth west of Crete, as well as the shear velocity contrasts. The inferred configuration is very similar to that found in WHSZ by Pearce et al. (2012), Sachpazi et al. (2016) and Halpaap et al. (2018). The depth distribution of hypocentral locations provide implications on the characteristics of the transition between continental and oceanic subduction. Earthquakes deeper than 100 km in Greece are known to occur roughly southerly of the 39°N latitude (e.g. Bocchini et al. 2018; Fig. 12). Halpaap et al. (2018), observed an abrupt deepening of seismicity that occurs south of a line which connects MFZ with CHSZ, suggesting this as the COB boundary. An interesting feature of our model, revealed in the cross-sections of Fig. 13, is the reversal of the anomalies’ polarity above the Moho of the downgoing African Plate, north and south of the proposed COB; negative anomalies to the north (Fig. 13, a1–a2, b1–b2, c1–c2, d1–d2) and positive anomalies to the south (Fig. 13, f1–f2) support continental and oceanic crust, respectively. Negative anomalies are related to a shallower, gently dipping ‘continental’ WBZ and positive ones to a deeper, steeply dipping ‘oceanic’ WBZ. Furthermore, this is a region where the resolution of the model is high, as indicated by both reconstructed checkerboard models (Fig. 6) and high resolution matrix diagonal element (RDE) values (Figs S4 and S17b), which indicates that these results are reliable. The recent earthquake of Mw = 6.7 that took place southwest of Zakynthos on the 25th of October 2018 (Papadimitriou et al. 2019) has likely occurred at the northernmost part of the oceanic slab's interface with the upper plate, in a region delimited by the F1 and F2 vertical, transverse slab faults proposed by Sachpazi et al. (2016). The VSV model's horizontal resolution does not permit an exact positioning of the COB, but the above observations support a COB similar to the one suggested by Halpaap et al. (2018), that is a SW–NE trending line, subparallel to the MFZ, connecting it with N. Euboea (Fig. 14), likely coinciding with the F1 fault of Sachpazi et al. (2016). Figure 14. Open in new tabDownload slide Summary of the main constraints deduced for the geometry of the Hellenic lithosphere. Blue and red arrows denote horizontal compressional and extensional principal strain-rate axes, respectively (Kreemer et al. 2003). Epicentres of selected seismicity (focal depth ≥ 40 km) are from the herein used catalogues. For more placenames discussed in the text refer to Fig. 1. Figure 14. Open in new tabDownload slide Summary of the main constraints deduced for the geometry of the Hellenic lithosphere. Blue and red arrows denote horizontal compressional and extensional principal strain-rate axes, respectively (Kreemer et al. 2003). Epicentres of selected seismicity (focal depth ≥ 40 km) are from the herein used catalogues. For more placenames discussed in the text refer to Fig. 1. COB is critical for the regional tectonics, in particular concerning trench-perpendicular strike-slip faulting across CLAB, E–W normal faulting along the Corinth rift to its north turning into N–S normal faulting towards its south (Kapetanidis & Kassaras 2019). The inferred feature is compatible with the combined results of mantle anisotropy (Endrun et al. 2011; Evangelidis 2017; Kaviris et al. 2018) and crustal stress-states (Kapetanidis & Kassaras 2019). Specifically, the largest angular difference (75°–90°) between the SKS fast polarization directions and the maximum horizontal crustal stress axis, SHmax, is observed along COB, possibly indicating trench-parallel mantle-flow, correlated with a slab tear that induces strain in a perpendicular direction. A noticeable feature in the crustal stress model of Kapetanidis & Kassaras (2019) is a reversal of stress-shape (Φ) values south of the inferred interface, from uniaxial extension in the north to uniaxial compression in the south, possibly probing increase of slab-pull forces towards the south (Fig. S5a). 7.2 Slab segmentation The issue of STEP fault(s) across HSZ is thereinafter examined based on the new data. We note that, given the limited capacity of the attained 3-D VSV model to provide accurate constraints, implications on the matter are supported by subcrustal seismicity patterns. 7.2.1 Vertical tearing For the depth range between 40 and 140 km, our 3-D model, tentatively, does not show evidence of a trench-perpendicular vertical slab tear in W. Greece, in agreement with Pearce et al. (2012), Halpaap et al. (2018) and Hansen et al. (2019). Prolongation of MFZ to the CHSZ towards the western termination of NAT (Serpetsidaki et al. 2014) is also not indicated by the arrangement of velocities. Nonetheless, subcrustal seismicity exhibits a particularly interesting feature beneath central and NW Peloponnese, down to the depth of ∼60 km, cross-cutting with the Moho of the downgoing slab (cross-sections e1–e2, f1–f2, i1–i2 of Fig. 13). Halpaap et al. (2019) found a similar pattern of intermediate depth seismicity in this locality, interpreting it as the effect of upward migration of earthquake-inducing fluids through vents at the slab interface. This distribution of earthquake hypocentres shows brittle deformation along the COB, herein suggested to occur along a line connecting MFZ with NAT. Provided that existing models suggest little to no coupling between the crust and upper-mantle (Jackson & McKenzie 1988; Shaw & Jackson 2010), or that strong coupling is limited to the upper part of the crust (Laigle et al. 2002), clustering of subcrustal earthquakes in the vicinity of MFZ likely indicates brittle deformation at the whole crust (Halpaap et al. 2018) reaching upper-mantle through a STEP. However, this issue, of particular significance for the region's geodynamics and consequently seismic hazard is a subject of further research. A different structure is implied between the western and eastern part of the Hellenic Arc, manifested by a steeper and deeper WBZ that reaches the depth of ∼180 km in EHSZ in the area of Kos (cross-section l1–l2 in Fig. 13). A sharp stopping of the intermediate depth seismicity east of Rhodes indicates the easternmost termination of the EHSZ (Bocchini et al. 2018; Fig. 12). A shallower (25°) and less active WBZ is indicated in the area of Crete (Fig. 13, k1–k2) along with a steeper (45°) and more seismically active one in the area of Karpathos-Rhodes (Fig. 13, l1–l2). As seen in all cross-sections, seismic velocities roughly follow the trend of intermediate seismicity. Inversion of focal mechanisms (Rontogianni et al. 2011), indicates stress differentiation between the western and eastern part of HSZ at least for the shallow part of the WBZ (50–90 km depth), where horizontal compression and triaxial stress are observed, respectively. On the other hand, vertical compression indicated at depths of 90–180 km could account for increased slab-pull forces which in turn likely explain the larger dip of the eastern slab segment. The eastern termination of EHSZ is a region of apparent complexity, where a change of the slab geometry is indicated, separating EHSZ into a western and an eastern part. Cross-sections m1–m2, n1–n2 in Fig. 13 show a northward dipping WBZ at 30° in the west beneath Crete and another dipping NW at 40° in the east, beneath the area of Karpathos-Rhodes. While velocity amplitudes follow the trend of seismicity, we observe an LVZ beneath the WBZs (Fig. 13, m1–m2, n1–n2). Velocities correlate well with a hot mantle, explained as the result of convection from the asthenosphere through a N–S vertical tear separating the HSZ from the Western Cyprus Subduction Zone (e.g. Piromallo & Morelli 2003; van Hinsbergen et al. 2010; Biryol et al. 2011; Legendre et al. 2012; Govers & Fichtner 2016), or a horizontal tear (Faccenna et al. 2014). There is no indication for a horizontal slab tear or slab break-off (Faccenna et al. 2014) given the continuously distributed seismicity, in contrast to Papazachos et al. (2005), who observed lack of seismicity at depths between 90 and 140 km. Our compiled catalogue does not exhibit lack of seismicity at these depths, but on the contrary, it demonstrates that the deep part of the slab, between 100 and 150 km, hosts the major part of seismicity, especially at the easternmost part of the HSZ. Although the resolved velocities for this part of the 3-D model should be considered with caution, since the associated phase-velocity relative error is not negligible (6–9 per cent), the observed configuration complies with an LVZ that is present in most tomographic models in the region (e.g. Piromallo & Morelli 2003; Biryol et al. 2011; Salaün et al. 2012; Gessner et al. 2013), interpreted as asthenospheric flow through a first-order slab tear beneath SW Anatolia producing high-temperature metamorphism and felsic magmatism (Doglioni et al. 2002; Agostini et al. 2007; Dilek & Altunkaynak 2008; Salaün et al. 2012; Jolivet et al. 2015). It is moreover consistent with shear wave splitting measurements in EHSZ between Crete and Rhodes, interpreted as flow in the mantle wedge through a slab tear (Evangelidis 2017). Taking into account the above implications, a N–S vertical slab tear beneath SW Anatolia could likely explain the geometry of the eastern termination of the EHSZ without considering segmentation east of Crete as suggested by Bocchini et al. (2018). Asthenospheric flow along a N–S axis could explain the inferred similar slab orientation, while it likely facilitates slab's sink within a hot, buoyant mantle, increasing its dip, due to increased slab-pull forces (Rontogianni et al. 2011). 7.2.2 Horizontal tearing Since variations of intermediate depth seismicity can be the effect of slab-pull forces, lack of intermediate depth seismicity north of N39° (Bocchini et al. 2018) could be explained by diminishing of slab-pull forces related to a totally detached slab along a trench-parallel tear beneath mainland Greece that was observed in tomographic images at depths between ∼150 and 500 km (Spakman et al. 1988; Bijwaard et al. 1998; Koulakov et al. 2009; Zhu et al. 2015; Hansen et al. 2019). In this respect, deepening of seismicity towards the south could be justified by a partially detached southward migrating horizontal slab tear, whose southernmost extent is yet unclear. Zhu et al. (2015) suggest that the tear only affects northern Greece, whereas Spakman et al. (1988), Meijer & Wortel (1996) argue that it continues into southern Peloponnese. Meighan et al. (2013) notice that if one would expect a southward propagation of the tear, it should correlate with a simultaneous migration of shallow deformation related to earthquake clusters and swarms. Interestingly, the NE tip of MFZ marks an N–S alignment of abundant shallow seismicity that crosses the whole Peloponnese from the western Corinth gulf in the north, down to Messinia (e.g. Kassaras et al. 2014) in the south (Fig. 12). Southward propagation of the tear could be the cause of strong tectonic instability in the region of CLAB, marked by the occurrence of several moderate-to-strong earthquakes, as well as swarms, which have been monitored in detail by the HUSN and temporary deployments during the last decade (e.g. Kassaras et al. 2014; Kapetanidis et al. 2015). 7.3 The lithosphere–asthenosphere boundary (LAB) Variations in the topography of the LAB provide an essential constraint for understanding mantle deformation in response to tectonic forces, volcanism processes and vice versa (e.g. Fischer et al. 2010). The lithospheric thickness beneath the Hellenic region and W. Anatolia has been measured by teleseismic P- and S-receiver functions (Sodoudi et al. 2006, 2015; Kind et al. 2015). For the Greek region, Sodoudi et al. (2006) found the LAB of the subducted African Plate at 100 km depth beneath Crete and S. Aegean and at about 225 km depth beneath SAAVA, whereas the LAB beneath mainland Greece is indicated near 150 km depth. Lithospheric thickness of the subducted African Plate could not be herein identified, given the maximum depth (140 km) and the limited resolution of the 3-D model. However, the new 3-D model provides some interesting implications on the LAB of the Aegean and Anatolian plates. Low velocity patterns are systematically indicated in the Aegean and Anatolia upper plates that could likely be attributed to the asthenospheric interface and/or asthenospheric bodies within the upper mantle, which are linked to crustal kinematics, heat-flow and volcanism. A pattern widely seen in Fig. 13 is that of low velocities of the overriding plate(s) at shallow depths throughout the region, indicating a thin lithosphere, of the order of ∼50 km. Although this finding should be taken with caution since low period dispersion data are not available, it is in contradiction with the 150 km depth suggested by Sodoudi et al. (2006) for continental Greece, but is in agreement with Kind et al. (2015), who suggest a laminated mantle lithospheric structure beneath the broader region of eastern Mediterranean, whose bottom consists of several neighbouring discontinuities of decreasing velocity with depth, and with a vertical extent of the laminated region of up to 50 km. Moreover, it is consistent with a shallow LAB found beneath Anatolia at an average depth of about 70 km by receiver functions (Sodoudi et al. 2006, 2015; Kind et al. 2015), attenuation observations (Gök et al. 2003), results obtained from body and surface wave tomography studies (Biryol et al. 2011; Bakirci et al. 2012; Salaün et al. 2012; Fichtner et al. 2013) and body-wave velocity anomalies (e.g. Al-Lazki et al. 2004; Maggi & Priestley 2005; Gök et al. 2007) indicative for lithospheric removal in the region (Kind et al. 2015). Such a thin lithosphere/laminated mantle lithospheric structure beneath Aegean/W. Anatolia is relevant to asthenospheric convection (e.g. Jolivet et al. 2013), as pointed out also by the large amplitude of mantle shear- (Hatzfeld et al. 2001; Evangelidis 2017; Kaviris et al. 2018) and Rayleigh-wave anisotropy (Endrun et al. 2011).The existence of a shallow LAB/laminated mantle affects volcanism and increased heat-flow which is widely seen throughout the Aegean and W. Anatolia (Pe-Piper & Piper 2002; Serpen et al. 2010). More specifically, shallow, low velocity anomalies in NW Greece (cross-sections a1–a2, b1–b2, c1–c2 of Fig. 13) are correlated with low-enthalpy (e.g. Paiko volcano, Andritsos et al. 2011) and remnant volcanism (thermal springs) in N. Aegean and NW Anatolia (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015). Southerly, the non-deforming central Aegean lithosphere displays higher seismic velocities with respect to the deforming lithosphere of N. Aegean. This pattern, particularly evidenced by Rayleigh-wave anisotropy, is explained by viscosity contrast between the two regions due to temperature and/or fluids content alteration that leads to modification of the lithosphere (Endrun et al. 2011). SAAVA, a prominent feature in our tomographic images in S. Aegean, exhibits fast anomalies in its western part and slow velocities in the eastern part (Fig. 11). The latter is compatible with different eruptive behaviour and magmatic compositions between the two segments (Pe-Piper & Piper 2002; Francalanci et al. 2005) and different stress states (Kapetanidis & Kassaras 2019). Effusive volcanism in the western part has been explained by the presence of water-rich fluids or melts associated with high Vp/Vs ratios (Halpaap et al. 2018), a suggestion consistent with seismicity of the WBZ of WHSZ, the largest part of which occurs at a shallower depth (80–100 km) with respect to the EHSZ (100–150 km). The abovementioned observations could be linked to variation of lithospheric thickness beneath the western and eastern segment of SAAVA, as it is observed from the cross-sections of Fig. 13, where LAB (or laminated upper mantle) appears to be shallower beneath the eastern segment (<50 km depth) with respect to the western one (∼50 km depth). 8 CONCLUSIONS In this work, we have investigated the shear wave structure of the lower part of the Hellenic lithosphere (>40 km depth) using surface wave tomography. We have studied phase-velocity dispersion by analysing several thousand fundamental-mode teleseismic Rayleigh wave recordings using a cross-correlation scheme that takes into account multipath/backscattering effects. The inversion of phase-velocities yielded a 3-D shear wave velocity model for the region down to the depth of 140 km. The main outcome of this work can be summarized in the following: Surface wave analysis has proven sufficient in exhibiting and discriminating the gross geometrical characteristics of lithospheric plates in complex convergence zones. The array-based methodology is found to be robust, yielding a better resolution with respect to previous surface wave studies based on the two-station approach. Phase-velocity maps have been constrained with a varying resolution, the best of which has been achieved in the area of central Greece. Both single-event and stacked Eikonal and Helmholtz phase-velocity maps exhibit consistent features which in turn are relevant to gross geological structures of the lithosphere, such as the slab beneath the Hellenides in N. Greece and the slab in the HSZ extending from northern Peloponnese to the Turkish coast, the mantle wedge above the slab along the forearc and a hot/buoyant upper mantle beneath the Aegean and W. Anatolia related to increased heat flow. Lack of low velocities south of SW Anatolia supports that vertical tear is limited to regions beneath continental Anatolia. Both continental and oceanic subduction occur along the HSZ, as particularly supported by both phase velocity maps and the 3-D shear velocity model. The COB is suggested to exist along a trench-perpendicular line that connects MFZ with N. Euboea (Fig. 14). Resolution restraints of the 3-D VSV model do not permit the discrimination of trench-perpendicular vertical tear of the slab in the WHSZ within the resolved depth range 40–140 km. Nonetheless, reliable intermediate seismicity in the region of N. Peloponnese, tentatively supports a vertical lithospheric tear (STEP fault) along the inferred COB. Resolved seismic velocities indicate a different structure between the western and eastern part of the Hellenic Arc, manifested in particular by the steeper and deeper WBZ in the EHSZ, dipping to the NW. This differentiation can been likely explained by a NW vertical slab tear restricted onshore SW Anatolia as phase- and shear-velocity distribution implies. The thickness of the subducted African lithosphere cannot be resolved given the limited depth of the model and the low resolution. Implications of shear velocity anomalies for a shallow LAB and/or a laminated mantle in the overriding plate(s) are in-line with the spatial distribution of volcanism (present and remnant) throughout the region. SUPPORTING INFORMATION Figure S1. Map of epicentres of the four seismic catalogues used in this study. Figure S2. Statistics of the NKUA-SL seismic catalogue. (a) RMS error, (b) Location errors in the three dimensions, (c) azimuthal gap and (d) focal depth. Figure S3. Left-hand panel: map of crustal thickness (CRUST1.0 model), Right-hand panel: map of sediments thickness (Laske & Masters 1997; https://igppweb.ucsd.edu/∼gabi/crust1.html). Figure S4. Vertical resolution of the 3-D model, in terms of resolution matrix diagonal element (RDE) values. Best resolution is achieved between 40 and 100 km depth. Figure S5. (a) Distribution of SHmax axes and stress shape at crustal depths in the broader area of Greece (modified after Kapetanidis & Kassaras 2019). The area is divided in 8 subregions which are used in panel (b), Fig. S6 and Tables S1 and S2. (b) Average 1-D VSV models for the 8 subregions of panel (a). Figure S6. Average dispersion curves for the eight subregions of Fig. S5(a). Figure S7. Cross-sections along the profile lines of Fig. 12 (in the paper) as velocity perturbations relative to the mean VSV value of the respective depth. Purple contour lines represent the relative error of phase-velocities, translated at depth. Seismicity (this study, Karakonstantis 2017; Bocchini et al. 2018) as well as volcanoes and thermal springs (Siebert & Simkin 2002; Serpen et al. 2010; Erkan 2015) are also plotted. Black dashed lines denote the Moho interface of the subducted lithosphere (after Halpaap et al. 2018). The white dashed line corresponds to the CRUST1.0 model of the upper plate Moho surface. VSV velocities are smoothened horizontally within circles of radius 100 km on each grid cell over the whole map, prior to the generation of cross-sections. See also Fig. 13 for cross-sections of relative VSV perturbations. Figure S8. Maps indicating the effect of the amplitude correction term on the apparent phase velocities for T = 40 s (top panel) and 70 s (bottom panel). Left-hand panel: the differential distribution between apparent and structural phase velocities and right-hand panel: the respective amplitude correction term applied to the Eikonal tomography. Figure S9. Map of epicentres of the compiled seismicity catalogue (a) for all focal depths (same as Fig. S1), (b) for depths ≥ 40 km (as in Figs 12 and 14), with different colours representing the four seismic catalogues used in the present study. (1) Grey circles: NKUA-SL manual analysis 2012–2017, (2) EGELADOS catalogue (Friederich & Meier 2008; Brüstle 2012): red circles, (3) LIBNET catalogue (Becker et al. 2010): blue circles and (4) CYCNET catalogue (Bohnhoff et al. 2004, 2006; Brüstle 2012): green circles. Figure S10. Maps of average amplitude ratio used to correct biases before the calculation of the correction term, (a) for T = 40 s, (b) for T = 60 s. Triangles show measurements at individual stations while the background map shows the respective minimum curvature surface fit, similar to the one applied on the measured amplitudes. Figure S11. Trade-off curves between model variance and data misfit for the Eikonal maps, used to select appropriate smoothing parameters (value in red) for each period. Figure S12. Single-event measurements at a period of T = 60 s for an Mw = 6.1 earthquake that occurred on 11 May 2013 02:08:08 UTC in S. Iran, (a) phase-delays relative to one of the stations (zero relative travel-time) with valid measurements, (b) propagation direction, (c) differences between propagation direction and great-circle path and (d) apparent velocity map through the Eikonal equation. Contour spacing corresponds to ½ of the respective wavelength while vectors denote the propagation direction (normal to the contours). Figure S13. Single-event measurements at a period of T = 60 s for the same event as in Fig. S12. (a, d) Amplitude minimum curvature surface fit, (b, e) amplitude correction term before smoothing, (c, f) amplitude correction term after smoothing, (h, i) structural phase-velocity (Helmholtz) map, (a, b, c, h) before and (d, e, f, i) after correction for amplitude-ratio biases. The respective Eikonal map (g; same as Fig. S12d) is displayed for comparison. Figure S14. Single-event measurements at a period of T = 80 s for an Mw = 6.8 earthquake that occurred on 12 May 2015 21:12:59 UTC in Japan, (a) phase-delays relative to one of the stations (zero relative travel-time) with valid measurements, (b) propagation direction, (c) differences between propagation direction and great-circle path and (d) apparent velocity map through the Eikonal equation. Contour spacing corresponds to ½ of the respective wavelength while vectors denote the propagation direction (normal to the contours). Figure S15. Single-event measurements at a period of T = 80 s for the same event as in Fig. S14. (a, d) amplitude minimum curvature surface fit, (b, e) amplitude correction term before smoothing, (c, f) amplitude correction term after smoothing, (h, i) structural phase-velocity (Helmholtz) map, (a, b, c, h) before and (d, e, f, i) after correction for amplitude-ratio biases. The respective Eikonal map (g; same as Fig. S14d) is displayed for comparison. Figure S16. Diagrams showing the azimuthal distribution of mean deviations between local propagation directions estimated from the apparent phase gradients and the epicentre-station great-circle path for T = 30 s (top panel) and 70 s (bottom panel). Figure S17. Right-hand panel: resolution matrices for selected nodes a–d. Left-hand panel: marked locations of nodes a–d on the map over the reconstruction of checkerboard pattern for T = 40 s. Figure S18. VSV resolution kernels for different layers at node ‘a’ of Fig. S17. Figure S19. VSV resolution kernels for different layers at node ‘b’ of Fig. S17. Figure S20. VSV resolution kernels for different layers at node ‘c’ of Fig. S17. Figure S21. VSV resolution kernels for different layers at node ‘d’ of Fig. S17. Table S1. Average 1-D VSV models for the subregions of Fig. S5(a). Table S2. Average phase velocity per period for the subregions of Fig. S5(a). Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. ACKNOWLEDGEMENTS We acknowledge support of this study by the project ‘HELPOS—Hellenic Plate Observing System’ (MIS 5002697) which is implemented under the Action ‘Reinforcement of the Research and Innovation Infrastructure’, funded by the Operational Programme ‘Competitiveness, Entrepreneurship and Innovation’ (NSRF 2014–2020) and co-financed by Greece and the European Union (European Regional Development Fund). We thank EIDA for providing data from regional networks HL (doi: 10.7914/SN/HL), HT (doi: 10.7914/SN/HT), HA (doi: 10.7914/SN/HA), HP (doi: 10.7914/SN/HP), HC (doi: 10.7914/SN/HC), CQ (doi: 10.7914/SN/CQ), Z3 (doi: 10.14470/M87550267382), and CL (doi: 10.15778/RESIF.CL). We thank Dr G. Kaviris and I. Spingos for their help with collecting the teleseismic earthquake data and Dr G.-M. Bocchini for making available seismicity catalogue data. We are greatly indebted to Ge Jin for clarifications on the ASWSM suite, N. Accardo for providing the checkerboard test plugin, A. Ganas for discussions. Maps and cross-sections were drawn using the Generic Mapping Tools (GMT) software (Wessel & Smith 1998). We would like to thank Ge Jin and the other two anonymous reviewers, whose constructive comments have greatly helped to improve the quality and clarity of this paper. Authors Contributions: I.K. conceived the idea for this study, I.K. and V.K. collected the data, performed data processing and analysis, A.K. provided the NKUA-SL seismicity catalogue, P.P. supervised the project, all co-authors discussed the results and contributed to the writing of the paper. REFERENCES Accardo N.J. et al. . , 2017 . Surface wave imaging of the weakly extended Malawi Rift from ambient-noise and teleseismic Rayleigh waves from onshore and lake-bottom seismometers , Geophys. J. Int. , 209 , 1892 – 1905 . 10.1093/gji/ggx133 Google Scholar Crossref Search ADS WorldCat Crossref Agostini S. , Doglioni C. , Innocenti F. , Manetti P. , Tonarini S. , Savascin M.Y. , 2007 . The transition from subduction-related to intraplate Neogene magmatism in the Western Anatolia and Aegean area , in Cenozoic Volcanism in the Mediterranean Area , Vol. 418 , pp. 1 – 16 ., eds Beccaluva L. , Bianchini G. , Wilson M. , Geological Society of America Special Paper . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Al-Lazki A. I. , Seber D. , Sandvol E. , Turkelli N. , Mohamad R. , Barazangi M. , 2004 . Tomographic Pn velocity and anisotropy structure beneath the Anatolian Plateau (eastern Turkey) and the surrounding regions , Geophys. Res. Lett. , 30 , 8043 . OpenURL Placeholder Text WorldCat Andritsos N. , Dalabakis P. , Karydakis G. , Kolios N. , Fytikas M. , 2011 . Characteristics of low-enthalpy geothermal applications in Greece , Renew. Ener. , 36 , 1298 – 1305 . 10.1016/j.renene.2010.10.008 Google Scholar Crossref Search ADS WorldCat Crossref Armijo R. , Meyer B. , King G. , Rigo A. , Papanastassiou D. , 1996 . Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean , Geophys. J. Int. , 126 , 11 – 53 . 10.1111/j.1365-246X.1996.tb05264.x Google Scholar Crossref Search ADS WorldCat Crossref Avallone A. , Briole P. , Agatza-Balodimou A. , 2004 . Analysis of eleven years of deformation measured by GPS in the Corinth Rift laboratory area . C. R. Acad. Sci. Geosci. , 336 , 301 – 311 . 10.1016/j.crte.2003.12.007 Google Scholar Crossref Search ADS WorldCat Crossref Bakirci T. , Yoshizawa K. , Özer M.F. , 2012 . Three-dimensional S wave structure of the upper mantle beneath Turkey from surface wave tomography , Geophys. J. Int. , 190 , 1058 – 1076 . 10.1111/j.1365-246X.2012.05526.x Google Scholar Crossref Search ADS WorldCat Crossref Becker D. , Meier T. , Bohnhoff M. , Harjes H.P. , 2010 . Seismicity at the convergent plate boundary offshore Crete, Greece, observed by an amphibian network , J. Seismol. , 14 , 369 – 392 . 10.1007/s10950-009-9170-2 Google Scholar Crossref Search ADS WorldCat Crossref Biryol C.B. , Beck S.L. , Zandt G. , Özacar A.A. , 2011 . Segmented African lithosphere beneath the Anatolian region inferred from teleseismic P-wave tomography , Geophys. J. Int. , 184 , 1037 – 1057 . 10.1111/j.1365-246X.2010.04910.x Google Scholar Crossref Search ADS WorldCat Crossref Bijwaard H. , Spakman W. , Engdahl E.R. , 1998 . Closing the gap between regional and global travel time tomography , J. geophys. Res. , 103 , 30 055 – 30 078 . 10.1029/98JB02467 Google Scholar Crossref Search ADS WorldCat Crossref Bocchini G.M. et al. . , 2018 . Tearing, segmentation, and backstepping of subduction in the Aegean: new insights from seismicity , Tectonophysics , 734–735 , 96 – 118 . 10.1016/j.tecto.2018.04.002 Google Scholar Crossref Search ADS WorldCat Crossref Bodin T. , Maupin V. , 2008 . Resolution potential of surface wave phase velocity measurements at small arrays , Geophys. J. Int. , 172 , 698 – 706 . 10.1111/j.1365-246X.2007.03668.x Google Scholar Crossref Search ADS WorldCat Crossref Bohnhoff M. , Makris J. , Papanikolaou D. , Stavrakakis G. , 2001 . Crustal investigations of the Hellenic subduction zone using wide aperture seismic data , Tectonophysics , 343 , 239 – 262 . 10.1016/S0040-1951(01)00264-5 Google Scholar Crossref Search ADS WorldCat Crossref Bohnhoff M. , Rische M. , Meier T. , Becker D. , Stavrakakis G. , Harjes H.P. , 2006 . Microseismic activity in the Hellenic Volcanic Arc, Greece, with emphasis on the seismotectonic setting of the Santorini–Amorgos zone , Tectonophysics , 423 , 17 – 33 . 10.1016/j.tecto.2006.03.024 Google Scholar Crossref Search ADS WorldCat Crossref Bohnhoff M. , Rische M. , Meier T. , Endrun B. , Becker D. , Harjes H.P. , Stavrakakis G. , 2004 . CYC-NET: a temporary seismic network on the Cyclades (Aegean Sea, Greece) , Seismol. Res. Lett. , 75 , 352 – 359 . 10.1785/gssrl.75.3.352 Google Scholar Crossref Search ADS WorldCat Crossref Bourova E. , Kassaras I. , Pedersen H.A. , Yanovskaya T. , Hatzfeld D. , Kiratzi A. , 2005 . Constraints on absolute S velocities beneath the Aegean Sea from surface wave analysis , Geophys. J. Int. , 160 , 1006 – 1019 . 10.1111/j.1365-246X.2005.02565.x Google Scholar Crossref Search ADS WorldCat Crossref Briole P. et al. . , 2000 . Active deformation of the Corinth rift Greece: Results from repeated Global Positoning System surveys between 1990 and 1995 , J. geophys. Res. , 105 , 25 605 – 25 625 . 10.1029/2000JB900148 Google Scholar Crossref Search ADS WorldCat Crossref Brüstle A. , 2012 , Seismicity of the eastern Hellenic subduction zone , PhD thesis , University of Bochum , Germany . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Channell J.E.T. , Kozur H.W. , 1997 . How many oceans? Meliata, Vardar and Pindos oceans in Mesozoic Alpine paleogeography , Geology , 25 , 183 . 10.1130/0091-7613(1997)025%3c0183:HMOMVA%3e2.3.CO;2 Google Scholar Crossref Search ADS WorldCat Crossref Chousianitis K. , Ganas A. , Evangelidis C.P. , 2015 . Strain and rotation rate patterns of mainland Greece from continuous GPS data and comparison between seismic and geodetic moment release , J. geophys. Res. , 120 , 3909 – 3931 . 10.1002/2014JB011762 Google Scholar Crossref Search ADS WorldCat Crossref Cianetti S. , Gasperini P. , Giunchi C. , Boschi E. , 2001 . Numerical modelling of the Aegean-Anatolian region: geodynamical constraints from observed rheological heterogeneities , Geophys. J. Int. , 146 , 760 – 780 . 10.1046/j.1365-246X.2001.00492.x Google Scholar Crossref Search ADS WorldCat Crossref Clément C. , Sachpazi M. , Charvis P. , Graindorge D. , Laigle M. , Hirn A. , Zafiropoulos G. , 2004 . Reflection-refraction seismics in the Gulf of Corinth: hints at deep structure and control of the deep marine basin , Tectonophysics , 391 , 85 – 95 . 10.1016/j.tecto.2004.07.010 Google Scholar Crossref Search ADS WorldCat Crossref Dercourt J. , Ricou L. E. , Vrielynck B. , eds., 1993 . Atlas Tethys Palaeoenvironmental Maps , Gauthier-Villars , 260 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Di Luccio F. , Pasyanos M. , 2007 . Crustal and upper-mantle structure in the eastern Mediterranean from the analysis of surface wave dispersion curves , Geophys. J. Int. , 169 , 1139 – 1152 . 10.1111/j.1365-246X.2007.03332.x Google Scholar Crossref Search ADS WorldCat Crossref Dilek Y. , Altunkaynak S. , 2008 . Geochemical and temporal evolution of Cenozoic magmatism in western Turkey: mantle response to collision, slab breakoff, and lithospheric tearing in an orogenic belt , in Collision and Collapse at the Africa-Arabia-Eurasia Subduction Zone , Vol. 311 , pp. 213 – 233 ., eds Van Hinsbergen D.J.J. , Edwards M.A. , Govers R. , Geological Society, London Special Publications . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Doglioni C. , Agostini S. , Crespi M. , Innocenti F. , Manetti P. , Riguzzi F. , Savascin Y. , 2002 . On the extension in western Anatolia and the Aegean Sea , in Reconstruction of the Evolution of the Alpine-Himalayan Orogen , Vol. 8 , pp. 169 – 183 ., eds Rosenbaum G. , Lister G.S. , Journal of the Virtual Explorer . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Durand V. , Bouchon M. , Floyd M.A. , Theodulidis N. , Marsan D. , Karabulut H. , Schmittbuhl J. , 2014 . Observation of the spread of slow deformation in Greece following the breakup of the slab , Geophys. Res. Lett. , 41 , 7129 – 7134 . 10.1002/2014GL061408 Google Scholar Crossref Search ADS WorldCat Crossref Erkan K. , 2015 . Geothermal investigations in western Anatolia using equilibrium temperatures from shallow boreholes , Solid Earth , 6 , 103 – 113 . 10.5194/se-6-103-2015 Google Scholar Crossref Search ADS WorldCat Crossref Endrun B. , Lebedev S. , Meier T. , Tirel C. , Friederich W. , 2011 . Complex layered deformation within the Aegean crust and mantle revealed by seismic anisotropy , Nat. Geosci , 4 , 203 – 207 . 10.1038/ngeo1065 Google Scholar Crossref Search ADS WorldCat Crossref Endrun B. , Meier T. , Lebedev S. , Bohnhoff M. , Stavrakakis G. , Harjes H. P. , 2008 . S velocity structure and radial anisotropy in the Aegean region from surface wave dispersion , Geophys. J. Int. , 174 , 593 – 616 . 10.1111/j.1365-246X.2008.03802.x Google Scholar Crossref Search ADS WorldCat Crossref Evangelidis C.P. , 2017 . Seismic anisotropy in the Hellenic subduction zone: effects of slab segmentation and subslab mantle flow , Earth planet. Sci. Lett. , 480 , 97 – 106 . 10.1016/j.epsl.2017.10.003 Google Scholar Crossref Search ADS WorldCat Crossref Faccenna C. et al. . , 2014 . Mantle dynamics in the Mediterranean , Rev. Geophys. , 52 , 283 – 332 . 10.1002/2013RG000444 Google Scholar Crossref Search ADS WorldCat Crossref Faccenna C. , Jolivet L. , Piromallo C. , Morelli A. , 2003 . Subduction and the depth of convection in the Mediterranean mantle , J. geophys. Res. , 108 , doi:10.1029/2001JB001690 . 10.1029/2001JB001690 OpenURL Placeholder Text WorldCat Crossref Fichtner A. , Saygin E. , Taymaz T. , Cupillard P. , Capdeville Y. , Trampert J. , 2013 . The deep structure of the North Anatolian Fault Zone , Earth planet. Sc. Lett. , 373 , 109 – 117 . 10.1016/j.epsl.2013.04.027 Google Scholar Crossref Search ADS WorldCat Crossref Finetti I. , Del Ben A. , 2005 . Crustal Tectono-Stratigraphy of the ionian sea from new integrated CROP seismic data , in CROP Project: Deep Seismic Exploration of the Central Mediterranean and Italy , pp. 447 – 470 ., ed. Finetti I.R. , Elsevier . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Fischer K.M. , Ford H.A. , Abt D.L. , Rychert C.A. , 2010 . The lithosphere asthenosphere boundary , Annu. Rev. Earth planet. Sci. , 38 , 551 – 575 . 10.1146/annurev-earth-040809-152438 Google Scholar Crossref Search ADS WorldCat Crossref Foti S. et al. . , 2018 . Guidelines for the good practice of surface wave analysis: a product of the InterPACIFIC project , Bull. Earthq. Eng. , 16 , 2367 – 2420 . 10.1007/s10518-017-0206-7 Google Scholar Crossref Search ADS WorldCat Crossref Francalanci L. , Vougioukalakis G.E. , Perini G. , Manetti P. , 2005 . A West-East Traverse along the magmatism of the south Aegean volcanic arc in the light of volcanological, chemical and isotope data , Dev. Volcanol. , 7 , 65 – 111 . 10.1016/S1871-644X(05)80033-6 Google Scholar Crossref Search ADS WorldCat Crossref Friederich W. , Meier T. , 2008 . Temporary seismic broadband network acquired data on Hellenic subduction zone , Eos , 89 , 378 – 378 . 10.1029/2008EO400002 Google Scholar Crossref Search ADS WorldCat Crossref Fytikas M. , Innocenti F. , Manetti P. , Mazzuoli R. , Peccerillo A. , Villari L. , 1985 . Tertiary to Quaternary evolution of volcanism in the Aegean region: the Geological Evolution of the Eastern Mediterranean , Spec. Publ. Geol. Soc. , 17 , 687 – 699 . 10.1144/GSL.SP.1984.017.01.55 Google Scholar Crossref Search ADS WorldCat Crossref Ganas A. , Serpelloni E. , Drakatos G. , Kolligri M. , Adamis I. , Tsimi C. , Batsi E. , 2009 . The Mw 6.4 SW-Achaia (Western Greece) earthquake of 8 June 2008: seismological, field, GPS observations, and stress modeling , J. Earthq. Eng. , 13 , 1101 – 1124 . 10.1080/13632460902933899 Google Scholar Crossref Search ADS WorldCat Crossref Gee L.S. , Jordan T.H. , 1992 . Generalized seismological data functionals , Geophys. J. Int. , 111 , 363 – 390 . 10.1111/j.1365-246X.1992.tb00584.x Google Scholar Crossref Search ADS WorldCat Crossref Gesret A. , Laigle M. , Diaz J. , Sachpazi M. , Charalampakis M. , Hirn A. , 2011 . Slab top dips resolved by teleseismic converted waves in the Hellenic subduction zone , Geophys. Res. Lett. , 38 , L20304 . 10.1029/2011GL048996 Google Scholar Crossref Search ADS WorldCat Crossref Gessner K. , Gallardo L.A. , Markwitz V. , Ring U. , Thomson S.N. , 2013 . What caused the denudation of the Menderes Massif: review of crustal evolution, lithosphere structure, and dynamic topography in southwest Turkey , Gondwana Res. , 24 , 243 – 274 . 10.1016/j.gr.2013.01.005 Google Scholar Crossref Search ADS WorldCat Crossref Gök R. , Sandvol E. , Törkelli N. , Seber D. , Barazangi M. , 2003 . Sn attenuation in the Anatolian and Iranian plateau and surrounding regions , Geophys. Res. Lett. , 30 , 8042 . 10.1029/2003GL018020 OpenURL Placeholder Text WorldCat Crossref Gök R. , Pasyanos M. , Zor E. , 2007 . Lithospheric structure of the continent-continent collision zone: eastern Turkey , Geophys. J. Int. , 169 , 1079 – 1088 . 10.1111/j.1365-246X.2006.03288.x Google Scholar Crossref Search ADS WorldCat Crossref Goldsworthy M. , Jackson J. , Haines A.J. , 2002 . The continuity of active fault systems in Greece , Geophys. J. Int. , 148 , 596 – 618 . 10.1046/j.1365-246X.2002.01609.x Google Scholar Crossref Search ADS WorldCat Crossref Govers R. , Fichtner A. , 2016 . Signature of slab fragmentation beneath Anatolia from full-waveform tomography , Earth planet. Sci. Lett. , 450 , 10 – 19 . 10.1016/j.epsl.2016.06.014 Google Scholar Crossref Search ADS WorldCat Crossref Govers R. , Wortel M.J.R. , 2005 . Lithosphere tearing at STEP faults: response to edges of subduction zones , Earth planet. Sci. Lett. , 236 , 505 – 523 . 10.1016/j.epsl.2005.03.022 Google Scholar Crossref Search ADS WorldCat Crossref Halpaap F. , Rondenay S. , Perrin A. , Goes S. , Ottemöller L. , Austrheim H. , Shaw R. , Eeken T. , 2019 . Earthquakes track subduction fluids from slab source to mantle wedge sink , Sci. Adv. , 5 , eaav7369 . 10.1126/sciadv.aav7369 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Halpaap F. , Rondenay S. , Ottemöller L. , 2018 . Seismicity, deformation and metamorphism in the western Hellenic subduction zone: new constraints from tomography , J. geophys. Res. , 123 , 3000 – 3026 . 10.1002/2017JB015154 Google Scholar Crossref Search ADS WorldCat Crossref Hansen S.E. , Evangelidis C.P. , Papadopoulos G.A. , 2019 . Imaging slab detachment within the western Hellenic subduction zone , Geochem. Geophys. Geosyst. , 20 , doi:10.1029/2018GC007810 . 10.1029/2018GC007810 OpenURL Placeholder Text WorldCat Crossref Hatzfeld D. et al. . , 2001 . Shear wave anisotropy in the upper mantle beneath the Aegean related to internal deformation , J. geophys. Res. , 106 , 30 737 – 30 753 . 10.1029/2001JB000387 Google Scholar Crossref Search ADS WorldCat Crossref Hearn T.M. , 1999 . Uppermost mantle velocities and anisotropy beneath Europe , J. geophys. Res. , 104 , 15 123 – 15 139 . Google Scholar Crossref Search ADS WorldCat Herrmann R. , 2013 . Computer programs in seismology: an evolving tool for instruction and research , Seismol. Res. Lett. , 84 , 1081 – 1088 . 10.1029/1998JB900088 10.1785/0220110096 Google Scholar Crossref Search ADS WorldCat Crossref Ho T. , Priestley K. , Debayle E. , 2016 . A global horizontal shear velocity model of the upper mantle from multimode Love wave measurements , Geophys. J. Int. , 207 , 542 – 561 . 10.1093/gji/ggw292 Google Scholar Crossref Search ADS WorldCat Crossref Hosa A.M. , 2008 , Imaging of the Hellenic subduction zone by seismic tomography , BSc thesis , Massachusetts Institute of Technology , Cambridge . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Jackson J. , McKenzie D.P. , 1988 . The relationship between plate motions and seismic moment tensors, and the rates of active deformation in the Mediterranean and Middle East , Geophys. J. , 93 , 45 – 73 . 10.1111/j.1365-246X.1988.tb01387.x Google Scholar Crossref Search ADS WorldCat Crossref Jin G. , Gaherty J.B. , 2015 . Surface wave phase-velocity tomography based on multichannel cross-correlation , Geophys. J. Int. , 201 , 1383 – 1398 . 10.1093/gji/ggv079 Google Scholar Crossref Search ADS WorldCat Crossref Jin G. , Gaherty J.B. , Abers G.A. , Kim Y. , Eilon Z. , Buck W.R. , 2015 . Crust and upper mantle structure associated with extension in the Woodlark Rift, Papua New Guinea from Rayleigh-wave tomography , Geochem., Geophys. Geosyst. , 16 , 3808 – 3824 . 10.1002/2015GC005840 Google Scholar Crossref Search ADS WorldCat Crossref Jolivet L. , Brun J.P. , 2010 . Cenozoic geodynamic evolution of the Aegean , Int. J. Earth Sci. , 99 , 109 – 138 . 10.1007/s00531-008-0366-4 Google Scholar Crossref Search ADS WorldCat Crossref Jolivet L. et al. . , 2013 . Aegean tectonics: strain localisation, slab tearing and trench retreat , Tectonophysics , 597–598 , 1 – 33 . 10.1016/j.tecto.2012.06.011 Google Scholar Crossref Search ADS WorldCat Crossref Jolivet L. et al. . , 2015 . The geological signature of a slab tear below the Aegean , Tectonophysics . 659 , 166 – 182 . 10.1016/j.tecto.2015.08.004 Google Scholar Crossref Search ADS WorldCat Crossref Jolivet L. , Rimmelé G. , Oberhänsli R. , Goffé B. , Candan O. , 2004 . Correlation of synorogenic tectonic and metamorphic events in the Cyclades, the Lycian Nappes and the Menderes Massif. Geodynamic implications , Bull. Soc. Geol. Fr. , 175 , 217 – 238 . 10.2113/175.3.217 Google Scholar Crossref Search ADS WorldCat Crossref Kapetanidis V. , Kassaras I. , 2019 . Contemporary crustal stress of the Greek region deduced from earthquake focal mechanisms , J. Geodyn. , 123 , 55 – 82 . 10.1016/j.jog.2018.11.004 Google Scholar Crossref Search ADS WorldCat Crossref Kapetanidis V. et al. . , 2015 . The 2013 earthquake swarm in Helike, Greece: seismic activity at the root of old normal faults , Geophys. J. Int. , 202 , 2044 – 2073 . 10.1093/gji/ggv249 Google Scholar Crossref Search ADS WorldCat Crossref Karagianni E.E. , Papazachos C.B. , 2007 . Shear velocity structure in the Aegean region obtained by joint inversion of Rayleigh and Love waves , in The Geodynamics of the Aegean and Anatolia , Vol. 291 , pp. 159 – 181 ., eds Taymaz T. , Yilmaz Y. , Dilek Y. , Geological Society, London Special Publications . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Karagianni E.E. , Papazachos C.B. , Panagiotopoulos D.G. , Suhadolc P. , Vuan A. , Panza G.F. , 2005 . Shear velocity structure in the Aegean area obtained by inversion of Rayleigh-waves , Geophys. J. Int. , 160 , 127 – 143 . 10.1111/j.1365-246X.2005.02354.x Google Scholar Crossref Search ADS WorldCat Crossref Karakonstantis A. , 2017 , 3-D simulation of crust and upper mantle structure in the broader Hellenic area through seismic tomography , PhD dissertation , National and Kapodistrian University of Athens (in Greek) . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Karakonstantis A. , Papadimitriou P. , 2010 . Earthquake relocation in Greece using a unified and homogenized seismological catalogue , Bull. Geol. Soc. Greece , 43 , 2043 – 2052 . 10.12681/bgsg.11394 Google Scholar Crossref Search ADS WorldCat Crossref Kassaras I. , Makropoulos K. , Bourova E. , Pedersen H. , Hatzfeld D. , 2005 . Upper mantle structure of the Aegean derived from two-station phase velocities of fundamental mode Rayleigh waves , Dev. Volcanol. , 7 , 19 – 45 . 10.1016/S1871-644X(05)80031-2 Google Scholar Crossref Search ADS WorldCat Crossref Kassaras I. , Louis F. , Makropoulos K. , Magganas A. , Hatzfeld D. , 2009 . Elastic-Anelastic properties of the aegean Lithosphere-Asthenosphere inferred from long period Rayleigh-waves , in The Lithosphere: Geochemistry, Geology and Geophysics , pp. 383 , eds Anderson J.E. , Coates R.W. , Nova Publishers . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kassaras I. , Kapetanidis V. , Karakonstantis A. , 2016 . On the spatial distribution of seismicity and the 3D tectonic stress field in western Greece , Phys. Chem. Earth, Parts A/B/C , 95 , 50 – 72 . 10.1016/j.pce.2016.03.012 Google Scholar Crossref Search ADS WorldCat Crossref Kassaras I. et al. . , 2014 . Constraints on the dynamics and spatio-temporal evolution of the 2011 Oichalia seismic swarm (SW Peloponnesus, Greece) , Tectonophysics , 614 , 100 – 127 . 10.1016/j.tecto.2013.12.012 Google Scholar Crossref Search ADS WorldCat Crossref Kaviris G. , Fountoulakis I. , Spingos I. , Millas C. , Papadimitriou P. , Drakatos G. , 2018 . Mantle dynamics beneath Greece from SKS and PKS seismic anisotropy study , Acta Geophys ., 66 , 1341 – 1357 . 10.1007/s11600-018-0225-z Google Scholar Crossref Search ADS WorldCat Crossref Kennett B.L.N. , Engdahl E.R. , Buland R. , 1995 . Constraints on seismic velocities in the earth from travel times , Geophys. J. Int. , 122 , 108 – 124 . 10.1111/j.1365-246X.1995.tb03540.x Google Scholar Crossref Search ADS WorldCat Crossref Kind R. et al. . , 2015 . Thickness of the lithosphere beneath Turkey and surroundings from S-receiver functions , Solid Earth , 6 , 971 – 984 . 10.5194/se-6-971-2015 Google Scholar Crossref Search ADS WorldCat Crossref Kissel C. , Laj C. , 1988 . The Tertiary geodynamical evolution of the Aegean arc: a paleomagnetic reconstruction , Tectonophysics , 146 , 183 – 201 . 10.1016/0040-1951(88)90090-X Google Scholar Crossref Search ADS WorldCat Crossref Klein F.W. , 1989 . HYPOINVERSE, a program for VAX computers to solve for earthquake locations and magnitudes , U.S. Geological Survey Open-File Report: 89–314, doi:10.3133/ofr89314 OpenURL Placeholder Text WorldCat Kokinou E. , Kamberis E. , Vafidis A. , Monopolis D. , Ananiadis G. , Zelilidis A. , 2005 . Deep seismic reflection data from offshore western Greece: a new crustal model for the Ionian Sea , J. Pet. Geol. , 28 , 185 – 202 . 10.1111/j.1747-5457.2005.tb00079.x Google Scholar Crossref Search ADS WorldCat Crossref Kokinou E. , Papadimitriou E. , Karakostas V. , Kamberis E. , Vallianatos F. , 2006 . The Kefalonia Transform Zone (offshore Western Greece) with special emphasis to its prolongation towards the Ionian Abyssal Plain , Mar. Geophys. Res. , 27 , 241 – 252 . 10.1007/s11001-006-9005-2 Google Scholar Crossref Search ADS WorldCat Crossref Konstantinou K.I. , Melis N.S. , 2008 . High-frequency shear-wave propagation across the Hellenic subduction zone , Bull. seism. Soc. Am. , 98 , 797 – 803 . 10.1785/0120060238 Google Scholar Crossref Search ADS WorldCat Crossref Koulakov I. , Kaban M.K. , Tesauro M. , Cloetingh S. , 2009 . P and S velocity anomalies in the upper mantle beneath Europe from tomographic inversion of ISC data , Geophys. J. Int. , 179 , 345 – 366 . 10.1111/j.1365-246X.2009.04279.x Google Scholar Crossref Search ADS WorldCat Crossref Kreemer C. , Holt W.E. , Haines A.J. , 2003 . An integrated global model of present-day plate motions and plate boundary deformation , Geophys. J. Int. , 154 , 8 – 34 . 10.1046/j.1365-246X.2003.01917.x Google Scholar Crossref Search ADS WorldCat Crossref Laigle M. , Hirn A. , Sachpazi M. , Clément C. , 2002 . Seismic coupling and structure of the Hellenic subduction zone in the Ionian Islands region , Earth planet. Sci. Lett. , 200 , 243 – 253 . 10.1016/S0012-821X(02)00654-4 Google Scholar Crossref Search ADS WorldCat Crossref Laske G. , Masters G. , 1996 . Constraints on global phase velocity maps from long‐period polarization data , J. geophys. Res. , 101 , 16059 – 16075 . 10.1029/96JB00526 Google Scholar Crossref Search ADS WorldCat Crossref Laske G. , Masters G. , Ma Z. , Pasyanos M. , 2013 . Update on CRUST1.0—a 1-degree Global Model of Earth's Crust , Geophys. Res. Abstracts , 15 , Abstract EGU2013 – 2658 . OpenURL Placeholder Text WorldCat Legendre C.P. , Meier T. , Lebedev S. , Friederich W. , Viereck-Götte L. , 2012 . A shear wave velocity model of the European upper mantle from automated inversion of seismic shear and surface waveforms , Geophys. J. Int. , 191 , 282 – 304 . 10.1111/j.1365-246X.2012.05613.x Google Scholar Crossref Search ADS WorldCat Crossref Levshin A.L. , Ratnikova L.I. , Berger J. , 1992 . Peculiarities of surface wave propagation across the Central Eurasia , Bull. seism. Soc. Am. , 82 , 2464 – 2493 . OpenURL Placeholder Text WorldCat Li C. , van der Hilst R.D. , Engdahl E.R. , Burdick S. , 2008 . A new global model for 3-D variations of P-wave velocity in the Earth's mantle , Geochem. Geophys. Geosyst. , 9 , Q05018 , doi:10.1029/2007GC001806 . 10.1029/2007GC001806 OpenURL Placeholder Text WorldCat Crossref Li X. et al. . , 2003 . Receiver function study of the Hellenic subduction zone: imaging crustal thickness variations and the oceanic Moho of the descending African lithosphere , Geophys. J. Int. , 155 , 733 – 748 . 10.1046/j.1365-246X.2003.02100.x Google Scholar Crossref Search ADS WorldCat Crossref Lin F.C. , Ritzwoller M.H. , 2011 . Helmholtz surface wave tomography for isotropic and azimuthally anisotropic structure , Geophys. J. Int. , 186 , 1104 – 1120 . 10.1111/j.1365-246X.2011.05070.x Google Scholar Crossref Search ADS WorldCat Crossref Lin F.-C. , Ritzwoller M.H. , Snieder R. , 2009 . Eikonal tomography: surface wave tomography by phase front tracking across a regional broadband seismic array , Geophys. J. Int. , 177 , 1091 – 1110 . 10.1111/j.1365-246X.2009.04105.x Google Scholar Crossref Search ADS WorldCat Crossref Louvari E. , Kiratzi A.A. , Papazachos B.C. , 1999 . The Cephalonia Transform Fault and its extension to western Lefkada Island (Greece) , Tectonophysics , 308 , 223 – 236 . 10.1016/S0040-1951(99)00078-5 Google Scholar Crossref Search ADS WorldCat Crossref Lŭ L. , Wang M. , Laib C.-H. , 2008 . Image denoise by Fourth-Order PDE based on the changes of laplacian , J. Algorithm. Comput. Technol. , 2 , 99 – 110 . 10.1260/174830108784300295 Google Scholar Crossref Search ADS WorldCat Crossref Maggi A. , Priestley K. , 2005 . Surface waveform tomography of the Turkish–Iranian Plateau , Geophys. J. Int. , 160 , 1068 – 1080 . 10.1111/j.1365-246X.2005.02505.x Google Scholar Crossref Search ADS WorldCat Crossref Maupin V. , 2016 . Surface Wave Inversion , Encyclopedia of Earthquake Engineering , Springer-Verlag , doi:10.1007/978-3-642-36197-5_377-1 . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC McClusky S. et al. . , 2000 . Global Positioning System constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus , J. geophys. Res. , 105 , 5695 – 5719 . 10.1029/1999JB900351 Google Scholar Crossref Search ADS WorldCat Crossref Meier T. , Becker D. , Endrun B. , Rische M. , Bohnhoff M. , Stöckhert B. , Harjes H.-P. , 2007 . A model for the Hellenic subduction zone in the area of Crete based on seismological investigations , in The Geodynamic of the Aegean and Anatolia , pp. 183 – 199 ., eds Taymaz T. , Yilmaz Y. , Dilek Y. , Geological Society, London Special Publications . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Meier T. , Dietrich K. , Stöckhert B. , Harjes H.-P. , 2004 . One-dimensional models of shear wave velocity for the eastern Mediterranean obtained from the inversion of Rayleigh wave phase velocities and tectonic implications , Geophys. J. Int. , 156 , 45 – 58 . 10.1111/j.1365-246X.2004.02121.x Google Scholar Crossref Search ADS WorldCat Crossref Meighan H.E. , Brink U. , Pulliam J. , 2013 . Slab tears and intermediate‐depth seismicity , Geophys. Res. Lett. , 40 , 4244 – 4248 . 10.1002/grl.50830 Google Scholar Crossref Search ADS WorldCat Crossref Meijer P.T. , Wortel M.J.R. , 1996 . Temporal variation in the stress field of the Aegean region , Geophys. Res. Lett. , 23 , 439 – 442 . 10.1029/96GL00380 Google Scholar Crossref Search ADS WorldCat Crossref Mindevalli O.Y. , Mitchell B.J. , 1989 . Crustal structure and possible anisotropy in Turkey from seismic surface wave dispersion , Geophys. J. Int. , 98 , 93 – 106 . 10.1111/j.1365-246X.1989.tb05516.x Google Scholar Crossref Search ADS WorldCat Crossref Mitropoulos P. , Tarney J. , 1992 . Significance of mineral composition variations in the Aegean Island Arc , J. Volc. Geotherm. Res. , 51 , 283 – 303 . 10.1016/0377-0273(92)90104-L Google Scholar Crossref Search ADS WorldCat Crossref Mitropoulos P. , Tarney J. , Saunders A.D. , Marsh N.G. , 1987 . Petrogenesis of cenozoic volcanic rocks from the Aegean volcanic arc , J. Volc. Geotherm. Res. , 32 , 177 – 193 . 10.1016/0377-0273(87)90043-6 Google Scholar Crossref Search ADS WorldCat Crossref Olive J.-A. , Pearce F. , Rondenay S. , Behn M.D. , 2014 . Pronounced zonation of seismic anisotropy in the Western Hellenic subduction zone and its geodynamic significance , Earth planet. Sci. Lett. , 391 , 100 – 109 . 10.1016/j.epsl.2014.01.029 Google Scholar Crossref Search ADS WorldCat Crossref Papadimitriou E.E. , Sykes L.R. , 2001 . Evolution of the stress field in the Northern Aegean Sea (Greece) , Geophys. J. Int. , 146 , 747 – 759 . 10.1046/j.0956-540x.2001.01486.x Google Scholar Crossref Search ADS WorldCat Crossref Papadimitriou P. , Kapetanidis V. , Karakonstantis A. , Spingos I. , Pavlou K. , Kaviris G. , Kassaras I. , Voulgaris N. , 2019 . The 25th October, 2018 Zakynthos Earthquake , Bull. Geol. Soc. Greece , Sp. Pub. 7, Ext. Abs. GSG2019-224 . OpenURL Placeholder Text WorldCat Papazachos B.C. , Papazachou C.B. , 1997 . The Earthquakes of Greece , Ziti Editions , Thessaloniki , 304 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Papazachos B.C. , Dimitriadis S.T. , Panagiotopoulos D.G. , Papazachos C.B. , Papadimitriou E.E. , 2005 . Deep structure and active tectonics of the Southern Aegean volcanic arc , Dev. Volcanol. , 7 , 47 – 64 . 10.1016/S1871-644X(05)80032-4 Google Scholar Crossref Search ADS WorldCat Crossref Papazachos C.B. , Kiratzi A.A. , 1996 . A detailed study of the active crustal deformation in the Aegean and surrounding area , Tectonophysics , 253 , 129 – 153 . 10.1016/0040-1951(95)00047-X Google Scholar Crossref Search ADS WorldCat Crossref Papazachos C.B. , Nolet G. , 1997 . P and S deep structure of the Hellenic area obtained by robust nonlinear inversion of travel times , J. geophys. Res. , 102 , 8349 – 8367 . 10.1029/96JB03730 Google Scholar Crossref Search ADS WorldCat Crossref Pasyanos M.E. , Walter W.R. , 2002 . Crust and upper-mantle structure of North Africa, Europe, and the Middle East from inversion of surface waves , Geophys. J. Int. , 149 , 463 – 481 . 10.1046/j.1365-246X.2002.01663.x Google Scholar Crossref Search ADS WorldCat Crossref Paul A. , Karabulut H. , Mutlu A.K. , Salaün G. , 2014 . A comprehensive and densely sampled map of shear-wave azimuthal anisotropy in the Aegean-Anatolia region , Earth planet. Sci. Lett. , 389 , 14 – 22 . 10.1016/j.epsl.2013.12.019 Google Scholar Crossref Search ADS WorldCat Crossref Pearce F.D. , Rondenay S. , Sachpazi M. , Charalampakis M. , Royden L.H. , 2012 . Seismic investigation of the transition from continental to oceanic subduction along the western Hellenic Subduction Zone , J. geophys. Res. , 117 , B07306 . 10.1029/2011JB009023 Google Scholar Crossref Search ADS WorldCat Crossref Pedersen H.A. , Coutant O. , Deschamps A. , Soulage M. , Cotte N. , 2003 . Measuring surface wave velocities beneath small broad-band arrays: test of an improved algorithm and application to the French Alps , Geophys. J. Int. , 154 , 903 – 912 . 10.1046/j.1365-246X.2003.02013.x Google Scholar Crossref Search ADS WorldCat Crossref Pe-Piper G. , Piper D.J.W. , 2002 . The Igneous Rocks of Greece. The Anatomy of an Orogen , Beitrage zur Regionalen Geologie der Erde (Series) , xvi , 573 . OpenURL Placeholder Text WorldCat Piromallo C. , Morelli A. , 2003 . P wave tomography of the mantle under the Alpine-Mediterranean area , 108 , 2065 . OpenURL Placeholder Text WorldCat Pollitz F.F. , Snoke J.A. , 2010 . Rayleigh-wave phase-velocity maps and three-dimensional shear velocity structure of the western US from local non-plane surface wave tomography , Geophys. J. Int. , 180 , 1153 – 1169 . 10.1111/j.1365-246X.2009.04441.x Google Scholar Crossref Search ADS WorldCat Crossref Reilinger R. et al. . , 2006 . GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions , J. geophys. Res. , 111 , B05411 . 10.1029/2005JB004051 Google Scholar Crossref Search ADS WorldCat Crossref Reilinger R. , McClusky S. , Paradissis D. , Ergintav S. , Vernant P. , 2010 . Geodetic constraints on the tectonic evolution of the Aegean region and strain accumulation along the Hellenic subduction zone , Tectonophysics , 488 , 22 – 30 . 10.1016/j.tecto.2009.05.027 Google Scholar Crossref Search ADS WorldCat Crossref Robertson A.H.F. , Dixon J.E. , 1984 . Introduction: aspects of the geological evolution of the eastern Mediterranean , in The Geological Evolution of the Eastern Mediterranean , Vol. 17 , pp. 1 – 74 ., eds Dixon J.E. , Robertson A.H.F. , Geological Society, London Special Publications . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Robertson A.H.F. , 2004 . Development of concepts concerning the genesis and emplacement of Tethyan ophiolites in the Eastern Mediterranean Tethyan region , Lithos , 65 , 1 – 67 . 10.1016/S0024-4937(02)00160-3 Google Scholar Crossref Search ADS WorldCat Crossref Rontogianni S. , Konstantinou K.I. , Melis N.S. , Evangelidis C.P. , 2011 . Slab stress field in the Hellenic subduction zone as inferred from intermediate-depth earthquakes , Earth Planets Space . 63 , 139 – 144 . 10.5047/eps.2010.11.011 Google Scholar Crossref Search ADS WorldCat Crossref Royden L.H. , Papanikolaou D.J. , 2011 . Slab segmentation and late Cenozoic disruption of the Hellenic arc , Geochem. Geophys. Geosyst. , 12 , Q03010 , doi:10.1029/2010GC003280 . 10.1029/2010GC003280 Google Scholar Crossref Search ADS WorldCat Crossref Sachpazi M. et al. . , 2007 . Moho topography under central Greece and its compensation by Pn time-terms for the accurate location of hypocenters: the example of the Gulf of Corinth 1995 Aigion earthquake , Tectonophysics , 440 , 53 – 65 . 10.1016/j.tecto.2007.01.009 Google Scholar Crossref Search ADS WorldCat Crossref Sachpazi M. et al. . , 2016 . Segmented Hellenic slab rollback driving Aegean deformation and seismicity , Geophys. Res. Lett. , 43 , 651 – 658 . 10.1002/2015GL066818 Google Scholar Crossref Search ADS WorldCat Crossref Salaün G. et al. . , 2012 . High-resolution surface wave tomography beneath the Aegean-Anatolia region: Constraints on upper-mantle structure , Geophys. J. Int. , 190 , 406 – 420 . 10.1111/j.1365-246X.2012.05483.x Google Scholar Crossref Search ADS WorldCat Crossref Schmid C. , van der Lee S. , Giardini D. , 2006 . Correlated shear and bulk moduli to 1400 km beneath the Mediterranean region , Phys. Earth planet. Inter. , 159 , 213 – 224 . 10.1016/j.pepi.2006.07.003 Google Scholar Crossref Search ADS WorldCat Crossref Serpen U. , Aksoy N. , Tahir Ö. , 2010 , 2010 present status of geothermal energy in Turkey . Proceedings of the 35th Workshop on Geothermal Reservoir Engineering , 1–3 February 2010, SGP-TR-188, Stanford University, Stanford, CA . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Serpetsidaki A. et al. . , 2014 . New constraints from seismology and geodesy on the Mw = 6.4 2008 Movri (Greece) earthquake: evidence for a growing strike-slip fault system , Geophys. J. Int. , 198 , 1373 – 1386 . 10.1093/gji/ggu212 Google Scholar Crossref Search ADS WorldCat Crossref Shapiro N.M. , Ritzwoller M.H. , 2002 . Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle , Geophys. J. Int. , 151 , 88 – 105 . 10.1046/j.1365-246X.2002.01742.x Google Scholar Crossref Search ADS WorldCat Crossref Shaw B. , Jackson J. , 2010 . Earthquake mechanisms and active tectonics of the Hellenic subduction zone , Geophys. J. Int. , 181 , 966 – 984 . OpenURL Placeholder Text WorldCat Siebert L. , Simkin T. , 2002 . Volcanoes of the World: an Illustrated Catalog of Holocene Volcanoes and their Eruptions , Smithsonian Institution , Global Volcanism Program Digital Information Series, GVP-3 (http://www.volcano.si.edu) . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sodoudi F. , Brüstle A. , Meier T. , Kind R. , Friederich W. , 2015 . Receiver function images of the Hellenic subduction zone and comparison to microseismicity , Solid Earth , 6 , 135 – 151 . 10.5194/se-6-135-2015 Google Scholar Crossref Search ADS WorldCat Crossref Sodoudi F. et al. . , 2006 . Lithospheric structure of the Aegean obtained from P and S receiver functions , J. geophys. Res. , 111 , doi:10.1029/2005JB003932 . 10.1029/2005JB003932 OpenURL Placeholder Text WorldCat Crossref Spakman W. , Nolet G. , 1988 . Imaging algorithms, accuracy and resolution in delay time tomography , in Mathematical Geophysics: A Survey of Recent Developments in Seismology and Geodynamics , pp. 155 – 187 ., eds Vlaar N.J. , Nolet G. , Wortel M.J.R. , Cloetingh S.A.P. , Springer . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Spakman W. , Van Der Lee S. , Van Der Hilst R.D. 1993 . Travel-time tomography of the European–Mediterranean mantle down to 1400 km , Phys. Earth planet. Inter. , 79 , 3 – 74 . 10.1016/0031-9201(93)90142-V Google Scholar Crossref Search ADS WorldCat Crossref Spakman W. , Wortel M.J.R. , Vlaar N.J. , 1988 . The Hellenic subduction zone: a tomographic image and its geodynamic implications , Geophys. Res. Lett. , 15 , 60 – 63 . 10.1029/GL015i001p00060 Google Scholar Crossref Search ADS WorldCat Crossref Suckale J. , Rondenay S. , Sachpazi M. , Charalampakis M. , Hosa A. , Royden L.H. , 2009 . High-resolution seismic imaging of the western Hellenic subduction zone using teleseismic scattered waves , Geophys. J. Int. , 178 , 775 – 791 . 10.1111/j.1365-246X.2009.04170.x Google Scholar Crossref Search ADS WorldCat Crossref Taymaz T. , Jackson J. , Mckenzie D.P. , 1991 . Active tectonics of the north and central Aegean Sea , Geophys. J. Int. , 106 , 433 – 490 . 10.1111/j.1365-246X.1991.tb03906.x Google Scholar Crossref Search ADS WorldCat Crossref Tiberi C. et al. . , 2000 . Crustal and upper mantle structure beneath the Corinth rift (Greece) from a teleseismic tomography study , J. geophys. Res. , 105 , 28 159 – 28 171 . 10.1029/2000JB900216 Google Scholar Crossref Search ADS WorldCat Crossref Tzanis A. , Efstathiou A. , Chailas S. , Stamatakis M. , 2018 . Evidence of recent plutonic magmatism beneath Northeast Peloponnesus (Greece) and its relationship to regional tectonics , Geophys. J. Int. , 212 , 1600 – 1626 . 10.1093/gji/ggx486 Google Scholar Crossref Search ADS WorldCat Crossref van der Meer D.G. , van Hinsbergen D.J.J. , Spakman W. , 2018 . Atlas of the underworld: slab remnants in the mantle, their sinking history, and a new outlook on lower mantle viscosity , Tectonophysics , 723 , 309 – 448 . 10.1016/j.tecto.2017.10.004 Google Scholar Crossref Search ADS WorldCat Crossref van Hinsbergen D. , Kaymakci N. , Spakman W. , Torsvik T.H. , 2010 . Reconciling the geological history of western Turkey with plate circuits and mantle tomography , Earth planet. Sci. Lett. , 297 , 674 – 686 . 10.1016/j.epsl.2010.07.024 Google Scholar Crossref Search ADS WorldCat Crossref Van Hinsbergen D.J.J. , Edwards M.A. , Govers R. , 2009 . Geodynamics of collision and collapse at the Africa –Arabia –Eurasia subduction zone – an introduction , Geol. Soc., Lond., Spec. Publ. , 311 , 1 – 7 . Google Scholar Crossref Search ADS WorldCat van Hinsbergen D.J.J. , Hafkenscheid E. , Spakman W. , Meulenkamp J.E. , Wortel M.J.R. , 2005 . Nappe stacking resulting from subduction of oceanic and continental lithosphere below Greece , Geology , 33 , 325 – 328 . 10.1130/G20878.1 Google Scholar Crossref Search ADS WorldCat Crossref von Huene R. , Reston T. , Kukowski N. , Dehghani G. , Weinrebe W. , 1997 . A subducting seamount beneath the Mediterranean Ridge , Tectonophysics , 271 , 249 – 261 . 10.1016/S0040-1951(96)00241-7 Google Scholar Crossref Search ADS WorldCat Crossref Wessel P. , Smith W.H.F. , 1998 . New, improved version of generic mapping tools released , EOS, Trans. Am. Geophys. Un. , 79 , 579 . 10.1029/98EO00426 Google Scholar Crossref Search ADS WorldCat Crossref Westaway R. , 1991 . Continental extension on sets of parallel faults: observational evidence and theoretical models , in The Geometry of Normal Faults , Vol. 56 , pp. 143 – 169 ., eds Roberts A.M. , Yielding G. , Freeman B. , Geological Society London, Special Publication . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Wortel M.J.R. , Spakman W. , 2000 . Subduction and slab detachment in the Mediterranean–Carpathian region , Science , 290 , 1910 – 1917 . 10.1126/science.290.5498.1910 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Zelt B.C. , Taylor B. , Sachpazi M. , Hirn A. , 2005 . Crustal velocity and Moho structure beneath the Gulf of Corinth, Greece , Geophys. J. Int. , 162 , 257 – 268 . 10.1111/j.1365-246X.2005.02640.x Google Scholar Crossref Search ADS WorldCat Crossref Zhu H. , Bozdağ E. , Tromp J. , 2015 . Seismic structure of the European upper mantle based on adjoint tomography , Geophys. J. Int. , 201 , 18 – 52 . 10.1093/gji/ggu492 Google Scholar Crossref Search ADS WorldCat Crossref Laske G. , Masters G.A. , 1997 . Global Digital Map of Sediment Thickness , EOS Trans. AGU , 78 , F483 . OpenURL Placeholder Text WorldCat APPENDIX: SINGLE-EVENT PHASE-VELOCITY MAPS As described in Section 3, phase velocity variation maps are constructed for each period of interest by inverting the derived phase-delays between station-pairs for each teleseismic event, by applying the Eikonal and Helmholtz equations together with quality restrictions regarding the CC and dispersion measurements. Single-event phase-velocity Eikonal maps provide, in general, a rough image of the lithospheric structure and the margins of the major plates of Africa, Anatolia and the Aegean microplate. However, their quality is limited by the inhomogeneous distribution of stations, along with other site-dependent factors that may bias the results. Consequently, lateral inconsistencies may occur, likely the effect of multipathing related to the Rayleigh-wave propagation and the fact that station-pairs baselines may be largely offset from the great-circle path. These influences are treated by the Helmholtz correction term (see Section 3). To obtain an understanding of the Helmholtz correction, a series of amplitude maps per event/period are generated from the minimum-curvature surface fit, the preliminary correction term and the smoothed correction term. Moreover, the influence of focusing and defocusing of the Rayleigh wavefield that should produce amplitude biases is investigated by mapping the differences between local propagation directions estimated from the apparent phase gradients and the epicentre-station great-circle path (Fig. A1). Figure A1. Open in new tabDownload slide Single-event measurements at a period of T = 50 s for an earthquake that occurred on 30 January 2016 03:25:12 UTC in Kamchatka Peninsula (ray path shown in red on the map at the lower-left corner), (a) phase-delays relative to one of the stations (zero relative traveltime) with valid measurements, (b) propagation direction, (c) differences between propagation direction and great-circle path and (d) apparent velocity map through the Eikonal equation. Contour spacing corresponds to ½ of the respective wavelength while vectors denote the propagation direction (normal to the contours). Figure A1. Open in new tabDownload slide Single-event measurements at a period of T = 50 s for an earthquake that occurred on 30 January 2016 03:25:12 UTC in Kamchatka Peninsula (ray path shown in red on the map at the lower-left corner), (a) phase-delays relative to one of the stations (zero relative traveltime) with valid measurements, (b) propagation direction, (c) differences between propagation direction and great-circle path and (d) apparent velocity map through the Eikonal equation. Contour spacing corresponds to ½ of the respective wavelength while vectors denote the propagation direction (normal to the contours). Deviations from the predicted great-circle path were found to significantly vary per event and period, in a range ±25°. When the Helmholtz amplitude correction term is applied in the single-event maps, the bias from multipathing effects is significantly reduced and the shape of the phase-velocity anomalies becomes more relevant to the subsurface geology (Fig. A2). Figure A2. Open in new tabDownload slide Single-event measurements at a period of T = 50 s for the same event as in Fig. S12. (a, d) amplitude minimum curvature surface fit, (b, e) amplitude correction term before smoothing, (c, f) amplitude correction term after smoothing, (h, i) structural phase-velocity (Helmholtz) map, (a, b, c, h) before and (d, e, f, i) after correction for amplitude-ratio biases. The respective Eikonal map (g; same as Fig. S12d) is displayed for comparison. Figure A2. Open in new tabDownload slide Single-event measurements at a period of T = 50 s for the same event as in Fig. S12. (a, d) amplitude minimum curvature surface fit, (b, e) amplitude correction term before smoothing, (c, f) amplitude correction term after smoothing, (h, i) structural phase-velocity (Helmholtz) map, (a, b, c, h) before and (d, e, f, i) after correction for amplitude-ratio biases. The respective Eikonal map (g; same as Fig. S12d) is displayed for comparison. Interestingly, largest deviations from the predicted great-circle path are concentrated in the azimuths of roughly N150° and N220°, with reversal of the anomaly polarity occurring at N200° (Fig. S14 in Supporting Information). For more examples, the reader is referred to the Supporting Information. © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Deep structure of the Hellenic lithosphere from teleseismic Rayleigh-wave tomography JF - Geophysical Journal International DO - 10.1093/gji/ggz579 DA - 2020-04-01 UR - https://www.deepdyve.com/lp/oxford-university-press/deep-structure-of-the-hellenic-lithosphere-from-teleseismic-rayleigh-1uoNHlCtyk SP - 205 VL - 221 IS - 1 DP - DeepDyve ER -