The Slab Puzzle of the Alpine‐Mediterranean Region: Insights From a New, High‐Resolution, Shear Wave Velocity Model of the Upper Mantle

Mediterranean tectonics since the Lower Cretaceous has been characterized by a multiphase subduction and collision history with temporally and spatially variable, small‐scale plate configurations. A new shear wave velocity model of the Mediterranean upper mantle (MeRE2020), constrained by a very large set of over 200,000 broadband (8–350 s), interstation, Rayleigh wave, phase velocity curves, illuminates the complex structure and fragmentation of the subducting slabs. Phase velocity maps computed using these measurements were inverted for depth‐dependent, shear wave velocities using a stochastic particle‐swarm‐optimization (PSO) algorithm. The resulting three‐dimensional (3‐D) model makes possible an inventory of slab segments across the Mediterranean. Fourteen slab segments of 200–800 km length along‐strike are identified. We distinguish three categories of subducted slabs: attached slabs reaching down to the bottom of the model; shallow slabs of shorter length in downdip direction, terminating shallower than 300 km depth; and detached slab segments. The location of slab segments are consistent with and validated by the intermediate‐depth seismicity, where it is present. The new high‐resolution tomography demonstrates the intricate relationships between slab fragmentation and the evolution of the relatively small and highly curved subduction zones and collisional orogens characteristic of the Mediterranean realm.

