Two subduction-related heterogeneities beneath the Eastern Alps and the Bohemian Massif imaged by high-resolution P-wave tomography

We present high-resolution tomographic images of the upper mantle beneath the Eastern Alps and the adjacent Bohemian Massif (BM) in the north based on recordings from the AlpArray-EASI and AlpArray seismic networks. The tomography locates the Alpine high-velocity perturbations between the Periadriatic Lineament and the Northern Alpine Front. The northward-dipping lithosphere is imaged down to ∼ 200–250 km of depth, without signs of delamination. The small amount of crustal shortening compared to that in the Western Alps and the bimodal character of the positive perturbations with a separation beneath the Tauern Window indicate a dual source of the velocity heterogeneity, most probably formed by a mixture of a fragment of detached European plate and the Adriatic plate subductions. A detached high-velocity heterogeneity, sub-parallel to and distinct from the Eastern Alps heterogeneity, is imaged at ∼ 100–200 km beneath the southern part of the BM. We associate this anomaly with the western end of a SW–NE-striking heterogeneity beneath the south-eastern part of the BM, imaged in models of larger extent. The strike, parallel with the Moldanubian–Brunovistulian mantle–lithosphere boundary in the BM and with the westernmost part of the Carpathian front, leads us to consider potential scenarios relating the heterogeneity to (1) a remnant of the delaminated European plate, (2) a piece of continental-and-oceanic lithosphere mixture related to the building of the BM, particularly to the closure of the old Rheic ocean during the MD–BV collision, or (3) a lithospheric fragment going through to the NW between the Eastern Alps and Western Carpathians fronts in a preceding subduction phase. The study is dedicated to our outstanding and respected colleague Vladislav Babuška, who coined innovative views on the European lithosphere and died on 30 March 2021.

1 Introduction

Teleseismic body-wave tomography represents a powerful tool to study regional velocity structure of the upper mantle and to image velocity anomalies, particularly those related to subducted plates in collision zones. The Alps have developed at a collision zone of the Eurasian and Adriatic plate since the Variscan orogeny (Fig. 1). The classical concept assumed that the European lithospheric slab subducted south-eastward–southward along the entire Alps (Laubscher, 1970; Mueller, 1982) without any segmentation. However, interactions of the European lithosphere with the translating and rotating Adriatic plate and several micro-plates involved in the collision, tearing, and retreat of the slabs resulted in the bent (arcuate) shape of the western Alpine mountain range on the surface and in the complicated geometry of the rigid lithosphere penetrating the ductile mantle.

Figure 1 Elevation map of the western part of the Bohemian Massif and the Eastern Alps, with the EASI and AASN seismological stations (a). The region is marked in light brown along with Variscan massifs (dark grey patches) in (b), also presenting the large-scale tectonic context. The grid for the modelling structure of the region (c) shows inverted and non-inverted nodes filled in dark green and in yellow, respectively. Labels of main tectonic faults and units: PAL – Periadriatic Lineament, NAF – Northern Alpine Front, SAF – Southern Alpine Front, SEMP – Salzach–Ennstal–Mariazell–Puchberg fault, TW – Tauern Window (dashed), DF – Danube fault, Id – Idria fault, MD – Moldanubian part of the Bohemian Massif (BM).

Early tomographic models of Europe (Aki et al., 1977; Spakman, 1990) have been replaced with more advanced models of segmented Alpine slab during the last decades (Wortel and Spakman, 2000; Kissling et al., 2006; Malusa et al., 2021, for reviews). For the first time, Babuška et al. (1990) imaged two segments of Alpine slab, one dipping SE–S beneath the Western Alps and another dipping N beneath the Eastern Alps, with a gap between them. The density of stations and teleseismic rays enabled resolving the high-velocity heterogeneities only on a 1.5^∘ × 1.5^∘ grid, but the bent shape of the south-eastward-dipping subduction of the European plate in the Western Alps and steep northward-dipping lithosphere beneath the Eastern Alps, with a gap between the two Alpine slabs, were evident. More recent tomographic studies from data recorded in regional passive seismic experiments with densely spaced stations resolved the Alpine subductions at finer grids and confirmed the suggested change in polarity between the distinct western and eastern Alpine lower lithosphere roots (Lippitsch et al., 2003; Mitterbauer et al., 2011; Karousová et al., 2013; Zhao et al., 2016; Hua et al., 2017). Standard isotropic regional velocity tomography (e.g. Piromallo and Morelli, 2003; Koulakov et al., 2009) based on pre-AlpArray data (Hetényi et al., 2018a) imaged the south-eastward-dipping curved slab of the Eurasian lithosphere in the Western Alps and the northward-dipping plate beneath the Eastern Alps of similar sizes (geometry), though interpretations are different. Dando et al. (2011) interpret high-velocity heterogeneities at the bottom of their regional tomographic model of the upper mantle beneath the Alpine–Pannonian region (Lombardi et al., 2009; Hetényi et al., 2009) as a graveyard of relict downwelling of continental lithosphere.

While resolution of early tomography of the Alpine region allowed researchers to deal only with the large and most distinct heterogeneities in the upper mantle, accumulation of high-quality data from the dense AlpArray network motivated us to search finer images of the upper mantle and to search for answers on segmentation of the Alps, dip directions of the subductions, their relevance to the European or Adriatic plate, extent of slab delamination, and particularly to elucidate the smaller size heterogeneity beneath the Bohemian Massif north of the Eastern Alps subduction. We concentrate in this paper on imaging the Europe–Adria plate collision across a transect in the Eastern Alps (EASI), aiming at understanding orogen-forming processes, and present high-resolution tomographic images of the upper mantle beneath the Eastern Alps and the adjacent Bohemian Massif (BM) in the north. The upper mantle beneath the BM is slow in general (e.g. Piromallo and Morelli, 2003; Amaru, 2007; Koulakov et al., 2009; Fichtner and Villaseñor, 2015). Based on data from the AlpArray-EASI (Hetényi et al., 2018b) and AlpArray Seismic (AASN) networks (AlpArray Seismic Network, 2014, 2015; Hetényi et al., 2018a), our tomography shows a high-velocity heterogeneity of small size at ∼ 100–200 km beneath the south-eastern part of the BM (referred to as HV-BM throughout the paper), sub-parallel to and distinct from the Eastern Alps high-velocity heterogeneity (labelled HV-EA). Kästle et al. (2018) identified in their surface-wave tomography approximately 1^∘ to the north of the Periadriatic Fault a high-velocity heterogeneity, similar to HV-BM, but further to the west. Considering the NE continuations of the HV-BM as imaged in body-wave tomography of a larger extent (e.g. Karousová et al., 2013; Paffrath et al., 2021), the heterogeneity strikes with the SW–NE trend, in parallel with the boundary of the Moldanubian (MD) and Brunovistulian (BV) mantle–lithosphere in the BM and the westernmost part of the Carpathian front. Finding a secondary slab and/or thickened and cold lithosphere north of the Alps strengthens our motivation to determine its proper location and discuss hypotheses regarding its origin. Besides linking the shallow most probably dual-source heterogeneity beneath the Eastern Alps to subduction of the Adriatic plate, we also present and discuss in our paper three potential scenarios of the origin of the positive heterogeneity located beneath the southern Bohemian Massif (HV-BM).