The counterclockwise rotation of Africa since the Cretaceous is a key indicator for the larger amount of subducted oceanic lithosphere in the east compared to the west. At present, the eastern and western Mediterranean are separated by the about 300 km wide, continental remnant of the Adriatic microplate that, on its general northwestward path (e.g., Platt & Vissers, 1989), has been colliding with Eurasia to form the Alps, the Dinarides, and northern Hellenides (Figure 1; Carminati & Doglioni, 2012;Handy et al., 2015;Király et al., 2018;Rosenbaum & Lister, 2004;van Hinsbergen et al., 2019).
The Mediterranean (Figure 1) is well suited to study the interaction of tectonic processes like fragmented subduction, back-arc extension, collision, plate deformation and plate creation. It is worth noting that subduction itself has been first proposed by Ampferer (1906) as part of what he called the "Unterströmungstheorie" to explain missing oceanic strata in the Alps. The eastern Mediterranean also presents spectacular examples of the concept of slab rollback (Spakman et al., , 1993, and related surface deformation (see also the review of Jolivet et al., 2013). Intermediate-depth (70-300 km) and deep (300-700 km) seismicity is the classical way to identify the geometry of subducted slabs in the upper mantle. Figure 2 shows the distribution of the intermediate-depth and deep seismic events in the Mediterranean realm for the time period 1990-2019 according to the ISC catalogue (Bondár & Storchak, 2011). However, only the Calabrian slab segment below southern Italy, subducting northwestward beneath the Calabrian Arc, can be identified clearly in the Mediterranean upper mantle using seismicity: An almost continuous Wadati-Benioff zone is present down to the 660 km discontinuity (e.g., Faccenna et al., 2005;Faccenna et al., 2014;Giardini & Velonà, 1991;Selvaggi & Chiarabba, 1995). Only 10 deep earthquakes (∼7) in the westernmost Mediterranean indicate a slab in the mantle transition zone beneath an area extending from the Gibraltar Arc along the Betic Cordillera to the Balearic Islands. The intermediate-depth seismic activity beneath the Alboran Sea, down to a depth of 150 km, is sparse (e.g., Grevemeyer et al., 2015).  Müller et al. (2008) and the main tectonic features after Faccenna et al. (2014). The white arrows show the absolute plate motion with respect to Eurasia after Becker et al. (2015). In the eastern Mediterranean, slabs in the Aegean and Antalya subduction zones can be outlined by intermediate-depth seismicity down to a maximum depth of about 200 km (Bocchini et al., 2018;Caputo et al., 1970;Knapmeyer, 1999;Le Pichon & Angelier, 1979;Papazachos & Comninakis, 1971;Papazachos et al., 2005). A patch of intermediate-depth seismicity in the Vrancea area in Romania (e.g., Hurukawa et al., 2008) indicates a slab segment beneath the southeastern Carpathians. Sparse intermediate-depth seismicity is also found below the Hellenides (e.g., Kilias, 2018) and the northern Apennines (e.g., Chiarabba et al., 2005). Because of the overall low intermediate-depth and deep seismic activity in most of the Mediterranean subduction zones, the outlines of the slabs are difficult to infer using seismicity alone.
Seismic tomographic imaging of the Mediterranean deep structure is thus essential to better understand and constrain the geometry of the slabs, and further the driving forces of the tectonic processes. However, because of the small size of the subduction zones and their spatially highly variable nature this remains a challenging task. Furthermore, the seismological stations are not evenly distributed ( Figure 3; see also Table S1 in the supporting information and references therein), even though for many stations, digital data are now however available for more than 25 years. The tomographic imaging of the Mediterranean slab structure has been pioneered by Spakman et al. (1988Spakman et al. ( , 1993, Spakman (1990), Bijwaard et al. (1998), and Kárason and Van Der Hilst (2000). Wortel and Spakman (2000) related subducting slabs to the retreating arcs and pointed out the importance of slab detachment or horizontal tearing for the terminal stage of subduction. Furthermore, they showed that the slab in the Hellenic Subduction Zone is subducting into the lower mantle, whereas beneath the western Mediterranean and the Pannonian Basin, the slabs are resting on the 660 km discontinuity. Faccenna et al. (2014) reviewed the tectonics and deformation in the Mediterranean and compared the P wave velocity models by Piromallo and Morelli (2003) and by Li et al. (2008) with the S wave velocity models by Zhu et al. (2012) and Auer et al. (2014). They concluded that the large-scale features of these models with wavelength greater than about 300 km agree well but smaller scale features including the about 100 km thick slabs vary considerably due to limited and spatially variable resolution. Recently, van der Meer et al. (2018) identified and described several slabs in the Mediterranean down to the lower mantle based on P wave traveltime tomography (Amaru, 2007), waveform tomography (Schaeffer & Lebedev, 2013), and a global S wave velocity model (Ritsema et al., 2011). They focused, however, on the slab structure below 250 km depth. The fine-scale slab structure in the Mediterranean upper mantle remains uncertain and a matter of heated debates. This concerns, for example, the shallow structure of the Cyprus Slab (e.g., van der Meer et al., 2018), the geometry of the Antalya and Aegean Slabs (e.g., Biryol et al., 2011;Govers & Fichtner, 2016;Salaün et al., 2012), the slab structure and the presence of tears in the transition from the Aegean to the Hellenides and Dinarides (e.g., Handy et al., 2019;Hansen et al., 2019;Šumanovac et al., 2017), the slab configuration and subduction polarity in the Alps (e.g., Hua et al., 2017;Lippitsch et al., 2003;Mitterbauer et al., 2011;Zhao et al., 2016), the shape and presence of slab segments beneath the Apennines and in the transition toward Sicily and the Kabylides in Northern Africa (e.g., Faccenna et al., 2004Faccenna et al., , 2005Faccenna et al., , 2014Fichtner & Villaseñor, 2015;Lucente & Margheriti, 2008), and the slab configuration in the westernmost Mediterranean beneath the Alboran-Betic subduction zone (e.g., Alpert et al., 2013;Bezada et al., 2013;Fichtner & Villaseñor, 2015;Spakman & Wortel, 2004). In summary, a number of slab segments in the Mediterranean upper mantle have been imaged but their configuration, especially down to 300 km depth, remains inconsistently resolved by seismic tomography.
To understand the driving forces of plate motions and deformation in the Mediterranean, high-resolution seismic imaging of the upper mantle is required. Surface waves are well suited to study the three-dimensional shear wave velocity structure in the shallower upper mantle down to about 300 km depth (Biryol et al., 2011;Boschi et al., 2004;Chang et al., 2010;Legendre et al., 2012;Levshin et al., 1989;Marone et al., 2004;Montagner & Tanimoto, 1990;Schivardi & Morelli, 2011;Salaün et al., 2012;Schaeffer & Lebedev, 2013;van der Lee & Nolet, 1997). It remains, however, a challenge to resolve the small slab segments in the Mediterranean.
Here, we use Rayleigh wave tomography to image the entire Mediterranean upper mantle down to about 300 km. An automated procedure to measure phase velocities in a broad period range from 8 to 350 s  is applied to all available data around the Mediterranean for the time period from 1990 to 2015 including, for the first time, data from the Egyptian National Seismological Network (ENSN). Rayleigh wave phase velocity maps for the entire Mediterranean are determined and inverted by a newly elaborated stochastic inversion algorithm, Particle Swarm Optimization (PSO), to obtain a 3-D shear wave velocity model. We focus on the identification of slab segments in the Mediterranean upper mantle and distinguish attached slab segments reaching down to the bottom of the model from slab segments that terminate at shallower depth and slab segments that are detached from the lithosphere above.

Rayleigh Wave Tomography of the Mediterranean Upper Mantle
First, surface wave measurements in the Mediterranean were carried out by Papazachos (1969), Payo and Ruiz de la Parte (1974), Nolet (1977), and Mueller and Sprecher (1978). Calcagnile (1982) presented first rough group velocity maps and noticed the difference between the thinner lithosphere in the western Mediterranean compared to the eastern Mediterranean. Snieder (1988) performed a waveform inversion of surface waves and found indications for the Aegean Slab. Since then group and phase velocity measurements using earthquake data were carried out by Ritzwoller and Levshin (1998), Boschi et al. (2004), Panza et al. (2007, Schivardi and Morelli (2011), Bakırcı et al. (2012), Greve et al. (2014), Palomeras et al. (2014), and Nita et al. (2016). Phase velocity maps determined by Boschi et al. (2004) show indications for slabs beneath the Aegean, the Hellenides, the Dinarides, and beneath Gibraltar.
Furthermore, many local studies were carried out. Bakırcı et al. (2012) suggested a subhorizontal subduction of the north dipping Cyprus Slab. Salaün et al. (2012) imaged with a higher resolution the northward dipping Aegean Slab. Greve et al. (2014) resolved the fine structure of uppermost mantle beneath the Tyrrhenian Sea down to about 160 km depth. Palomeras et al. (2014Palomeras et al. ( , 2017 imaged the Gibraltar Slab between ∼65 and ∼250 km depth. Furthermore, Lu et al. (2018) used ambient noise cross correlations to measure Rayleigh wave group velocities in the period range from 5-150 s and imaged slabs in the Alps, the northern Apennines, and the Dinarides. Surface wave studies using ambient noise focusing on the crust and the uppermost mantle have also been performed in local studies of the Alpine region (Delph et al., 2017;Kästle et al., 2018;Li & Michelini, 2010;Lu et al., 2018;Schippkus et al., 2018;Stehly et al., 2009;Verbeke et al., 2012), to the Pannonian Basin and the Carpathians (Ren et al., 2012), and to the central Mediterranean (Molinari et al., 2015). A surface wave study to resolve the slab configuration in the entire Mediterranean upper mantle has however not been performed yet. Soomro et al. (2016) proposed an algorithm for the automated measurement of interstation phase velocities well suited to analyze large and heterogeneous data sets as those resulting from the inhomogeneous station distribution in and around the Mediterranean (Figure 3). Here, this method is applied to more than 3.5 millions of waveforms recorded by approximately 4,500 broadband seismic stations for ∼3,900 seismic events in the time period from 1990 to 2015. Recordings of the vertical component for events with epicentral distances between 2.7°and 120°and a minimum magnitude varying in dependence on the epicentral distance between 4 M w to 6 M w have been downloaded from the European Integrated Data Archive (http:/www.webdc.eu) and the IRIS data center using Arclink. For references to the seismic networks see Table S1. These recordings were combined with waveform data from 25 broadband seismic stations operated by the Egyptian National Seismological Network (ENSN, https://doi.org/ 10.7914/SN/EY). In Figure 3, epicenters of the seismic events and the station distribution are shown. The seismic stations are color coded according to the number of available waveforms per station. The resulting interstation raypath coverage is given in Figure S1.
For each station pair, approximately located on the same great-circle path, the recorded waveforms are cross correlated and dispersion curves of fundamental modes are calculated from the phase of the cross-correlation functions weighted in the time-frequency plane according to Meier et al. (2004) andSoomro et al. (2016. The cross-correlation function is filtered with a set of frequency-dependent Guassian band-pass filters and windowed in the time domain to minimize the effects from other signals. The automated selection of acceptable phase velocity measurements for each single event is performed in the frequency domain based on a number of fine-tuned quality criteria including the deviation from a path-specific reference model, smoothness criteria, and a criterion for the minimum length of the picked dispersion curve (for details see Soomro et al., 2016). The smooth parts of the single-event dispersion measurements are averaged to obtain reliable and accurate path-average phase velocities that show narrower sensitivity kernels than single-event measurements (de Vos et al., 2013). Between 5 and 2,000 single-event dispersion measurements are averaged per interstation path. The phase velocity measurements are subject to a rigorous quality control . This results in a large number (∼200,000) of broadband measurements spanning the period range from 8 to 350 s with standard deviations mostly lower than 2% and standard errors mostly lower than 0.5% (Figures S1 and S2).

Phase Velocity Maps
Following Deschamps et al. (2008), the obtained interstation phase velocities are used to simultaneously calculate isotropic and azimuthally anisotropic Rayleigh wave phase velocity maps at each individual period. Isotropic Rayleigh wave phase velocity maps at four different periods (15, 30, 60, and 100 s) are depicted in Figure 4. Phase velocity perturbations are plotted with respect to a regional average given on top of each map. Additional phase velocity maps are shown in Figure S3 in the supporting information. Their depth sensitivity is period dependent: Longer periods penetrate deeper and are sensitive to the mantle lithosphere and the asthenosphere (100 s period, ∼150-200 km), whereas short periods are particularly sensitive to the crust (15 s period, ≤15-30 km) (e.g., Laske & Widmer-Schnidrig, 2007;Lebedev et al., 2013).
To estimate the lateral resolution of the phase velocity maps, we perform checkerboard resolution tests. The checkerboard consists of blocks of alternating positive and negative velocity anomalies (±200 m/s) discretized on block sizes of 75, 100, 150, and 200 km, respectively ( Figure S4). Random Guassian noise with a standard deviation of 40 m/s corresponding to 20% of the checkerboard perturbations has been added to the synthetic data. The reconstructed models are shown for 30 and 100 s period in Figures S5 and S6, respectively. They show that the anomalies are overall well retrieved in areas with good path coverage, both in geometry and amplitude, especially in the central Europen-Mediterranean region. Due to the spatially variable raypath density, the maps show however a variable lateral resolution ( Figures S5 and S6). At short periods (≤30 s) small-scale features as small as ∼75 km are well resolved throughout the core of the studied region. At 100 s period, the lateral resolution amounts to about 100-150 km in the northern and central Mediterranean. The retrieval of the anomalies decreases to ∼200 km in areas of poor raypath coverage, especially toward the southern, western, and eastern edges of the maps.
The phase velocity maps ( Figure 4) indicate strong lateral variability of the crust and upper mantle structure in the Mediterranean. At short periods of 15 s (Figure 4a), very low phase velocities are due to sedimentary basins like the Po Basin (PoB: Figure 1), the accretionary wedge of the Calabrian Arc (Cernobori et al., 1996;de Voogd et al., 1992;Finetti et al., 2005;Lenci et al., 2004), or the marine Alboran Basin beneath the Alboran Sea (AS: Figure 1). High phase velocities in the western Mediterranean at 15 s ( Figure 4a) indicate thin crust. At 30 s period (Figure 4b), a narrow and continuous band of low Rayleigh wave phase velocities is observed that extends from the Atlas Mountains, to Sicily, to the Apennines and the Alpine Arc, the Carpathians, Dinarides, and Anatolia. It is related to a continuous belt of deep crustal roots (Laske et al., 2013). The sharp velocity contrast parallel to the coastline of the Italian peninsula reflects an abrupt increase of the crustal thickness from the Tyrrhenian Sea toward the Apennines (crustal models by, e.g., Greve et al., 2014;Molinari et al., 2015).
The phase velocity map at 60 s period ( Figure 4c) is mainly sensitive to lateral variations in shear wave velocity at depth between about 70 and 120 km related to variable properties of the lithosphere-asthenosphere system. Pronounced high-velocity anomalies indicate the East European Craton (EEC), the Western European Continental Mantle Lithosphere (Averbuch & Piromallo, 2012;Legendre et al., 2012), and oceanic lithosphere in the eastern Mediterranean from the Calabrian Arc (Cal: Figure 1) and the Malta Escarpment (MES: Figure 1) in the west to the Eratosthenes Seamount (ESM: Figure 1) in the east. Embedded in low phase velocity anomalies, moderately sized high-velocity anomalies are pointing to the presence of smaller slab segments in the Mediterranean, for example, in the northern Apennines, beneath the Dinarides, or in the central Alps (Kästle et al., 2018). At 100 s (Figure 4d), phase velocities are in general similar to phase velocities at 60 s. In the area of the Calabrian Arc, a narrow high-velocity anomaly is present that may point to the Calabrian Slab. Interestingly, anomalously low velocities beneath the central Apennines and the Dinarides are visible indicating a slab windows in these regions (Faccenna et al., 2014;Greve et al., 2014;Kästle et al., 2018;Spakman & Wortel, 2004;Verbeke et al., 2012).

Inversion for 3-D Isotropic Rayleigh Wave Velocity Model
From phase velocity maps at 100 different periods between 8 and 350 s, local phase velocity dispersion curves are constructed at geographical grid nodes with an average grid spacing of 30 km. Examples for local dispersion curves are shown in Figure S7. In order to relate those local phase velocities quantitatively to shear wave velocities as a function of depth, an efficient stochastic inversion algorithm based on the particle-swarmoptimization technique (Eberhart & Kennedy, 1995;Wilken & Rabbel, 2012) has been implemented and applied. Random perturbations of a background model within specified depth-dependent velocity ranges are tested with respect to the resulting misfit between the measured and synthetic dispersion curves. Nodewise 1-D background models for the inversion are constructed using CRUST1.0 (Laske et al., 2013) for the crust and the isotropic average of PREM (Dziewonski & Anderson, 1981) for the upper mantle. Topography and bathymetry are taken into account. The convergence is speeded up by allowing for a random local search at each position of the particles in the model space. Flexible depth-dependent parameterization ( Figure S8a and S8b) and regularization strategies are applied to velocity perturbations as well as to the depths of discontinuities. The shear wave velocity in the crust is described by four to five layers with variable layer thicknesses. Furthermore, the perturbation of the Moho depth is a parameter of the inversion. In the upper mantle, a cubic perturbation in the shear wave velocity is considered from the Moho down to the 410 km discontinuity. A quadratic perturbation is applied to the mantle transition zone followed by a linear perturbation below the 660 km discontinuity and downward ( Figure S8). Theoretical phase velocities for elastic, isotropic, 1-D model are calculated using the Thomson-Haskell matrix formalism (Knopoff, 1964;Schwab & Knopoff, 1972). Around 8,000 forward calculations are sufficient for the applied parametrization to obtain a good fit and to explore the model space.
Two examples of the 1-D inversion are shown in Figure 5 for the eastern Alps ( Figure 5, left) and Anatolia ( Figure 5, right), respectively. The location of these two model points are given in Figure S7. The blue colors in Figure 5 indicate acceptable models with root-mean-square (rms) misfits lower than an empirically chosen uncertainty threshold. It is calculated by enlarging the misfit of the global minimum by 0.5, which equals half of the rms standard deviation. Thus, we allow for significant increase of the misfit but the main properties of the dispersion curve are still explained by the acceptable models. The width of the bundle of acceptable isotropic Vs models provides information on the depth-dependent model uncertainties. Estimates of the depth-dependent uncertainties are determined by dividing the width of acceptable models at a given depth by 2. The gray colors indicate the range of tested models. The thick dashed line gives an optimal centroid model calculated from the accepted values of the inversion parameters. An example for the determination of the centroid value is given for the velocity at uppermost mantle depth in the eastern Alps in Figure S8c. The black dots represent the misfit for all considered models as a function of the tested parameter values. The blue line outlines the envelope of the misfit. The blue circle indicates the model with the lowest rms misfit, whereas the red triangle indicates the centroid value of the misfit envelope in the range of acceptable models.
Beneath the eastern Alps ( Figure 5, left), a pronounced high shear wave velocity anomaly is found down to a depth of approximately 200 km likely related to subducted mantle lithosphere (Kästle et al., 2018;Lippitsch et al., 2003;Mitterbauer et al., 2011), whereas anomalously thin lithosphere is detected underneath eastern Anatolia ( Figure 5, right) overlying a shallow asthenosphere as indicated by a strong decrease in the shear wave velocities below about 70 km depth. The 3-D shear wave velocity model (MeRE2020) beneath the Mediterranean region down to 300 km depth is constructed by assembling the centroid 1-D Vs profiles resulting from the phase velocity depth inversions at all grid nodes. In Figures 6 and 7 horizontal map views of MeRE2020 are shown for different depths (75, 100, 200, and 250 km). Additional horizontal slices at 80, 125, 230, and 300 km depth are given in Figure S11 in the supporting information. Figures 8 and 9 show vertical cross sections of the shear wave velocity perturbations along eight selected profiles. The velocity perturbations are plotted with respect to the 1-D depthdependent average velocity model shown in Figure 8. Depth-dependent uncertainty estimates of the shear wave velocity are shown in the vertical cross sections in Figure S12 along three selected profiles: The model is well constrained and characterized by mainly low uncertainties between about 80 and 200 km (0.1-0.15 km/s). The uncertainties increase toward shallower depths due to trade-offs with the crustal structure as well as toward greater depth because of the decreasing sensitivity of the Rayleigh waves. Some lateral variability in the model uncertainties is also present. For example, velocities in the southwesternmost Mediterranean are less well determined as evident from profile I-I' in Figure S12.

10.1029/2020GC008993
Geochemistry, Geophysics, Geosystems the Aegean, Anatolia, and the Middle East. In the same depth range, the model is also characterized by the presence of localized, well-resolved, small-scale high-velocity anomalies in subduction zones and beneath orogens.

Discussion of the Slab Segments in the Mediterranean Upper Mantle
In this section we discuss the spatial distribution of the different subducted slab segments in the Mediterranean upper mantle down to 300 km depth as revealed from MeRE2020 in comparison to previous studies. In a global study, van der Meer et al. (2018)   The upper 300 km of the mantle are however not the focus of the study by van der Meer et al. (2018). Therefore, slabs with maximum depth of less than 300 km, shallow slab gaps and slab tears are not identified. Using the horizontal cross sections in Figures 6, 7, and S9, and specifically the vertical cross sections in Figures 8 and 9, we revisit the presence of slab segments in the Mediterranean upper mantle down to about 300 km depth. A slab segment is defined if the corresponding high-velocity anomaly meets the following conditions: (i) its velocity is 2% greater than the depth-dependent average velocity model ( Figure 8); (ii) it is strongly dipping; (iii) it is located beneath or near a subduction zone, (iv) it is elongated laterally following, at least partly, the strike of a subduction zone; (v) it has a width between about 80-250 km, (vi) if the lateral resolution indicated by the checkerboard tests ( Figures S4, S5, and S6) is better than 150 km in the region of the anomaly; and (vii) if present, intermediate-depth seismicity at the location of the high-velocity anomaly is used to identify the top of the slab segment.
In order to account for the complexity of the slab segments, we distinguish between short slabs if the base of the slab segment is found at depth lower than 200 km, attached slabs that reach the bottom of the model, and detached slabs that do not show a connection to the lithosphere in the foreland. Following the notation suggested by van der Meer et al. (2018), we discuss our findings starting from the east referring to horizontal and  Figure 10 the location as well as the dip direction of the attached and shallow slab segments at 100 km depth in the Mediterranean upper mantle is shown on a topographic map. The orange and blue colors indicate African/Adriatic and Eurasian slab segments, respectively. The approximate maximum depth of the shallow slabs as well as the top of the detached slabs are also given on the map in Figure 10. In the following, we discuss the individual slab segments in greater detail.

Eastern Mediterranean 5.1.1. Cyprus Slab (Cy)
At least down to approximately 150 km depth, MeRE2020 does not show evidence for a slab segment east of Cyprus (Figures 6 and 7). Rather low velocities down to about 150 km that are also characterizing continental lithosphere of the area of the Levant Basin southeast of Cyprus. According to van der Meer et al. (2018), it remains an open question if the Cyprus Slab, that is reaching down to 1,000 km depth in their model, is attached or detached. The low velocities in MeRE2020 down to about 150 km support the suggestion of a slab break-off to the east of Cyprus (Biryol et al., 2011;Faccenna et al., 2006;Portner et al., 2018;Wong A Ton & Wortel, 1997). At 200 km depth, a weak high-velocity anomaly can be found in MeRE2020 east of Cyprus (Figure 7) that may indicate the top of the detached Cyprus slab. At depths larger than 250 km depth, the Cyprus slab may also contribute to the imaged widespread high-velocity anomaly beneath central Anatolia (Figure 7). The lateral resolution of MeRE2020 is however limited to about 200 km in this area because of a lack of stations. In Figure 10, no slab segment is indicated east of Cyprus.

Antalya Slab (Ant)
To the west of Cyprus and beneath the Bay of Antalya, where oceanic lithosphere is entering the subduction zone (Granot, 2016), we image a northeastward-dipping, continuous high-velocity anomaly (profile H-H' in Figure 9), which is well delineated by earthquake foci down to about 130 km depth (Figure 2; Bocchini et al., 2018). We interpret this anomaly as the Antalya Slab (Ant) subducting in the area of western Cyprus to the city Antalya at the southern Anatolian coast. The Antalya Slab is separated by a vertical tear from the Aegean Slab (Biryol et al., 2011;De Boorder et al., 1998;Govers & Fichtner, 2016;Legendre et al., 2012;. The geometry of the tear is visible in the horizontal cross section at depth 250 km depth but is more clearly defined by intermediate-depth seismicity between the Antalya and the Aegean Slabs (Figure 2; e.g., Bocchini et al., 2018;Brüstle et al., 2012). At 250 km ( Figure 7) the Antalya Slab is found beneath southwest Anatolia. This is in accordance with tomographic studies using P wave traveltimes (Biryol et al., 2011;Portner et al., 2018). Note Interestingly, an almost north-south striking high-velocity anomaly is found beneath central Anatolia at depth greater than 200 km (Figure 7). We refer to it as Central Anatolian Anomaly (CAA). In the vertical cross section H-H' in Figure 9, it becomes evident that at shallow depth it is intersected by low velocities due to the shallow asthenosphere beneath Anatolia. This leads to the southward dip of the high-velocity anomaly in profile H-H' that is oriented almost parallel to the strike of this anomaly (Figure 9). It remains an open question if this anomaly is related to Cretaceous subduction in the area (van Hinsbergen et al., 2016).

Aegean Slab (Aeg)
The Aegean Slab is generally the most prominent high-velocity anomaly in tomographic models of the Mediterranean that is interpreted as African oceanic lithosphere subducting beneath the Hellenic Arc (Bijwaard et al., 1998;Biryol et al., 2011;Chang et al., 2010;De Boorder et al., 1998;De Jonge et al., 1994;Hansen et al., 2019;Piromallo & Morelli, 1997Portner et al., 2018;Spakman et al., 1988;van der Meer et al., 2018;Zhu et al., 2012). MeRE2020 shows a continuous high-velocity anomaly down to approximately 300 km depth (Figures 6 and 7, and profile G-G' in Figure 9). The anomaly gets weaker with depth likely due to lower resolution at larger depth. At shallower depths (≤125 km, Figures 6 and S11), the imaged high-velocity anomaly shows a concave curvature following the shape of the Hellenic Arc. Note that maximum 400-500 km of the African oceanic lithosphere have been subducted south of Crete (Meier et al., 2004;Meulenkamp et al., 1988;Ten Veen & Kleinspehn, 2002;Ten Veen & Meijer, 1998;. Taking the dip angle of the subduction into account, lithosphere subducted south of Crete would reach only down to about 200 km depth. The high-velocity anomaly at depth larger than 200 km is thus supporting the single-slab model suggesting that the long Aegean Slab reaching down to the lower mantle is composed of oceanic lithosphere of several oceanic domains as well as of continental mantle lithosphere of intervening microcontinents (Faccenna et al., 2003;Jolivet & Brun, 2010;Meier et al., 2004;van der Meer et al., 2018;van Hinsbergen et al., 2005a).
Toward the northwest, high velocities indicating the Aegean Slab are decreasing in amplitude but are continuing toward the Hellenides (horizontal cross sections at 150 and 200 km depth in Figures 7 and 8). The Aegean Slab is further characterized by strong intermediate-depth seismicity down to about 180 km depth but not deeper (Figure 2). The intermediate-depth seismicity is also decreasing to the north beneath the Hellenides. These changes toward the Hellenides are likely due to the transition between oceanic and continental subduction (Bocchini et al., 2018). MeRE2020 (horizontal cross section at 200 km depth in Figure 7) is supporting the suggestion that the Kefalonia Transform Fault (KTF: Figure 1) is located above the slab edge of the Aegean Slab (Bocchini et al., 2018). The resolution of MeRE2020 is however not sufficient to resolve if a vertical slab tear is present between the Aegean Slab to the south and the slab beneath the Hellenides to the north as suggested by Govers and Wortel (2005), Suckale et al. (2009), andHansen et al. (2019) but opposed by Halpaap et al. (2018). The presence of a vertical tear is however likely because of considerably larger amount of subducted lithosphere south of the Kefalonia Transform Fault in the last about 15 Ma (Bocchini et al., 2018).
To the south of the Hellenic Subduction Zone, old Triassic oceanic lithosphere is present in the eastern Mediterranean (Granot, 2016). MeRE2020 depicts the distribution of oceanic lithosphere in the eastern Mediterranean (EMOML: profiles F-F', G-G', and H-H' in Figure 9) as high shear wave velocity anomaly from the Malta Escarpment in the west toward the Eratosthenes Seamount along the western margin of the Levant Basin in the east with however a variable thickness throughout the region (horizontal cross sections at 100-200 km depth in Figures 6 and 7). The maximum thickness of the oceanic lithosphere of about 200 km is found south of the Hellenic Subduction Zone (vertical cross section F-F' in Figure 9). The thickness is however reduced to about ∼150 km. To the east, a clear abrupt change in the shear wave velocities is imaged along the western margin of the Levant Basin, which reflects the transition from oceanic to stretched continental lithosphere south of Cyprus (vertical cross section F-F' in Figure 9). The presence of thick oceanic lithosphere in the eastern Mediterranean has been noticed previously, for example, by Boschi et al. (2009) and Legendre et al. (2012).

Southern Carpathian Slab (CpS)
Beneath the southeastern Carpathians, MeRE2020 shows a nearly vertical, high-velocity anomaly in the area of the Vrancea zone down to about 350 km depth (Figures 6 and 7, and profile G-G' in Figure 9). Interestingly, this high-velocity anomaly correlates with the distribution of intermediate-depth seismicity beneath the southeastern Carpathians (Figure 2 and profile K-K' in Figure S12). In accordance with previous studies we relate this anomaly to subducting oceanic lithosphere that is still attached to the southeastern foreland Koulakov et al., 2010;Martin et al., 2006;Weidle et al., 2005;Wortel & Spakman, 2000). MeRE2020 reveals that the high-velocity anomaly may extend to the west beyond the area of intermediate-depth seismicity (horizontal cross section at 250 km in Figure 7 and at 230 and 300 km depth in Figure S11). Another high-velocity anomaly is observed beneath the central Pannonian Basin at depth of about 230 km and downward. This anomaly could represent a subducted fragment of the Alpine-Tethyan margin beneath the generally low velocities in the Pannonian Basin uppermost mantle. Its relation to the high-velocity anomaly in the mantle transition zone beneath the Pannonian Basin remains an open question (Bijwaard et al., 1998;Dando et al., 2011;Ren et al., 2012;Spakman, 1991;Spakman et al., 1993;Wortel & Spakman, 2000;van der Meer et al., 2018;Zhu et al., 2012).

Central Mediterranean 5.2.1. Hellenides Slab (He)
Down to depths larger than 300 km, a clear east dipping high-velocity anomaly beneath the Hellenides is imaged (He: Figures 6 and 7). It is also well presented along the vertical cross section E-E' in Figure 9. We interpret this high-velocity anomaly as Adriatic continental lithosphere subducting eastward beneath the Hellenides. van der Meer et al. (2018) do not mention slabs beneath the Hellenides and the Dinarides explicitly. However, a slab beneath the Hellenides has been suggested by Handy et al. (2019) based on the tomographic studies by Piromallo and Morelli (2003), Spakman and Wortel (2004), and Amaru (2007). They favor the presence of a continuous slab beneath the Hellenides. Using scattered waves, Pearce et al. (2012) image continental lithosphere down to about 220 km depth. In contrast to a long continuous slab beneath the Hellenides, a tear propagating horizontally southward beneath the Hellenides has been suggested by Wortel and Spakman (2000). Also Hansen et al. (2019) propose a horizontal tear beneath the downdip limit of intermediate-depth seismicity at about 180 km depth (Figure 2). No evidence for this horizontal tear is found in MeRE2020 down to a depth of about 300 km.

Dinarides Slab (Di)
Along the entire Dinarides MeRE2020 shows a narrow, shallow, northeast dipping high-velocity anomaly however with variable depth extent. In the northern and the central Dinarides, this anomaly extends down to approximately 150 km depth (Figures 6 and 10 and profile B-B' in Figure 8), whereas it reaches depths of 300 km in the southern Dinarides. In the northern and central Dinarides, low velocities are found beneath about 150 km depth. A northeast dipping high-velocity anomaly beneath the southern Dinarides has been previously imaged (Bijwaard & Spakman, 2000;Koulakov et al., 2009;Schmid et al., 2008). It has been attributed to underthrusting of Adriatic lithosphere beneath the Dinarides (Handy et al., 2015(Handy et al., , 2019. In the northern and central parts of the Dinarides Bijwaard and Spakman (2000) reported, however, a large low-velocity anomaly that extends laterally to the eastern Alps. Handy et al. (2015) attributed this anomaly to a slab gap beneath the northern and central Dindarides due to a break-off 10.1029/2020GC008993 Geochemistry, Geophysics, Geosystems of the Adriatic Slab at the onset of continental collision (∼30 Ma ago) and invoked the upwelling of asthenospheric material to explain this low-velocity anomaly. This is however contradicting the results by Šumanovac and Dudjak (2016) and Šumanovac et al. (2017), who imaged a northeast steeply dipping high-velocity anomaly throughout the entire Dinarides with a variable depth extent from (∼450 km) in the south to (∼250 km) in the north with a shallow slab down to only about 100 km in-between. In the northern and central Dinarides our model (see profile B-B' in Figure 8) suggests the presence of both a shallow slab down to about 150 km depth (Di) and a slab gap below (DiG) in accordance with Kästle et al. (2020). In the southern Dinarides, the slab extends to larger depth (horizontal cross sections at 250 km depth in Figure 7and at 230 and 300 km depth, Figure S11). The transition to the slab beneath the Hellenides is however not well resolved.
Furthermore, at approximately 250 km depth, another high-velocity anomaly is depicted beneath the Adriatic Sea (horizontal cross section at 300 km, Figure S11). This high-velocity anomaly has been imaged before by Giacomuzzi et al. (2011). We interpret this anomaly as a detached part of the Adriatic lithosphere. If a sinking rate of about 1 cm/a is adopted (e.g., Handy et al., 2015) and the counterclockwise rotation of Adria since the slab break-off in the Dinarides is taken into account, a current location of the detached part of the Adriatic slab beneath the Adriatic Sea is likely. In Figure 8, this high-velocity anomaly is well presented along the vertical cross sections B-B' and D-D' where it is marked by DiD?.

Alpine Slabs (EA, CA, and WA)
The Alps represent a rather small, highly curved orogen formed by collision of Adria and Eurasia along the current northern margin of Adria. The presence of slab segments beneath the Alps has been suggested by a number of tomographic studies (Bijwaard & Spakman, 2000;Dando et al., 2011;Hua et al., 2017;Kästle et al., 2018;Koulakov et al., 2009;Lippitsch et al., 2003;Mitterbauer et al., 2011;Spakman et al., 1993;Zhao et al., 2016). A high-velocity anomaly in MeRE2020 beneath the central Alps indicates southward subduction of Eurasia (CA: Figures 6 and 7 and profile C-C' in Figure 8), in agreement with previous studies (Dando et al., 2011;Hua et al., 2017;Kästle et al., 2018;Lippitsch et al., 2003;. The south dipping high-velocity anomaly is partly found beneath the northern foreland. This is expected because of the thickness of the slab of about 100 km. Our results also show that the dip of the slab in the central Alps can be imaged by surface wave tomography, although the detailed shape of the south dipping anomaly remains to be resolved. Beneath the western Alps, MeRE2020 shows high velocities down to 100 km (vertical cross section B-B' in Figure 8) interpreted as Eurasian mantle lithosphere (WA) in accordance with Kästle et al. (2018). At larger depths than 100 km, a low-velocity anomaly is imaged (horizontal cross sections at 125 km depth in Figure  S11 and 200 km depth in Figure 7, profile B-B' in Figure 8). This supports the hypothesis of a slab break-off as suggested by Lippitsch et al. (2003), and Kästle et al. (2018) in contrast to models of continuous, eastward dipping subduction beneath the western Alps (Koulakov et al., 2009;Zhao et al., 2016), Hua et al. (2017), andLyu et al. (2017). In the horizontal cross section at 150 km depth, the proximity to the Northern Apenninic Slab (NA) becomes evident. The E-W oriented vertical cross section B-B' in Figure 8 running nearly parallel to the Northern Apenninic deformation front shows the presence of a high-velocity anomaly beneath the southern Po Basin that is however not related to an east dipping Eurasian Slab beneath the Po Basin but to a south dipping slab beneath the northern Apennines.
In the eastern Alps, no continuous deep-reaching, high velocities are found in our model (profile D-D' in Figure 8). This supports the idea of slab break-off in the eastern Alps (Handy et al., 2015;Kissling et al., 2006;Lippitsch et al., 2003;Schmid et al., 2004). There is an ongoing debate on the presence and geometry of shallow slabs in the eastern Alps. The classical view of south dipping Eurasia (Lammerer, 1988;Lüschen et al., 2006) has been opposed by Lippitsch et al. (2003), Schmid et al. (2004), Kissling et al. (2006), andHandy et al. (2015) and a north dipping Adriatic Slab has been proposed to explain the high-velocity anomaly down to about 250 km depth beneath the northeastern foreland of the Alps. Mitterbauer et al. (2011) favor a Eurasian origin of the north dipping high-velocity anomaly in their model. In MeRE2020, the almost continuous extension of a high-velocity anomaly from the Dinarides toward the eastern Alps (horizontal cross section at 100 km, Figure 6 and vertical cross section D-D' in Figure 8) hints at the presence of a shallow Adriatic slab beneath the eastern Alps down to about 150 km depth. In addition, another high-velocity anomaly is found beneath the northern Alps and the Molasse Basin (vertical cross section D-D' in Figure 8). It is interpreted as southward subducting Eurasian mantle lithosphere. Thus, in the 10.1029/2020GC008993 Geochemistry, Geophysics, Geosystems eastern Alps subducting mantle lithosphere of both Eurasian and Adriatic origin may be present. However, a clear understanding of the subducted slabs beneath the eastern Alps remains to be further investigated within the scope of the AlpArray project (Hetényi et al., 2018).

Northern Apenninic Slab (NA)
In the horizontal cross sections at 100 to 200 km depth, there is a pronounced high-velocity anomaly beneath the northern Apennines and the southern Po Basin (Figures 6 and 7). This anomaly extends down to at least 250 km depth and is well imaged on the vertical cross sections B-B' and C-C' in Figure 8. We relate this anomaly to Adriatic mantle lithosphere subducting southward beneath the northern Apennines (NA). The presence of this slab has been indicated by a number of previous tomographic studies (Benoit et al., 2011;Giacomuzzi et al., 2011;Kästle et al., 2018;Koulakov et al., 2009;Lucente et al., 1999;Rosenbaum et al., 2008;Spakman & Wortel, 2004;van der Meer et al., 2018;Zhao et al., 2016). Its lateral and vertical extent is however debated. Spakman and Wortel (2004) suggests that the slab extends down to 300-400 km depth, whereas Giacomuzzi et al. (2011) imaged a continuous high-velocity anomaly only down to a depth of approximately 140 km. Our model supports the hypotheses of a longer Northern Apenninic Slab.
Although the intermediate-depth seismic activity is low beneath the northern Apennines (Chiarabba et al., 2005), the top of the Northern Apenninic Slab can be identified by a few intermediate-depth events. At about 70 to 100 km depth, the top of the slab is located beneath the northern Apennines ( Figure 2) as expected for a nearly vertically, south dipping slab. Also, in our horizontal and vertical cross sections (Figures 6-8), a high-velocity anomaly is oriented parallel to the northern Apennines. It extends to the north beneath the Po Basin by about 100 km because of the thickness of the nearly vertically dipping slab (e.g., vertical cross section B-B', Figure 8). Because the lateral resolution in the area is about 75 km, Rayleigh wave has a stronger sensitivity to higher velocities at the center of the slab.
Beneath the central Apennines, MeRE2020 indicates a low-velocity anomaly from 70 km down to about ∼250 km depth supporting the presence of the Central Apenninic Slab Gap (CAG: Figure 6, profile E-E' in Figure 9) in accordance with the results of the majority of regional studies (Amato et al., 1993;Lippitsch et al., 2003;Rosenbaum et al., 2008;Spakman & Wortel, 2004;Zhao et al., 2016). At 250 km depth (Figure 7), a high-velocity anomaly is detected below the central Apennines that can be related to the edge of the northwestward dipping Calabrian Slab (Lucente & Margheriti, 2008;Lucente & Speranza, 2001;Wortel & Spakman, 2000). The vertical cross section E-E' in Figure 9 shows the asymmetry of the western and eastern Adriatic margins: A slab gap is found beneath the central Apennines along the western margin whereas beneath the Hellenides a deep reaching slab is imaged.

Calabrian Slab (Cl)
MeRE2020 shows a continuous, steeply northwestward dipping, high-velocity anomaly beneath the Calabrian Arc reaching down to depths larger than 300 km beneath the southeastern Tyrrhenian Sea (Cl: horizontal cross sections in Figures 6 and 7 and vertical cross section F-F' in Figure 9). It points to the Calabrian Slab representing subducting oceanic lithosphere of the Ionian Sea. The slab geometry is also well defined by the distribution of intermediate-depth and deep seismicity down to 550 km (Figure 2). A high-velocity anomaly correlates well with the Wadati-Benioff zone below the Calabrian Arc (profile J-J' in Figure S12). Intermediate-depth seismicity occurs within slabs close to the upper plate interface. Therefore, the seismicity is located at the upper margin of the dipping anomaly not at the center where the largest amplitudes of the anomaly are found. At depth of 200 km, detached parts of the Calabrian Slab extend toward the northwest beneath the southern Apennines as well as to the west beneath Sicily (horizontal cross section in Figure 7) in accordance with tomographic studies by Bijwaard and Spakman (2000), Piromallo and Morelli (2003), and Chiarabba et al. (2005). The slab gap beneath the central Apennines represents thus a slab tear of the Calabrian Slab. According to our model, the slab is attached to the lithosphere in the Ionian Sea in the area between Sicily and Calabria.

Western Mediterranean 5.3.1. Kabylides Slab (Kb)
Along the African coast in the area of the Kabylides, MeRE2020 depicts a high-velocity anomaly down to depth of 250 km (Figures 7 and S9). The high-velocity anomaly beneath the Kabylides is interpreted as subducting African lithosphere (Kb) that caused rollback and back-arc extension since ∼45-30 Ma (Faccenna 10.1029/2020GC008993 Geochemistry, Geophysics, Geosystems Handy et al., 2015;Spakman & Wortel, 2004;van der Meer et al., 2018;van Hinsbergen et al., 2014). It is attached to continental Africa beneath the Kabylides but a part that is reaching toward the region of southern Sardinia seems to be detached (vertical cross section C-C' in Figure 8). This detached part is well separated from the Calabrian Slab in MeRE2020 down to about 200 km depth (Figure 7). At larger depth it merges with a high-velocity anomaly that may be related to the Calabrian Slab (horizontal cross section at 250 km, Figure 7). A detached part is also found to the west of the Kabylides along the northern African coast (vertical cross section I-I' in Figure S12). While the location of the Kabylides Slab (Kb) is better defined by MeRE2020 compared to previous studies Spakman & Wortel, 2004), the geometry of this slab remains poorly resolved in the upper mantle mainly because of the low station density in the area.

Gibraltar and Betic Slabs (Gib and Btc)
MeRE2020 reveals an arcuate, eastward dipping, high-velocity anomaly beneath the Gibraltar Arc that extends down to depth of about 200 km only (Figures 6 and 7, vertical cross section A-A' in Figure 8). Toward the northeast, another weak high-velocity anomaly is found at depth larger than about 230 km along the southeastern Iberian coast beneath the Valencia Trough (Figures 6, 7, and S9). To distinguish between these two slab segments, we refer to the first as Gibraltar Slab segment (Gib) and to the latter as Betic Slab segment (Btc). The intermediate-depth seismicity is confined to the Gibraltar Slab ( Figure 2). It occurred close to the upper plate interface of the Gibraltar Slab as evident from the vertical cross section I-I' in Figure S12. The Gibraltar Slab seems to subduct not deeper than about 200 km where a slab gap is found in accordance with Heit et al. (2017). The intermediate-depth seismicity seems however to be located within the Gibraltar Slab and not at the tear as suggested by Heit et al. (2017). We refer to the Gibraltar Slab segment as shallow slab.
A tear above the Betic slab segment (BtG) shown in the vertical cross section A-A' in Figure 8 has been suggested before by Spakman and Wortel (2004), Garcia-Castellanos andVillaseñor (2011), andde Lis Mancilla et al. (2015). Below the Betic Slab segment, a few isolated deep earthquakes occurred in the mantle transition zone at depth of ∼500-600 km (Figure 2) beneath the southwestern coast of Iberia. This might indicate that the imaged detached slab beneath the Betic penetrates almost vertically down to the mantle transition zone. The deeper portion of the slab in the mantle transition zone, outlined by a few seismic events in Figure 2, has been imaged by body wave tomography (Amaru, 2007;Spakman & Wortel, 2004;van der Meer et al., 2010). Deep events occurred also beneath the Gibraltar area including the 29 March 1954 event (640 km depth, 7.8 M w ). Likely, the Betic Slab segment extends toward the southwest at larger depth beneath 300 km depth below the bottom of our model as indicated by P wave tomography Spakman & Wortel, 2004;van der Meer et al., 2018). Moreover, the Gibraltar Slab is separated from the Kabylides Slab by a wide gap of low velocities down to depths of about 200 km down (horizontal cross sections in Figures 6 and 7 and vertical cross section I-I' in Figure S12). This slab gap might, at least in the east, be underlain by a detached part of subducted African lithosphere.
The discussion of the horizontal and vertical cross sections through MeRE2020 (Figures 6-9) shows that small slab segments as present in the Mediterranean upper mantle can be resolved by a shear wave velocity model obtained from surface wave tomography down to about 300 km depth. They have so far been imaged only partly by local studies or beneath the entire Mediterranean area by P wave tomography however with limited vertical resolution in this depth range. As a result, an inventory of slab segments in the entire Mediterranean is provided for the upper mantle down to about 300 km depth. In Figure 10, we summarize our findings. The strong slab fragmentation is related to the highly curved subduction zones, orogens, and small-scale back-arc basins in the Mediterranean. The slab anomalies in our tomographic model can be used for defining 3-D slab surfaces for geodynamic numerical simulations. It remains however a challenge for further tomographic studies to resolve finer details of the three-dimensional geometry as well as the internal properties of the slab segments.

Conclusions
New, high-resolution Rayleigh wave tomography of the Mediterranean upper mantle down to about 300 km depth reveals the fine-scale structure and strong fragmentation of the subducting slabs beneath the region. An automated technique for the determination of interstation fundamental mode phase velocities was applied to a large, heterogeneous set of data from all publicly available stations around the Mediterranean 10.1029/2020GC008993 Geochemistry, Geophysics, Geosystems in the time period from 1990 to 2015 and to data from the ENSN. A very large data set of over 200,000 interstation dispersion curves was obtained and inverted for a set of phase velocity maps, spanning a broad period range (8-350 s). The maps were then inverted, point by point, for a 3-D shear wave velocity model using a stochastic, particle-swarm-optimization algorithm. The resulting tomographic model (MeRE2020) shows the highly fragmented nature of the subducting slabs in the Mediterranean upper mantle. The location of the slab segments is revealed by high-velocity anomalies in the model, which we interpret together with the spatial distribution of intermediate-depth seismicity. Slab fragmentation and horizontal tearing result in a large number of slab segments that vary in their lateral length between ∼200 and ∼800 km. In total, 14 slab segments are identified. We distinguish attached slab segments reaching down to the bottom of the model from shallow slabs that terminate at shallower depths and from detached ones that are not connected to the lithosphere in the foreland. The lower and upper boundaries of the shallow and detached slabs, respectively, are formed by slab tearing.
There is no evidence for a slab segment to the east of Cyprus within the top 200 km of the mantle, whereas to the northwest of Cyprus a continuous slab segment (Ant) is subducting northeastward at least down to the bottom of the model. Moreover, a high-velocity anomaly elongated north-south and dipping southward is detected beneath Central Anatolia (CAA). In the southeastern Aegean, a northwest dipping slab segment is imaged in the area of Rhodes (Aeg). In the southwestern Aegean, the Aegean Slab is dipping to the northeast. An attached slab is imaged beneath the Hellenides (He), reaching down to at least 300 km depth, whereas beneath the central and northern Dinarides a short slab (Di) is found down to about 150-200 km depth, above a slab tear. We suggest that the top of a deeper, detached portion of this slab is located beneath the Adriatic Sea (DiD?).
A south dipping slab is imaged in the central Alps (CA). In the east, two different, weak, high-velocity anomalies are present, one in the south, beneath the southern part of the eastern Alps, and another, in the north, beneath the northern part of the eastern Alps and the Molasse Basin. Both anomalies extend only down to 150-200 km depth. We interpret this as the northward subduction of the Adriatic Slab in the south and the southward subduction of a Eurasian slab segment in the north. In the western Alps, a shallow eastward dipping Eurasian slab segment (WA) is in close proximity to an almost vertically dipping attached slab segment beneath the northern Apennines (NA) and the southern Po Basin. A slab gap is found in the central Apennines (CAG). The northwest dipping Calabrian Slab (Cl) is attached in the area beneath Calabria and Sicily but detached beneath the southern Apennines. A detached portion of the Calabrian Slab is also found beneath northern Sicily. Beneath Algeria, the Kabylides Slab (Kb) appears to be attached to the North African lithosphere along the African coast. Detached part of the Kabylides slab seems to extend toward southern Sardinia. At least down to about 250 km depth, the Kabylides Slab is clearly separated from the Calabrian Slab (Cl) to the northeast and from the Gibraltar-Betic Slab to the west, the latter comprising two segments: a shallow Gibraltar (Gib) and a detached Betic one (Btc).
The geometry of the suggested inventory of slab segments represents an improvement compared to previous models because it is based on a very large set of measurments and on data sampling of the region denser than available previously. With continuing growth of seismic networks, it will be further refined by tomographic studies in the future. The slab geometries constrained by tomography provide essential inputs for quantitative geodynamic modeling of plate kinematics and plate deformation in the Mediterranean realm.

Data Availability Statement
Waveform data used in this study are publically available except for the Egyptian data. The model MeRE2020 as well as the phase velocity maps used to compute it are available on our institutional web page (https://www.seismologie.ifg.uni-kiel.de). In addition, the model will also be made available at the webpage of the IRIS Earth models repository (EMC, https://ds.iris.edu/ds/products/emc-earthmodels/). Obspy routines (Beyreuther et al., 2010) have been used to download waveform data.