2 Data

High spatial density of stations involved in passive seismic experiments, high-quality recorded data, and dense ray coverage of the upper mantle under study are essential prerequisites for reliable high-resolution tomographic imaging. The AlpArray passive seismic experiment, realized in a broad European cooperation (Hetényi et al., 2018a), provided the necessary high-quality recordings for such a study. We collected recordings from stations of the AlpArray Seismic Network (AASN, https://doi.org/10.12686/alparray/z3_2015) and AlpArray-EASI network (https://doi.org/10.12686/alparray/xt_2014). All the AASN stations installed in a 200 km wide band (Fig. 1) along the densely spaced stations of the AlpArray-EASI complementary experiment (Hetényi et al., 2018b) were selected for this study. The N–S band of the EASI stations is oriented perpendicularly to the eastern Alpine chain (crest) and runs through the BM in the north to the Adriatic Sea in the south with a length of ∼ 540 km.

We have applied several procedures to check the data quality (Vecsey et al., 2017), particularly data completeness and correct timing, to eliminate periods with the uncorrected leap seconds or failures of clock synchronizations. In station metadata files we fixed, for example, wrong channel naming, station mislocation, and incorrect units for poles and zeros (mixing Hz and radians per second as units). We collected seismograms from 1920 earthquakes recorded at 240 temporary and permanent stations involved in the AlpArray experiments, from which we selected a subset of 201 top-quality earthquakes from epicentral distances greater than 30^∘ with as uniform as possible distribution relative to the region of the Eastern Alps (Fig. S1). From this dataset we picked coherently teleseismic P-wave arrival times with a fully automatic picker, TimePicker 2017 (Vecsey, 2021), developed in the ObsPy/Python platform (Krischer et al., 2015). The TimePicker 2017 is based on two-step signal cross-correlations and allows us to measure absolute arrival times (Fig. S2a). Instead of a subjective selection of a reference trace, we cross-correlate all pairs of traces in the first step. For each trace a time shift related to maxima of the cross-correlation function is determined. Traces at individual stations are then shifted by a weighted average of the time shifts gained from cross-correlations with all other stations which recorded the event. The low-noise beam trace, created as a stack of the shifted traces at all stations, forms the reference trace for performing new beam-forming and arrival-time picking of an event in the second cross-correlation step.

In the second step, we correlate traces of all stations with the reference low-noise beam trace. The waveforms are aligned according to the times related to the maxima of the cross-correlation functions. The new low-noise beam is computed as a median of all aligned traces, and the P-wave onset on the beam is determined automatically (see Fig. S2). Then, arrival times at individual stations are derived from differences between times of corresponding extremes. Arrival times on the station signals are measured by three different methods. The final times of individual extremes (green P1 and P2 in Fig. S2a) and their error estimates are computed from the normal distribution, which approximates a mixture of normal distributions of partial picks (for details see Fig. S2a). We have applied the picker to the full dataset of 1920 teleseismic events recorded by the AASN. Figure S2b shows uncertainties of the measured P-wave arrivals as well as means and medians for both the complete dataset and for events selected for tomography (see below).

Each earthquake in the subset of selected 201 earthquakes was recorded at least by 50 stations of the array, i.e. at least 20 % of all the stations in the area. In this dataset, 130 rays per event sampled the mantle on average. The conditions ensure sufficient stability of the reference level in computing the relative travel-time residuals. The north–south elongated shape of the region oriented across the Alpine structures and perpendicular to the strike of the presumed subduction does not guarantee the same resolution along and across the strike. To eliminate mapping effects of heterogeneities alongside the model into its internal part and to enhance the resolution in the direction of the subductions, we selected additional rays coming from the northern and southern 60^∘ wide azimuth bins. Only rays propagating through the model within the two azimuthal fans (see Fig. S1b) have been included in the final tomographic inversions, tested, discussed, and interpreted further. These data come from 244 earthquakes, each of them being recorded by 120 stations on average, i.e. by 50 % of stations in the region. Figure S5 demonstrates that potential bias coming from heterogeneities west and east of the N–S-oriented elongated array is weak and that the events from the 60^∘ cones better illuminate the Bohemian Massif (BM)–Eastern Alps (EA) structures.

Teleseismic data cannot resolve velocities in the crust itself due to their sub-vertical propagation at shallow depths. To avoid mapping effects from the crust into the velocity perturbations in the upper mantle (e.g. Karousová et al., 2012), one has to introduce crustal corrections. Unfortunately, up to now, there has been no uniform, sufficiently detailed model of the crust for Europe or for the AlpArray region. For our body-wave studies in Europe we have collected accessible information on the crust for each station from different sources: from Karousová et al. (2012, and references therein) for the BM mostly based on results of controlled source seismics; from e.g. Di Stefano et al. (2011), Hua et al. (2017), and Tesauro et al. (2008) for areas south of the BM; and from Hetényi et al. (2018b) along the EASI transect. The crust is characterized beneath each station by the depth of Moho, the dip angle, and dip direction if the Moho is not flat, the Moho jump if there is any, velocity in the crust, thickness of sediments, and sediment velocity (Fig. S3). Sometimes there are significant differences between different models in their overlapping parts; therefore, we do not attempt to create any kind of “fine”-gridded model of the crust. Instead, we tune the corrections individually beneath each station and correct the travel-time residuals along each individually traced ray for the difference in the “real” crust and the crust of the reference model. Carefully pre-processed P-wave travel-time residuals calculated relative to the IASPEI'91 velocity model (Kennett and Engdahl, 1991), corrected for the crust, normalized by the average residual per event, and cleaned from outliers serve as input to the inversion in which we do not invert for the crust. With this approach we gathered a high-quality and uniform dataset of travel-time residuals for a proper tomographic inversion to resolve structures in our target region and below the crust. Of course, all crustal models remain idealized models and are not 100 % correct. Therefore, small “uncorrected” perturbations remain in each teleseismic tomography model of the upper mantle, which one needs to consider when interpreting the results.

3 Methods

We retrieved the velocity perturbations in the upper mantle with the isotropic mode of the coupled anisotropic–isotropic tomographic code AniTomo (Munzarová et al., 2018a) derived from the broadly used Telinv code (e.g. Weiland et al., 1995; Arlitt et al., 1999; Lippitsch et al., 2003; Sandoval et al., 2004; Shomali et al., 2006; Eken et al., 2007; Karousová et al., 2012, 2013; Plomerová et al., 2016; Silvennoinen et al., 2016; Chyba et al., 2017).

Weak anisotropy with hexagonal symmetry, with the high-velocity a axis or with the low-velocity b axis generally oriented in 3D, are assumed in the coupled anisotropic–isotropic tomographic code AniTomo. Velocity at each point can be expressed as

$$\begin{array}{}\text{(1)}& v=\stackrel{\mathrm{‾}}{v}\left(\mathrm{1}+{\displaystyle \frac{k}{\mathrm{2}}}\cos \mathrm{2}\mathit{\alpha }\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{v}$ is the isotropic component of anisotropic velocity, k is strength of anisotropy, and α is an angle between the symmetry axis and wave propagation direction (for details see Munzarová et al., 2018a). The linearized relation between a travel-time residual Δt and perturbations of the four anisotropic parameters $\mathrm{\Delta }\stackrel{\mathrm{‾}}{v}$ (isotropic component of the anisotropic velocity), Δk (strength of anisotropy), Δλ, and Δθ (azimuth and inclination of the symmetry axis) at each grid node (indexed with i) attains the form

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{\Delta }t& ={\sum }_i{\left({\displaystyle \frac{\partial t}{\partial \stackrel{\mathrm{‾}}{v}}}\right)}_i\mathrm{\Delta }{\stackrel{\mathrm{‾}}{v}}^i+{\sum }_i{\left({\displaystyle \frac{\partial t}{\partial k}}\right)}_i\mathrm{\Delta }{k}^i\\ & +{\sum }_i{\left({\displaystyle \frac{\partial t}{\partial \mathit{\lambda }}}\right)}_i\mathrm{\Delta }{\mathit{\lambda }}^i+{\sum }_i{\left({\displaystyle \frac{\partial t}{\partial \mathit{\theta }}}\right)}_i\mathrm{\Delta }{\mathit{\theta }}^i.\end{array}\end{array}$$

Due to the elongated shape of the region (Fig. 1), which is not suitable for coupled anisotropic–isotropic tomography and is out of the scope of this paper, we apply to the AlpArray-EASI data only the isotropic mode of the AniTomo code, in the first step. Then Eq. (2) reduces to the relation between the travel-time residual and the isotropic velocity perturbations:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }t={\sum }_i{\left({\displaystyle \frac{\partial t}{\partial \stackrel{\mathrm{‾}}{v}}}\right)}_i\mathrm{\Delta }{\stackrel{\mathrm{‾}}{v}}^i.\end{array}$$

The system of linear equations is then solved with the standard damped least-square method (e.g. Menke, 1984):

$$\begin{array}{}\text{(4)}& \mathit{m}={\left({\mathbf{A}}^T{\mathbf{W}}_D\mathbf{A}+{\mathit{\epsilon }}^{\mathrm{2}}{\mathbf{W}}_{\mathbf{M}}\right)}^{-\mathrm{1}}{\mathbf{A}}^T{\mathbf{W}}_{\mathbf{D}}\mathit{d},\end{array}$$

where m is a vector of model parameters $\mathrm{\Delta }\stackrel{\mathrm{‾}}{v}$ at all nodes. Data vector d contains travel-time residuals Δt and matrix A stores the partial derivatives from Eq. (2) or (3). Errors of arrival-time measurements are considered in weighting matrix W_D. Damping factor ε² stabilizes the ill-posed problem. Horizontal smoothing of model parameters can be achieved via matrix W_M. There is no vertical smoothing in the code. The inverse in Eq. (4) is approximated by truncated singular value decomposition. A 3D ray-tracing bending technique called Simplex (Steck and Prothero, 1991), in which ray paths are distorted by sinusoidal signals, is applied. Reliability of the model parameters for a given ray distribution and inversion setup can be assessed with resolution matrix R:

$$\begin{array}{}\text{(5)}& \mathbf{R}={\left({\mathbf{A}}^T{\mathbf{W}}_D\mathbf{A}+{\mathit{\epsilon }}^{\mathrm{2}}{\mathbf{W}}_M\right)}^{-\mathrm{1}}{\mathbf{A}}^T{\mathbf{W}}_D\mathbf{A}.\end{array}$$

The area of about 400 000 km², centred at 13.3^∘ E, 48.5^∘ N, is approximated by a 30 × 30 km cell size horizontally. The images are calculated down to 435 km of depth with 30 km spacing. To minimize creating false perturbations, we invert for the velocity perturbations only in the central $\mathrm{5}\times \mathrm{25}\times \mathrm{13}$ cells, which are well-sampled by crisscrossing rays (Fig. S1), i.e. in nodes between depths of 60 and 420 km. The model covers the Eastern Alps and a core of the BM, an area of ca. 140 400 km² in total. Variance reduction of the final model for the chosen damping parameter reaches 66 % (Fig. S4), and we kept the results after the second iteration (Fig. S4) as a further iteration has not resulted in any significant change in the velocity perturbations.

4 Results

The distinct, high-velocity, northward-dipping, ∼ 140 km broad perturbations related to the eastern Alpine root, imaged in the upper ∼ 250 km of the mantle by previous tomography (Babuška et al., 1990; Karousová et al., 2013; Hetényi et al., 2018b), had a tendency to split when we exploited data from the EASI experiment and nearby permanent stations (Plomerová et al., 2018). However, “only” adding data from the AASN led to the clear visualization of two separate sub-parallel high-velocity heterogeneities beneath the broader eastern Alpine region, both dipping to the north and each about 80 km thick (Figs. 2 and 3), with a low-velocity separation zone of ∼ 80–100 km extent.

A clear decrease in amplitudes with depth dominates in the horizontal depth slices (Fig. 2a) through the EASI-AA velocity perturbation model. They exceed ± 1 % only exceptionally below 220 km of depth (two deepest layers shown in Fig. 2). Negative perturbations in the two uppermost mantle layers concentrate along the Eger Rift (ER) and can be related to the lithosphere thinning in this region relative to the MD part of the BM (e.g. Plomerová and Babuška, 2010; Plomerová et al., 2016). At greater depth, the lower velocities dominate in the sub-lithospheric mantle beneath the whole BM (e.g. Amaru, 2007; Fichtner and Villasenor, 2015). The distinct positive velocity perturbations related to the eastern Alpine root are located north of the PAL at ∼ 47^∘ N and are distinct down to 225 km of depth. At greater depths both the positive and weak negative perturbations are arbitrarily mixed and do not indicate any continuous object. A relatively smaller-sized positive velocity heterogeneity, north of the strongest large one, lies beneath the southern BM at 100–200 km of depth.

Figure 2 Depth slices through the new EASI-AA velocity perturbation model (a) along with depth contours of the 1.5 % perturbations in map view (b). The arrow marks the northward dip of the central and eastern part of the eastern Alpine heterogeneity. Labels of main tectonic faults and units: PAL – Periadriatic Lineament, NAF – Northern Alpine Front, SAF – Southern Alpine Front, SEMP – Salzach–Ennstal–Mariazell–Puchberg fault, TW – Tauern Window (dashed), DF – Danube fault, Id – Idria fault, MD – Moldanubian and TB – Teplá–Barrandian units of the Bohemian Massif (BM), ER – Eger Rift.

To have a better sense of both lateral and vertical changes in the perturbations and thus the dip of the high-velocity heterogeneities, we contour the 1.5 % positive perturbations between 90 and 210 km of depth (Fig. 2b). The contours clearly mark the northward dip of the EA slab, particularly in its eastern part (east of 13.3^∘ E). A dip of the western rim of the slab is not clear from the contour curves only, as it appears to become steeper and thins significantly with depth. Similarly, it is difficult to judge a dip direction of the smaller-sized positive heterogeneity around ∼ 48.7^∘ N (HV-BM) in this visualization.

To image the dip directions of the heterogeneities, we present five N–S vertical cross sections through the EASI-AA model (Fig. 3a–c), perpendicular to the strike of the mountain belt. The images of slab geometry suggest changes along-strike of the Alpine orogen, even at short distances, but one has to keep in mind changes in resolutions toward the margins of the model as well (see sections with synthetic tests). In general, the positive perturbations reach down to ∼ 220 km. All perturbations below ∼ 250 km of depth are very weak without any clustering or evident association with the stronger heterogeneities above this level. The positive perturbations beneath the BM, mapped only between ∼ 100 and ∼ 220 km depths (marked as III in Fig. 3), are slightly weaker in comparison with those beneath the EA (marked as I and II in Fig. 3) and disconnected from shallow parts. On the other hand, the strongest heterogeneity in the south does not exhibit clear signs of detachment, and thus through its dip and connection to crustal levels an association with a continuous subduction of the Adriatic plate could be possible. The limits of the numerically obtained well-resolved area, derived from resolution matrix, are reported in the figures with grey shading.

Figure 3 Five north–south-oriented vertical cross sections through the EASI + AA velocity perturbation model: (a) east of the EASI, (b) along the EASI profile, and (c) west of the EASI. The location of the cross sections is marked in the map in the central part of the figure. The along-longitude cross sections run from 51.65 to 45.35^∘ N. Fault abbreviations are as earlier: Eu – Europe, Ad – Adria. Shading of the less well-resolved region was set according to the resolution matrix, ray coverage of grid cells, and derivative weighted sums. The contour follows a value of 0.15 of diagonal elements of the resolution matrix. Green dots mark Moho depths in the model used for calculation of crustal corrections (see also Fig. S3). Travel times are not inverted for velocity perturbations in the crust (nodes at 30 km of depth).

The positive velocity perturbations seem to penetrate the sub-lithospheric mantle in the northward direction, in general, at an apparent dip of ∼ 45^∘ or more. The general dip of these heterogeneities marked in the cross sections changes only slightly in direction toward the Central Alps. In the central cross sections, the HV-EA heterogeneity appears shorter than in the easternmost cross section, while the detached lithosphere fragment beneath the BM reaches slightly deeper. Cross sections through a 3D visualization of along-strike changes in the velocity perturbation in the EASI-AA model can be found in Fig. S9.

Looking at perturbations of the HV-EA heterogeneity at greater details, one can recognize its potential dual source. The positive heterogeneities I and II (Fig. 3) are separated just beneath the TW, where we modelled the European and Adriatic crust contact (Hetényi et al., 2018b). Dip directions of the heterogeneities slightly change – HV-EA-part II seems to immerse southward (Fig. 3a, b), whereas HV-EA-part I dips to the north. This allows us to argue for a mixing of a detached EU slab fragment and the shallow Adria slab connected at depth, as already suggested in Babuška et al. (1990) using a much coarser model inferred from P residuals. The length of the subduction imaged in our tomography beneath the Eastern Alps is similar to that in Lippitsch et al. (2003) and fits the upper limit of the amount of Adria–Europe Neogene convergence described in Le Breton et al. (2017).

5 Resolution tests

We have performed several synthetic tests to evaluate the resolution of the tomography results, particularly its ability to detect the two separate sub-parallel slabs beneath the Eastern Alps and BM as well as their dip direction. The polarity reversal of the northward subduction beneath the E Alps relative to the south-eastward subduction in the Western Alps is of particular importance. The polarity flip is still questioned by some authors (e.g. Kind et al., 2021, this issue) in spite of long-lasting various inferences speaking for the change in subduction polarity beneath the Western and Eastern Alps with a gap between them (e.g. Babuška et al., 1990; Lippitsch et al., 2003; Zhao et al., 2016; Paffrath et al., 2021, this issue, and references therein).

Test 1 (Figs. 4, S6a–d) was designed to compare data-retrieved perturbations with those resulting from one or two narrow vertical synthetic heterogeneities, without imposing any polarity of the subductions. The model with one 5 % velocity heterogeneity does not reproduce the velocity perturbations retrieved from real data. The model with two steep heterogeneities mimics the perturbations much better in both the central part of the model and its margins. The perturbations retrieved from the synthetic vertical heterogeneities remain vertical for the real ray geometry or with a weak southward-dipping tendency in the westernmost profile, which contradicts the northward dip of real perturbations. Evidently, there is no northward smearing due to the ray geometry.

Having documented the existence of two separate heterogeneities, one beneath the southern BM (HV-BM) and the second beneath the Eastern Alps (HV-EA), we have tested their relative orientations. Test 2 (Figs. 4, S6a–d) assumes two heterogeneities as in Test 1, but with 27^∘ southward dip. The resulting perturbations do not reproduce the northward dip of the real perturbations. On the other hand, Test 3 with the two heterogeneities dipping to the north at 27^∘ mimics the dip of real perturbations very well. Test 4 with two bivergently dipping heterogeneities (towards each other) matches the geometry of the real perturbations only at shallow depths above ∼ 150 km, but the deeper part of the northern heterogeneity is completely missing.

The synthetic tests performed above corroborate the fact that the data from the AlpArray and EASI networks are able to image two separate northward-dipping sub-parallel slabs beneath the Eastern Alps and southern rim of the BM. The two slabs are separated from each other, and the northern one is not connected with the shallow parts of the lithosphere (above ∼ 100 km). The difference between the dip directions of the subductions beneath the Eastern Alps and the subductions beneath the Western and Central Alps is undoubtedly real, and it is not produced by potential smearing due to ray geometry.

Besides the specific tests described above we also performed standard checkerboard tests to assess the resolution capability of the network (Fig. S7). The checkerboard test confirms that positive and negative perturbations are retrieved well down to at least 240 km with a weak vertical smearing (Fig. S7a – horizontal slices). Also, the vertical cross section through the central part of the model (Fig. S7b) images the synthetic perturbations reliably.

6 Imaging the high-velocity perturbations in different tomography models

The present shape of the Alpine mountain chain, characterized by the curved Western Alps and east–west-striking Central and Eastern Alps, reflects a multi-phase action of tectonic forces during the collision of the European and Adriatic plates, the AlCaPa micro-plate, and numerous lithosphere fragments in the regions, as well as the Piemont oceanic lithosphere. The processes are imprinted in the complex architecture of the broader Alpine region in both the crust (Handy et al., 2010; Rosenberg and Kissling, 2013; Schlunegger and Kissling, 2015) and the mantle (e.g. Kissling et al., 2006). Continuing debates on the exact setting of the Moho depth in the Alps and on the “gap” near the Tauern Window (Spada et al., 2013; Hetényi et al., 2018b; Brückl et al., 2010) document the complex structure of the Alpine orogen.

Though the upper mantle structure is less diverse in comparison with the crust, in general, ongoing studies of the Alpine upper mantle continue to reveal new and more detailed features in the geometry of the lower lithosphere, dip direction of the Alpine slabs, tears or detachments of the slabs, and interactions of the Alps with the Apennines and Dinarides. The current stage of knowledge from results of various disciplines – seismology, geology, petrology, tectonics, paleo-magnetism, geochemistry, and GPS studies, among others – reflects differences in the segmented slab responses to the acting forces. The complex structure of the fragmented Alpine slab(s) and the broader Europe–Adria collision zone is now visualized in tomography snapshots. The boundary between the central and eastern slab segments of the Alps is placed at the western end of the Tauern Window, i.e. at the northward prolongation of the Giudicarie–Brenner Fault, as already proposed by Handy et al. (2015) and other more recent authors to separate different crustal (lower crustal wedges) and mantle (slab gaps, slab dips) structures (e.g. Lippitsch et al., 2003; Schmid et al., 2004; Rosenberg et al., 2018). In recent studies, Paffrath et al. (2021) and Handy et al. (2021) document a change from a S-dipping slab to a N-dipping slab beneath the western part of the Eastern Alps (i.e. beneath the Tauern Window). This is in agreement with the previously proposed geometries of Babuška et al. (1990), Lippitsch et al. (2003), and Zhao et al. (2016). Mock et al. (2020) note a discordance between the slab geometry at depth and the boundary between the Eastern and Central Alps observed at the surface.

Figure 4 Synthetic tests of tomography capability to resolve one or two sub-parallel heterogeneities (TEST1) and their dip directions (TEST2, TEST3, and TEST4) in the central cross section along the EASI. The similar north–south cross sections parallel to the central EASI profile on its east and west, corresponding to profiles in Fig. 3, are in Supplement Fig. S6a–d. The along-longitude cross sections run from 51.65 to 45.35^∘ N. Fault labels are as in Fig. 2.

As early as 1990 Babuška et al. (1990) suggested segmentation of the Alpine slab into the western and eastern parts with opposite polarities, as well as a gap between them, in the tomography of central Europe, in which the authors inverted crust-corrected and source-side clustered travel-time residuals. This model (Fig. 5) has rather high variance reduction of 70 % but a resolution of only 1.5^∘ × 1.5^∘ due to station spacing available at that time. The Alpine slab segmentation has been confirmed in the upper mantle tomography of the Alpine orogen by Lippitsch et al. (2003). Their regional tomography with a finer lateral grid of 50 × 50 km has the highest resolution around 12^∘ E, just between the two strongest high-velocity perturbations beneath the Central and Eastern Alps, where the passive seismic network TRANSALP was deployed. Later tomography of the upper mantle, which included the Eastern Alps with data from regional passive experiments (Dando et al., 2011; Mitterbauer et al., 2011; Karousová et al., 2013), also retrieved the northward-dipping high-velocity heterogeneity of similar geometries (Fig. 5) in 250 km of the upper mantle. The imaged triangular shape of the LAB model (e.g. Babuška and Plomerová, 1992) beneath the thickened lithosphere of the Eastern Alps, which is detached from that beneath the Western and Central Alps, led Babuška et al. (1990) to suggest three main phases in building the Eastern Alps lithosphere root: (1) NW translation of the Adria and its thrusting over the (subducting?) European plate in the Alps, (2) fragmentation of the northern Adria along a deep-seated fault (possibly the Giudicarie Fault or at least a spatially nearby structure), and (3) counter-clockwise rotation of the Adria resulting in a start of its subduction after a collision with the European plate in the Eastern Alps (e.g. Ustaszewski et al., 2008), potentially above a delaminated European lithosphere residing at greater depths (Paffrath et al., 2021; Handy et al., 2021). The complexity of this region is reflected in ambiguous views on the crust, which is interpreted as a triple junction of three crustal terranes (Brückl et al., 2010), although the deformation style between the Eastern Alps and the Pannonian Basin is usually considered diffuse on the surface. Instead, e.g. Ustaszewski et al. (2008) and Horváth et al. (2015) assume that the Pannonian Basin is underlain by gradually thinned European and Adriatic crust, which has been recently imaged by Kalmár et al. (2021).

The most recent tomography of the entire Alps and surrounding regions (Paffrath et al., 2021; Handy et al., 2021 – both this issue) exploits data from the AlpArray seismic network and AlpArray complementary experiments. The large-scale AlpArray tomography and the EASI-AA model along the EASI with similar gridding exhibit remarkable coincidence of perturbation patterns. In both tomography images, the HV-EA is located between the PAL in the south and the NAF in the north, but the large-scale tomography (Paffrath et al., 2021, this issue) images the Eastern Alps subduction below 150 km of depth with low-velocity perturbations (red) above it. At this depth, horizontal slices show characteristic separation of the Western, Central, and Eastern Alps, which becomes more distinct deeper in the mantle. Similarly, Zhao et al. (2016) show only weak positive perturbations at 100 km beneath the Eastern Alps, but this region lies in a relatively less well-resolved part of their model.

Figure 5 Cross sections through various regional teleseismic P-wave tomography results showing segmentation of the Western and Eastern Alps and changes in the slab dip directions (a, b), the northward dip of the steep slab (HV-EA) beneath the Eastern Alps (b–f), and the small-sized HV-BM heterogeneity beneath the Bohemian Massif (c, f, g). Cross sections (a, b) from Babuška et al. (1990), with the crust in (b) from Hetényi et al. (2018b) and perturbation contours in (b) from Karousová et al et al. (2013). The profile locations in map view are in the upper right. The dark and light blue highlights refer to the location of the positive perturbations beneath the Eastern Alps and the Bohemian Massif (HV-EA and HV-BM), respectively. The full set of colour palettes can be found in the original publications.

A local discrepancy between models is present in the northern BM, where the large tomography returns positive perturbation, but the EASI-AA tomography maps low velocities there in accordance with other studies (e.g. Plomerová et al., 2016). However, the large-scale and EASI-AA tomography detects positive perturbations beneath the Moldanubian part of the southern BM. The perturbations continue further to the south-west in the EASI-AA in comparison with Paffrath (2021) tomography. Similar local high-velocity perturbations down to ∼ 250 km of depth, isolated from the EA heterogeneity, have also been detected in body-wave tomography by Dando et al. (2011) and Karousová et al. (2013) (Fig. 5) as well as in surface-wave tomography by Kästle et al. (2018, 2020), though the surface waves localize the similar high-velocity heterogeneity further south-west of the BM.

The north-oriented dip of the EA subduction, imaged in early tomography studies (Babuška et al., 1990; Aric et al., 1989), was questioned for a long time, until Lippitsch et al. (2003) clearly imaged the northward-dipping structures using TRANSALP data. Subsequent regional passive experiments provided body-wave data for tomography of the Eastern Alps and surrounding regions (namely in the Bohemian Massif and Pannonian Basin – BOHEMA, e.g. Karousová et al., 2013; Plomerová et al., 2016; ALPASS, e.g. Mitterbauer et al., 2011; and CBP, e.g. Dando et al., 2011, with passive experiments) and imaged the northward dip of the EA slab as well. Some of them agree in interpreting it to be of Adriatic plate origin (Karousová et al., 2013). However, Mitterbauer et al. (2011), and similarly in recent tomography based on the newest AlpArray data by Paffrath et al. (2021, this issue), the positive perturbations below 150 km are associated with a delaminated EU slab (see also Handy et al., 2021, this issue). Kästle et al. (2020) relate the HV-EA mainly with the European plate subductions as well and leave no role or only a minor role for the Adriatic subduction. The authors explain the northward-dipping subduction modelled beneath the Eastern Alps by imaging a combination of the short Adriatic and deep-delaminated, potentially overturned European slabs.

Synthetic tests from the EASI-AA dataset in our study proved that the EASI-AA model is capable imaging the two sub-parallel northward-dipping heterogeneities – the large and strong southern one beneath the Eastern Alps with connection to shallow depths and the weaker northern one beneath the southern BM that is apparently disconnected. To understand the positive perturbations beneath the southern BM, we compare it with results from the large-scale Paffrath et al. (2021) tomography and with the regional tomography of the BM (Karousová et al., 2013) of resolution similar to ours (Fig. 6). The strongest positive perturbations related to that heterogeneity overlap in the models, though they are of unrealistically large extent in the Paffrath et al. (2021, this issue) model. Paffrath et al. (2021) also show the strong positive perturbations beneath the north-western part of the BM (beneath the Eger Rift), where the lithosphere thins significantly (∼ 80 km). This is well-imaged by negative perturbations in the EASI-AA and in BOHEMA models (Karousová et al., 2013; Plomerova et al., 2016). Images of perturbations at 180 km of depth stress the importance of using data from stations in the northern part of the EASI array. Data from these stations capture the positive velocity perturbations much farther to the SW, beneath the southern rim of the BM, than the other two tomography studies. Cross sections through the BM regional tomography (see Fig. 5) locate the increased velocities within the low-velocity BM upper mantle (e.g. Amaru, 2007; Fichtner and Villasenor, 2015). The SW–NE elongated shape of the heterogeneity follows the boundary between the MD and BV mantle–lithosphere domains (Babuška and Plomerová, 2013).

Figure 6 Horizontal slices at 120 (a) and 180 km (b) depths through the velocity perturbation model EASI-AA (left) (this paper), the model by Karousová et al. (2013) (right), and the model by Paffrath et al. (2021, this issue) (centre), on which the HV-BM from the other two models is marked by orange contours. Labels are as in Fig. 2.

Different types of seismic waves propagate with different velocities and sample the mantle volume in different directions and wavelength, which affects the velocity–velocity perturbation images of the upper mantle and corresponding resolution. A cross section along 13.5^∘ E through the MeRE2020 model, which is a shear-wave velocity perturbation model from Rayleigh phase velocities (El-Sharkawy et al., 2020), runs parallel to the EASI transect in our P-wave model. The authors relate positive perturbations eastward of the cross section at ∼ 45^∘ N to the northernmost part of the Dinarides slab and image a short (only down to 150 km of depth) European slab, without any delamination. On the other hand, the body-wave tomography by Paffrath et al. (2021) sees the top of the high-velocity heterogeneity beneath the Eastern Alps at 150 km of depth. The heterogeneity extends down to ∼ 300 km of depth. The authors interpret it as the delaminated European plate lithosphere. There is an obvious contradiction between these two interpretations. To test whether a potential delamination can be detected in our EASI-AA tomography, we performed additional synthetic tests (Fig. S8). In each test the top of the steep 5 % high-velocity heterogeneity shallows from 150 to 60 km depths. Due to the leakage, which exists in each teleseismic P-wave tomography, one has to be careful and concentrate only on distinct features. Figure 8 demonstrates that the dark blue perturbations in the model with shallow heterogeneity reproduce the dark blue perturbations in the model from real data better than the model with heterogeneity below 150 km of depth. The tests demonstrate that slab detachment larger than the 30 km grid would be revealed in our EASI-AA model. Two northward-dipping slabs (Fig. 7), the northern small-sized delaminated one approximated with a 3 % velocity increase and the southern attached one with a 5 % velocity increase, mimic the velocity perturbations retrieved beneath the Bohemian Massif and beneath the Eastern Alps from real data well (Fig. 3). The leakage of the model perturbations of the HV-BM towards the surface is weak and acts against the negative perturbations imaged from real data. There is no surprise that a simple prism model cannot reproduce tiny details of the complex HV-EA heterogeneity.

Figure 7 Velocity perturbations along the five vertical along-longitude cross sections calculated from real data (upper row) and those (middle row) calculated from two synthetic 3 % and 5 % heterogeneities (bottom), plotted over the retrieved perturbations. Green dots mark Moho depths in the model used for calculation of crustal corrections (see also Fig. S3). The along-longitude cross sections run from 51.65 to 45.35^∘ N.

The inconsistency between the above-mentioned results in the upper ∼ 150 km of the upper mantle could come from different sources, e.g. a sub-optimal ray coverage in that area and depth range, but also can reflect differences in crustal corrections applied. A certain role can also be played by the difference in size of the arrays. Velocity perturbations in P-wave tomography are linked to averages within each layer (a plane of grid nodes), which might differ in a large array from that of a smaller-sized array. Data from large-scale arrays allow us to model heterogeneities down to greater depths than a smaller-sized regional tomography, which can, on the other hand, provide finer images. While tomography images presented here have the best resolution above 200 km and allow imaging a shallow slab beneath the Eastern Alps attached to Adriatic plate, the large-scale model from the entire AlpArray has the best resolution at depths 200–400 km (Paffrath et al., 2021, this issue) and is thus able to image the deep-delaminated European slab beneath the Eastern Alps.

The role of applying proper crustal corrections is significant in teleseismic regional tomography. Not applying any crustal corrections or applying inadequate ones can strongly affect velocity perturbations within the upper ∼ 100 km of the upper mantle (e.g. Karousová et al., 2013), which is the zone where the models discussed above differ. “Overcorrecting” the travel times due to the crustal model used, regardless of the method of correcting itself, can delete, erase, or substantially reduce positive perturbations in the upper 100 km if crust that is too slow and thick is considered and vice versa. From this point of view, developing a uniform, detailed, and reliable model of the European crust is urgently needed.

7 Potential scenarios for geneses of the high-velocity heterogeneities

Various evolution scenarios for the Eastern Alps slab exist, but there are none for the HV-BM beneath the southern BM. For detailed scenarios of the Alps subduction and the Europe–Adria plate collision we refer to Schmid et al. (2004), Handy et al. (2010, 2015), Le Breton et al. (2017), Kissling and Schlunegger (2018), Rosenberg et al. (2018), and Paffrath et al. (2021) as well as references therein. Their models of Alpine orogeny include subductions of the European and Adriatic plates, slab roll-back, tearing, break-offs and delamination, widespread intra-crustal and crust–mantle decoupling, and the NW translation and counter-clockwise rotation of the Adriatic plate. All these processes, accompanied by thermal erosion and/or deglaciation uplifts, are reflected in the different structure of the WA and EA. But how do we interpret the positive velocity perturbations within the low-velocity upper mantle beneath the southern BM (HV-BM)?

The elongated HV-BM is ∼ 80 km wide and extends, considering both the EASI-AA and Paffrath et al. (2021) model, approximately 12.5–16^∘ E over a length of ∼ 300 km on the surface. The elongated shape of this heterogeneity strikes in the SW–NE azimuth and extends from ∼ 100 km down to ∼ 200 km of depth. We estimate the total volume of the HV-BM at 1.5 million km³; these dimensions are comparable to a small lithospheric segment. The smaller part of the HV-BM imaged in the EASI-AA model extends at a low angle relative to the eastern Alpine front. But considering its full size, it is sub-parallel to the MD–BV contact in the BM as well as the Western Carpathians front. The high-velocity material hovers in the low-velocity upper mantle beneath the BM. We outline three potential scenarios as an explanation for he origin of this positive heterogeneity.

The simplest explanation would be to consider it as a fragment of the delaminated part of the European plate subductions, analogous to the scenario proposed by Handy et al. (2015) for the delamination of European lithosphere further to the south beneath the Eastern Alps (Fig. 8a). In their model, the second break-off or delamination of the European slab at ∼ 25 Ma opened space for the northward subduction of the Adriatic plate, pushed by Africa from the south and rotated by pulling due to subduction of its SE rim beneath the Hellenides. The delaminated piece of the continental lithosphere has continued sinking into the mantle in the model since then. However, the HV-BM is located at shallow depth (above ∼ 200 km, not really compatible with sinking) and is too far to the north (at ∼ 49^∘ N) from the Periadriatic Fault System. Therefore, an association of the HV-BM with the delaminated fragment of the EU subduction is not likely. Also, the clearly imaged separation (negative anomaly) between the subducting HV-EA and the HV-BM is a feature that would not be explained in this scenario.

Figure 8 Interpretation of the slab at shallow depth beneath the Eastern Alps (HV-EA), its Adriatic vs. European provenance, and potential scenarios for the origin of the heterogeneity beneath the south-eastern Bohemian Massif (HV-BM). (a) Cartoons of the European slab (EU) delamination at 20 Ma (upper left) and configuration of the Europe–Adria plate collision at present (lower right). Figure redrawn from Handy et al. (2015) (for the location of the profile see their Fig. 10) by overlaying the location of the high-velocity heterogeneities from the EASI-AA model in a cross section along 13.3^∘ E (Fig. 3b) and the model of Paffrath et al. (2021); (b) a scenario considering closure of the Rheic ocean and collision of the Brunovistulian micro-plate with the Moldanubian part of the BM in a portion of a schematic cartoon by Babuška and Plomerová (2013), with a piece of remnant lithosphere image as the HV-BM in this paper; (c) a scenario related to fragmentation of the Alpine and Carpathian front. Differences in the roll-back subductions of the Alps and the Carpathians (e.g. Royden and Burchfiel, 1989; Royden, 1993), northward push of Adria, and European slab delamination beneath the Eastern Alps (Handy et al., 2015) could have formed complex flows in the asthenosphere, which might “transport” lithosphere fragments toward the north-west. For more details see the main text. The high-velocity contours are drawn according to Fig. 3b.

The BM is an assemblage of fragments of continental lithosphere with their own large-scale anisotropic fabrics (e.g. Babuška and Plomerová, 2013, 2020). Changes in the fabrics delineate boundaries between the mantle–lithosphere domains. The location of the HV-BM, following the boundary of the MD–BV mantle–lithosphere domains, evokes a possible link to the MD–BV collision, which is related to the closure of the Rheic ocean in the late Devonian to the middle–late Carboniferous (Babuška and Plomerová, 2013, and discussion and references therein). The HV-BM is disconnected from the BM lithosphere in tomography cross sections (see e.g. Figs. 3, 5). The continual subduction of the oceanic plate due to negative buoyancy would have led to a removal of the denser materials from shallow depths since then. The Phanerozoic continental mantle–lithosphere may mechanically decouple from the upper or the full crust, since this lower lithosphere is denser than the asthenosphere into which it subducts. Densification of the lower crust through metamorphic reactions may enhance this process if the convergence is fast (e.g. 2 cm yr⁻¹ in Hetényi et al., 2011). The Phanerozoic continental mantle–lithosphere, composed of originally lighter rocks than those in the asthenosphere, becomes denser due to metamorphic phase changes as it subducts. This is the general process; however, at low convergence rates it is able to slowly return from the negative to positive buoyancy range (Boonma et al., 2019). Such a process could “stop” subductions and allow for long-term survival of this lithosphere material in the asthenosphere, with the high-velocity anomaly caused by chemical composition rather than temperature. Thus, the HV-BM could represent a remnant of the Rheic oceanic closure and/or a relic of the MD–BV collision, captured in the slow sub-lithospheric BM mantle (Fig. 8b).

A third scenario (Fig. 8c) could be some relatively light material brought to its current position beneath the BM not too long ago to survive there (i.e. not to sink rapidly). In this hypothesis the anomaly beneath the BM could be derived from subduction and lateral transport in the mantle of the oceanic embayment of the Alpine Tethys. To get such material there is difficult but not impossible considering the evolution of the European plate and the Carpathians. The shape of the originally linear, W–E-running Alpine–Carpathians front has changed since the time of the second European slab break-off, AlCaPa lateral escape, and beginning of the Adria subduction (Handy et al., 2015). The Carpathians front curved significantly and migrated to the north. Differences in the roll-back subductions of the Alps (e.g. Royden and Burchfiel, 1989; Royden, 1993) and the Carpathians, northward push of Adria, and European slab delamination beneath the EA could have formed complex flows in the asthenosphere (e.g. Vignaroli et al., 2008), which could “transport” a purely oceanic lithosphere or a mix of oceanic and continental lithosphere fragments through the open space between the eastern Alpine and Carpathian slabs north–north-westward into the mantle beneath the BM. The rotational displacement of the Adriatic plate indenter provided an additional driving force for modifications of the Alpine–Carpathian–Dinaric orogenic system (Ustaszewski et al., 2008). Remnants of the Penninic–Piemont oceanic lithosphere (Brückl et al., 2014 – Fig. 9 there), squeezed east of 14^∘ E between the AlCaPa and European plates in the Western Carpathians and pushed from the south by the Adria, can offer a possible explanation for the origin of the SW–NE elongated HV-BM striking sub-parallel to the Western Carpathian front.

Finding an unambiguous model of the complex Alpine orogeny and structure of the upper mantle in the broader surroundings of the Alps requires multi-method and interdisciplinary research that covers various spatial scales. The combination of gravity and seismic data represents one of such approaches (e.g. Lowe et al., 2021, this issue; Scarponi et al., 2021). Lowe et al. (2021, this issue) converted modified standard isotropic velocity perturbation models into velocity and then density models. Those are after that used to calculate the gravity signal, predicted up to 40 mGal for various slab configurations mimicking the Alps. The applied methods include severe simplifying assumptions. Nevertheless, neither including pre-defined slab geometries nor accounting for compositional and thermal variations with depth brings satisfactory results, which would allow them to distinguish between their two different slab configurations. The freshly compiled pan-Alpine surface-gravity database (Zahorec et al., 2021) will undoubtedly provide new impetus for structural investigations combining seismology and gravity. Regarding the gravity effect of the HV-BM, the expected signal is too weak to appear there clearly, as crustal effects are predominant in that area of the BM.

The anisotropic nature of the Earth has been proved to be a general characteristic in different seismological studies. Anisotropy influences mainly velocities and polarizations of seismic waves. Seismic anisotropy of the Earth's upper mantle carries key information for deciphering the tectonic history of the lithosphere–asthenosphere system (e.g. Babuška and Cara, 1991; Sobolev, 1999; Fouch and Rondenay, 2006; Long and Becker, 2010; Babuška and Plomerová, 2020, and references therein). However, effects of directional dependences of velocities are not considered in standard isotropic tomography images. Only long-wavelength shear-velocity models from surface waves traditionally include azimuthal and/or radial anisotropy in the mantle. Ignoring seismic anisotropy and assuming isotropic wave propagation or considering only azimuthal and/or radial anisotropy leads to significant isotropic and anisotropic imaging artefacts that may lead to spurious interpretations (van der Beek and Faccenda, 2021). In this study of the broader region around the Eastern Alps we have applied the isotropic mode of a coupled anisotropic–isotropic teleseismic P-wave tomography developed by Munzarová et al. (2018a). In spite of the general good agreement with the high-resolution large-scale isotropic tomography (Paffrath et al., 2021, this issue), the images can be biased due to seismic anisotropy (Eken et al., 2012; Qorbani et al., 2015, 2016; Bokelmann et al., 2021; Handy et al., 2021). Laterally varying anisotropy, which correlates with tectonics of the region, has been indicated in shear-wave splitting (e.g. Link and Rumpker, 2021). After collecting a sufficient amount of well-distributed high-quality data we will run the coupled anisotropic–isotropic mode of the code, which was successfully applied in northern Fennoscandia (Munzarová et al., 2018b). The resulting 3D anisotropic model of the region will map laterally and vertically varied anisotropy with symmetry axes oriented generally (i.e. inclined) in 3D. This further investigation may help decide among the drafted scenarios for the origin of the HV-BM or point to new ones.

8 Conclusions

The teleseismic P-wave tomography of the upper mantle beneath the Eastern Alps and the Bohemian Massif, as presented here, locates the Alpine high-velocity perturbations between the Periadriatic Lineament (PAL) and the Northern Alpine Front (NAF). The northward-dipping slab is imaged down to ∼ 200 km, without distinct signs of delamination from the Adriatic plate. The bimodal character of the positive perturbations with a separation beneath the Tauern Window indicates a dual source of the heterogeneity, most probably formed by a mixture of a fragment of detached European plate and the Adriatic plate subductions. The fine-gridded EASI-AA model of velocity perturbations images the individual high-velocity heterogeneity beneath the southern part of the Bohemian Massif at depths of ∼ 100–200 km. Its eastward continuation is visualized in other tomography results as well. Interpreting this heterogeneity as a remnant of the delaminated European plate seems unlikely. The SW–NE trend of the heterogeneity strike, in parallel with the Moldanubian–Brunovistulian mantle–lithosphere boundary in the Bohemian Massif or with the westernmost part of the Carpathian front, leads us to consider it to be a piece of a mixture of the continental and oceanic lithosphere related to building of the BM, particularly to the closure of the old Rheic ocean during the MD–BV collision, or a lithospheric fragment going through to the NW between the Eastern Alps and Western Carpathian fronts in a preceding subduction phase.