Constraints on Crustal Structure in the Vicinity of the Adriatic Indenter (European Alps) From Vp and Vp/Vs Local Earthquake Tomography

In this study, 3‐D models of P‐wave velocity (Vp) and P‐wave and S‐wave ratio (Vp/Vs) of the crust and upper mantle in the Eastern and eastern Southern Alps (northern Italy and southern Austria) were calculated using local earthquake tomography (LET). The data set includes high‐quality arrival times from well‐constrained hypocenters observed by the dense, temporary seismic networks of the AlpArray AASN and SWATH‐D. The resolution of the LET was checked by synthetic tests and analysis of the model resolution matrix. The small inter‐station spacing (average of ∼15 km within the SWATH‐D network) allowed us to image crustal structure at unprecedented resolution across a key part of the Alps. The derived P velocity model revealed a highly heterogeneous crustal structure in the target area. One of the main findings is that the lower crust is thickened, forming a bulge at 30–50 km depth just south of and beneath the Periadriatic Fault and the Tauern Window. This indicates that the lower crust decoupled both from its mantle substratum as well as from its upper crust. The Moho, taken to be the iso‐velocity contour of Vp = 7.25 km/s, agrees with the Moho depth from previous studies in the European and Adriatic forelands. It is shallower on the Adriatic side than on the European side. This is interpreted to indicate that the European Plate is subducted beneath the Adriatic Plate in the Eastern and eastern Southern Alps.

In the Central Alps, the lower crust and upper mantle of the European Plate were subducted southward beneath the Adriatic Plate (e.g., Pfiffner et al., 1997;Schmid et al., 1996). It has been proposed that the lower crust of the Adriatic Plate is detached from its upper crust and mantle and forms a narrow wedge within the Central Alpine orogenic crust, which overlies the European lower crust (Diehl et al., 2009;Kissling et al., 2006;Rosenberg & Kissling, 2013;Schmid et al., 1996). To the east of the Giudicarie Fault (GF), along the N-S oriented TRANSALP geophysical transect of the Eastern Alps at about 12° longitude, most studies indicate a south-dipping European Moho dipping beneath the shallower Adriatic Moho (Castellarin et al., 2006;Gebrande et al., 2006;Kummerow et al., 2004;Lüschen et al., 2006). Gebrande et al. (2006) suggested an Adriatic lower crustal wedge indenting the European orogenic crust, similar to what had been established for the Central Alps (see Castellarin et al., 2006 for an overview). Questions remain on the exact extent of this potential Adriatic lower crustal wedge beneath the Eastern Alps and how it relates to the one proposed in the Central Alps. East of the GF, teleseismic tomography of Lippitsch et al. (2003) indicated a positive Vp anomaly dipping NNE-ward to a depth varying from 150 km beneath the TW to 240 km further to the northeast at around 48°N/16°E. Together with the then-available models of Moho geometry (Waldhauser et al., 2002), this was interpreted as evidence for subduction of the Adriatic Plate beneath the European Plate. A first-order change in the subduction polarity somewhere beneath the transition from the Central to the Eastern Alps was thus proposed, which was initially supported by other seismological (Karousová et al., 2013;Kissling et al., 2006;Zhao et al., 2016) and geological studies (Handy et al., 2015;Kissling et al., 2006;Schmid et al., 2004). However, other studies, both before and since, do not favor a change in subduction polarity and postulate a single European slab beneath the Alps extending from west to east (Dando et al., 2011;Koulakov et al., 2009;Mitterbauer et al., 2011;Paffrath et al., 2021;Piromallo & Morelli, 2003). In recent interpretations, the part of the European slab beneath the Eastern Alps is largely detached from its orogenic lithosphere and dips southward in the Western Alps (∼9°E) and sub-vertically to steeply northward in the Eastern Alps (∼15°E; Paffrath et al., 2021).
Knowledge of the seismic velocity structure of the crust and upper mantle can improve tectonic models of the target area and provide information about the extent of the Adriatic lower crust and subduction polarity beneath the Central and Eastern Alps. The available seismic velocity models in the Alps are up to now either based on teleseismic tomography (e.g., Lippitsch et al., 2003;Mitterbauer et al., 2011;Paffrath et al., 2021) and thus cannot resolve the crustal structure or on the orogenic scale. Other models with local earthquakes as sources used less-dense seismic networks (e.g., Diehl et al., 2009;Solarino et al., 1997;Viganò et al., 2013) and thus have low resolution, hampering the detection of key crustal structures.
In this study, we present consistent 3-D Vp and Vp/Vs models of the Eastern and eastern Southern Alps derived from the inversion of a high-quality P and S travel-time data set based on Jozi . The high station density of the SWATH-D network (ca. 15 km spacing; Heit et al., 2017Heit et al., , 2021 allows for high-resolution images of the crustal structure in this part of the European Alps, even in a region of only moderate seismic activity. The seismic velocities (Vp and Vs) are influenced by the physical properties (e.g., density, porosity, fluid content, etc.) of the rocks. The Vp correlates positively with the density of the rock. Moreover, Vp and Vs are sensitive to the presence or absence of fluids, so that both Vp and Vp/Vs models provide information on the composition of the crust and mantle including possible partial melt.
The focus of our study is on the first-order features in the crust, especially on the extent of the lower crust beneath the Eastern Alps. Moreover, the crust-mantle discontinuity (Moho) in the target area is investigated in order to shed light on the subduction polarity beneath the Eastern Alps.

Local Earthquake Data
The waveform data used for this study was recorded by 161 stations of the temporary local SWATH-D network (Heit et al., 2017(Heit et al., , 2021  Jozi  analyzed a subset of this data (between September 2017 and the end of 2018) and generated a catalog of local earthquakes using the Markov chain Monte Carlo (McMC) approach, which uses a simultaneous inversion for hypocenters, station corrections, and the 1-D velocity structure (Ryberg & Haberland, 2019; see also ;Kissling, 1988). In the current study, we complemented the catalog of Jozi  by collecting more network data (from January to November 2019) and analyzing them using the same procedure. The final catalog consists of 499 local earthquakes which closely resembles the seismicity shown in Jozi their Figure 11).
The corresponding high-quality arrival-time data set was then utilized to perform the LET. The considered earthquakes have an azimuthal gap of less than 180° and at least 10 P and 5 S observations with a total of 16,484 P and 9,554 S picks (quality classes of 0-3; see Jozi  for more information). The average number of picks per event is 33 and 19 with average picking errors of 0.12 and 0.21 s for P and S observations, respectively (Table 1). More information on the number of P and S picks in each epicentral distance can be found in Figure A1.

Tomographic Inversion
We utilized the well-established SIMUL2000 algorithm (Eberhart-Phillips & Michael, 1993;Evans et al., 1994;Thurber, 1983), one of the most widely used approaches for LET. SIMUL2000 inverts the travel-time residuals simultaneously for 3-D velocity structure (Vp and Vp/Vs models) and earthquake parameters. The inversion is performed using an iterative damped least-squares technique to solve the linearized, coupled hypocenter-velocity problem. The tomography solutions strongly depend on the choice of initial parameters (i.e., velocities and hypocenters), model parametrization, and regularization (i.e., damping values).
The model is formed by nodes on a regular grid. The SIMUL2000 code requires a 3-D seismic grid with initial Vp and Vp/Vs values associated with the nodes. The velocity between grid nodes is linearly interpolated (Thurber, 1983). The code comprises the solution of the forward problem (for the travel-times) as well, which is accomplished using a combination of approximate ray tracing (Thurber, 1983) and a pseudo-bending scheme (Um & Thurber, 1987).

Initial 1-D Vp and Vp/Vs Models
For the initial Vp values at seismic grid nodes, we adopted the 1-D model by Jozi  for the crust which was calculated using the McMC approach for a subset of the data used here. Since our data set was retrieved from a set of local earthquakes predominantly located in the upper crust (depth <20 km) and the number of large epicentral distances is limited (due to the relatively short network time frame of 2 years), the deeper levels (particularly the upper mantle) are less well-resolved. Accordingly, we decided to use the 1-D velocity model of Diehl et al. (2009) for greater depth (i.e., depths below 45 km) which shows higher velocities, representative for the upper mantle. For the initial Vp/Vs model, constant homogeneous values of 1.72 for the crust and 1.76 for the upper mantle (retrieved from the Wadati diagram) were used as an initial estimate. The initial Vp and Vp/Vs models at depth layers of our selected grid are given in Table 2 and Figure 2.

Model Parametrization and Inversion Scheme
Model resolution can be strongly affected by the model parametrization which should thus be carefully selected. A coarse parameterization yields highly resolved model nodes due to a large number of rays passing through a block of the grid. However, it can potentially overlook small features (smaller than the grid node spacing) in the subsurface, whereas a too fine parameterization yields less well resolved nodes and gives sparse images of subsurface structures.
In the regions with variable ray sampling due to non-uniform earthquake and station distribution, some studies (Abers & Roecker, 1991;Bijwaard et al., 1998;Kissling, 1988;Thurber & Eberhart-Phillips, 1999) use uneven node spacing for enhancing ray sampling. These techniques give the possibility of using finer-scale parameterization where the data allows and improving the absolute model resolution by using a coarser grid where the ray coverage is sparser. However, the non-uniformity of the resolution within a given region could cause ambiguities in the interpretation of the velocity model, and velocity smearing between grid nodes might be interpreted as features.
Another strategy often used in regions with heterogeneous ray sampling is the graded inversion approach (Evans et al., 1994;Husen et al., 2000Husen et al., , 2003. In this approach, the inversion is first performed with a coarse parameterization and the obtained 3-D velocity model is then used as the initial model for another inversion with finer parameterization. This approach guarantees a smooth high-resolution velocity model throughout the volume  with dense earthquake and station distribution without introducing artifacts (Eberhart-Phillips, 1990;Haberland et al., 2009;Husen et al., 2003).
After several tests with various approaches and parameterizations, we found the most appropriate model (also based on synthetic tests) using a graded scheme. The inversion with a coarse grid of 30 × 30 km 2 (horizontally) was initialized using the 1-D Vp model as described above (Section 3.1 and Table 2). Subsequently, the output 3-D velocity model was linearly interpolated (Thurber, 1983) into a finer grid of 15 × 15 km 2 (horizontally) and used as input for the next finer grid inversion. To adapt the retrieved anomalies after the inversion of the coarse grid, we shifted the nodes of the fine grid 7.5 km toward east and north, so that the nodes of the coarse grid are located between the nodes of the fine grid (see Figure 3). However, taking into account the very dense ray coverage only in the central part of the region (Figure 3), the grid cells in the periphery of the network were kept the same as the coarse parameterization. In both coarse and fine grid inversions, the vertical layer spacing is 5 km in the upper crust, 10 km in the lower crust, and 15 km in the upper mantle (Table 2 and Figure 2).
As for example, pointed out by Diehl (2008), due to the trade-off between model parametrization, minimum resolvable velocity perturbation, and average data error, parametrization plays an important role in resolution. A velocity perturbation with amplitude ΔV (%) within the volume of a grid block with a background velocity of V0 can be resolved only if the average picking error is less than the travel-time residual between the background model and the perturbed model. Given the data error of 0.12 s and a volume with Vp = 6 km/s, a 5% velocity perturbation (a rough estimation of what is expected in the upper crust of our study region) needs a ray length of at least ∼15 km passing through the volume to be resolved (the reader is referred to the aforementioned publication for the mathematical formula). By assuming the ray length being almost equal to the volume dimension, our chosen parameterization of 15 × 15 × 10 km 3 can reliably resolve expected crustal structures with at least 5% velocity perturbations.

Regularization
The appropriate damping values of the inversions for Vp and Vp/Vs were determined by analyzing a series of single-iteration inversions and using a wide range of values. After inspection of the so-called trade-off curve (Eberhart-Phillips, 1986), the best compromise between data variance and model variance for Vp was obtained  with damping values of 500 for the coarse grid and 150 for the fine grid. These values were calculated by running the inversions with a fixed high Vp/Vs damping factor. Thereafter, in a similar procedure, the Vp/Vs damping factor was tested with the fixed chosen Vp damping and the final damping values of 700 and 500 produced the best compromise for coarse and fine grids, respectively.
The earthquake hypocenters were kept fixed in the first iteration to avoid any systematic shifts due to setting the station correction to zero. The earthquakes were relocated in the following iterations after the first one (Husen et al., 2003). After the relocation of the hypocenters using the 3-D velocity model, they shifted with a mean of 1 km and a standard deviation of ±1 km in the horizontal direction (mostly in the N-S direction). The values in the vertical direction are 0 and ±3 km, respectively. The inversion was terminated after a total of seven iterations for coarse and five iterations for fine grids. The final RMS misfit was reduced by 30% after the inversion with fine parameterization compared to the 1-D inversion . The final RMS misfit was 0.25 s which is still larger than the average picking error of 0.12 s. This discrepancy could indicate unresolved anomalies at a scale smaller than the node spacing which are not recoverable.

Resolution Assessment
The assessment of the solution quality is crucial to recognize the volumes of high, poor, or no resolution, smearing effects, and artifacts that might inadvertently cause improper interpretation. Only with the analysis of the solution quality, meaningful interpretations of the results can be made. Moreover, validation of the model parameterization and damping is essential.
The absolute resolution of a volume strongly depends on the density and geometry of the rays passing through that volume. Several parameters such as derivative weighted sum, hit count, resolution diagonal elements, and spread function (Michelini & McEvilly, 1991;Toomey & Foulger, 1989) are measures to assess the resolution. Here, we display the spread function values together with the smearing contours of the model resolution matrix (MRM; Eberhart-Phillips & Michael, 1998). Furthermore, we validate the chosen inversion strategy, parameterization, and damping using the classical checkerboard test as well as a synthetic test with a characteristic 3-D velocity model (expected real structures based on previous studies).

Spread Function
The spread function synoptically analyzes the information in each row of the MRM. Nodes with low spread values indicate good resolution, whereas the larger the spread value, the less well-resolved the particular node. High values of the spread function indicate smearing; however, it does not provide information on the direction of smearing. To show the directions, we plot the 70% smearing contour (following e.g., Diehl et al., 2009;Haberland et al., 2009;Haberland et al., 2014) for each node which is obtained by normalizing the corresponding row of the MRM by the value of the diagonal element. The well-resolved nodes with no smearing have circled contours around the center, whereas other nodes (mainly nodes with fewer crossing rays) tend to be smeared by neighboring nodes and show elongated contours (Eberhart-Phillips & Michael, 1998;Reyners et al., 1999;Toomey & Foulger, 1989).
Utilizing the spread value, we can distinguish the very well-resolved regions from those that are poorly resolved. However, since its range depends on the data set, parameterization, and regularization, there is no universally applicable threshold for the spread value (Toomey & Foulger, 1989).
Figures 4a and 4b show the spread values and the smearing contours of Vp and Vp/Vs, respectively. The Vp spread values from the surface down to 20 km depth are quite low and the well-resolved area is larger at the shallower layers. This resolution geometry is dominated by the seismicity distribution which is mostly restricted to the uppermost 20 km of the crust. Small spread values indicating excellent resolution are found in the volume with high seismicity. Reduced, nevertheless acceptable, resolution at a depth range of 30-45 km in the central volume of the region is indicated by moderate spread values. This is a consequence of diving Pn phases in these parts of the model. Below 50 km depth, no resolution is expected.
Likewise, the Vp/Vs model is well-resolved at the depth layers of 0, 5, and 10 km. At the 20 km depth layer, the Vp/Vs model is fairly resolved and below which the resolved area is significantly reduced. The Alpine area is one of the excellent examples, in which the direct Sg (Pg), the Moho-reflected SmS (PmP), and the Moho-refracted Sn (Pn) waves arrive closely spaced in time. However, our data set has only very few Sn arrivals (see e.g., Figure A1) mainly due to the short observational period of only 2 years.
The spread value of 2.7 was empirically defined as a reasonable threshold that fairly separates well-from poorly resolved volumes (also considering the recovery from synthetic tests; see below). The areas with spread values lower than this threshold illustrate minor smearing. However, in general, the E-W smearing is more considerable than N-S and vertical smearing, which we interpret as a consequence of the stations and earthquake geometry. Due to the short maximum distance between the stations and earthquakes in the N-S direction, the number of observed Pn arrivals in this direction is reduced. Therefore, the solution in the deeper parts is influenced mostly by the E-W Pn rays and consequently, the model is smeared in this direction.

Checkerboard Test
Checkerboard tests intuitively assess the solution quality and are quite common in tomographic inversions to analyze the blurring and smearing in the final model. Given the fine parameterization, we constructed a checkerboard model with the alternating anomalies spanning over two nodes in the x-direction and y-direction. Between the alternating anomalies, two nodes were left unperturbed to check the horizontal smearing. The alternating anomalies are placed at depths of 0, 10, 30, and 60 km with swapped polarities of the anomalies, and the other layers were left unperturbed to check the vertical smearing. The velocity perturbation includes ±10% of the 1-D model for Vp and ±5% for the Vp/Vs model ( Table 2). The synthetic travel times were computed using the forward calculation of the SIMUL2000 code and the same source-receiver distribution as the real data. Random noise was added to the synthetic travel times using the corresponding picking errors for each quality class (see Table 1). Similar to the real data inversion, we used the graded inversion approach and determined the damping values using the trade-off curves.
As can be inferred from Figure 5, the Vp anomalies are very well-resolved for the grid nodes with spread values less than 2.7. The high/low input anomalies are well recovered at depth layers of 0 and 10 km with minor vertical smearing at 5 and 20 km depth layers. The horizontal smearing is, however, very small down to a depth of 20 km. Below 20 km the vertical and horizontal smearing is higher, and the well-recovered area is smaller (similar to the resolved region indicated by small spread values; Figure 4). Below ∼50 km depth, no recovery of the anomalies in the checkerboard is expected. This was already inferred from the spread function (Figure 4) that the resolution below 50 km depth is low to non-existent. The geometry and amplitude of the Vp/Vs anomalies are well recovered down to 10 km depth with minor vertical smearing ( Figure 5) and fairly resolved at 20 km depth. Since the number of S arrivals from large epicentral distances in this data set is limited, in particular Sn arrivals, we do not expect a good recovery of the Vp/Vs below 20 km.
For further assessments, two additional checkerboard tests were conducted. In the first one, altering anomalies of Vp and Vp/Vs were placed at depths of 5, 20, and 45 km, while leaving the 0, 10, and 30 km layers unperturbed. This test had a similar result as the previous one ( Figure 5); that is, very good recovery of the main anomalies with minor smearing within the region with a spread value less than 2.7 and poor or no recovery outside ( Figure A2).
In the second test, the alternating high and low Vp anomalies were extended over one node with one unperturbed node between the alternating anomalies. The horizontal layers were placed at depths similar to the previous test ( Figure 5) and the forward and inverse problems were solved in the same manner. Using this test, we can check the recovery of the Vp anomalies in the size of one grid node spacing (15 km). Figure 6 shows that fine-scale Vp anomalies are geometrically very well recovered down to a depth of 20 km (where we have the most concentration of the crossing rays). However, the amplitudes are reduced by several percent, especially at deeper sections. In general, the bigger the size of the anomaly, the better its amplitude is recovered. Below 30 km, substantial smearing was observed. In a similar test for the Vp/Vs model, the geometry of the small-scale anomalies was resolved down to the depth of 10 km, however, the amplitudes were significantly reduced ( Figure A3).

Synthetic Test Using a Close-to-Realistic 3-D Velocity Model
A test attempting to reconstruct a synthetic Vp model is important because it allows for checking the potential of our data set in resolving the Moho topography. Our synthetic 3-D Vp model is aimed to mimic the realistic structure of the region similar to the synthetic Vp model of Jozi . It uses the Moho topography of the European and Adriatic Plates following Spada et al. (2013) without features such as the Moho overlap or gap (Figure 7a). The velocity gradually increases from the surface (6 km/s) down to the base of the crust, thereafter, it has a sharp velocity gradient from 7.0 to 8.0 km/s at the Moho depth. For more information on the initial velocity values, the reader is referred to Jozi .  The synthetic P and S arrival times were calculated for the earthquake and station geometry of the real data by using the 3-D FD Eikonal solver (Podvin & Lecomte, 1991;Tryggvason & Bergman, 2006). Synthetic noise was subsequently added to the travel-time data set (individually for each pick quality class; P and S picking uncertainty in Table 1). Ultimately, the synthetic travel times were inverted with the same initial velocities (as shown in Figure 2) and inversion strategy as in the case of real data (Section 3).
Figures 7c-7e show the recovered model after inversion through two N-S and one E-W cross-sections. We have to consider that the velocities in LET are interpolated between the nodes that are predefined by the model parameterization. Therefore, the sharp velocity contrasts (e.g., at the Moho) are broadened and can only be imaged with velocity gradients. It has proven to be useful to define a velocity contour line as a proxy for the velocity contrast expected at the Moho depth (see e.g., Diehl et al., 2009). In our synthetic test, a value of 7.7 km/s (dashed black line) in the inversion result showed the best fit to velocity step in the synthetic (input) model (dashed white line), however, the absolute depths depend on the absolute velocities and the initial model.
The recovered Moho proxy matches the Moho depth of the synthetic model very well in most parts of the resolved areas. The Moho (vertical) offsets are also recovered in terms of depth; however, they are significantly broadened and not sharply defined (due to the grid parameterization, interpolation, and smearing effects).
The low-velocity anomalies in the sedimentary basins of PoB, VFB, and MoB (see dashed red lines in Figure 7) cannot be observed in the recovered model. However, these sedimentary basins are located on the periphery of the network, and the model resolution (ray coverage) does not allow for thoroughly recovering of these anomalies.
No other crustal heterogeneities were included in the synthetic model. The recovered model is smooth without strong artifacts (or smearing) which provides further evidence for the high-resolution capability of our data set.

Results
As indicated by the resolution analysis (see above), the strength of our tomographic results lies in imaging the upper and lower crust. The 7.25 km/s iso-velocity line which we treat as a Moho proxy is fairly well recovered beneath the Giudicarie-Lessini region (GL), the TW, and the Inntal region. The Vp and Vp/Vs models as well as the seismicity pattern are shown in the horizontal slices in Figure 8 and vertical cross-sections in Figure 9. These sections illustrate that the velocity in the crust, in particular the upper crust, is both laterally and vertically highly heterogeneous. Several prominent anomalies in the crust can be identified: At the surface, reduced Vp values (down to 4.6 km/s) are observed in the resolved part of MoB and VFB and high Vp values (∼6.5 km/s) are seen in the Friuli-Venetia (FV) region down to a depth of 10 km (anomaly A; Figures 8a, 9c, and 9d). Further toward the east (along the Italian-Slovenian border), there are also elevated P velocities (6.5-6.9 km/s) found. This anomaly seems to be the continuation of anomaly A. However, its SE part extends down to 20 km. This region is also characterized by predominant high values of Vp/Vs (>1.85). In the upper crust, further increased Vp values (6-6.5 km/s) are observed northeast of GL and east of NGF at depth slices of 10 and 20 km (anomaly B ; Figures 8a, 9a, and 9f).
The TW region is characterized by a complex velocity structure featuring higher velocities (6-6.5 km/s) at the surface (anomaly C; Figure 8a At 30-km depth (Figure 8a), a prominent high Vp anomaly (>6.8 km/s; anomaly E) is observed in the central region of the study area. This anomaly extends from the junction of the SGF and NGF toward the eastern end of the TW (approximately at ∼13.5°E; depth slice of 30 km of Figure 8a). This anomaly is best expressed (thickest) south of and beneath the TW and beneath the PAF (see Figure 9). Figures 9b, 9c, 9d, and 9f reveal a low-velocity region (∼5.5 km/s) to the south of the PAF coinciding with the FV region (anomaly F). This implies that in this region, the velocity values corresponding to the upper crust (reddish colors) reach down to a depth of around 25 km, that is, the upper crust is thickened. This anomaly indicates a sharp lateral velocity transition beneath the PAF in the upper crust with higher velocity in the north of the PAF and lower values south of it. However, this transition disappears east of 13.5°E and west of 11.5°E (Figure 9f). Anomaly A with higher velocity (corresponding to the FV region) is located at depths of 0 and 5 km within the anomaly F.
The Vp/Vs model (Figure 8b) is resolved to a depth of only ∼10 km. This is mainly due to the limited number of Sn-picks, so that the deeper parts of the Vp/Vs model are not well resolved. The Vp/Vs values vary between 1.57 and 1.87 throughout the study area from the surface down to 10 km depth (Figure 8b). The most prominent feature in the Vp/Vs model is a zone of high values (>1.85) in the FV. Another zone of increased Vp/Vs values is observed in the NW of Engadine Fault (SE of MoB); however, this region is at the edge of the resolved volume. As seen in the depth slices of 0 and 5 km of Figure 8a, this region has lower Vp values (∼4.6 km/s).
LET velocity models are typically relatively smooth and sharp discontinuities (e.g., the Moho) are significantly broadened. Mantle velocities (>8 km/s) are observed at depths >45 km in our tomographic images, with an increased gradient between ∼30 and 45 km depth. Therefore, we use an iso-velocity contour of Vp = 7.25 km/s to define a proxy for the Moho in these gradual models (dashed black lines in sections of Figure 9; please see also Section 4.3 and Diehl et al., 2009). The Moho proxy is shallower on the Adriatic side (as shallow as 35 km depth) than on the European side (up to ∼55 km; see Figures 9a-9c).
The seismicity pattern is very similar to that of Jozi , with earthquakes concentrated mostly in the nappe stack south of the PAF. Seismic activity is observed in the Giudicarie-Lessini region, in the Engadine area, and in the Austroalpine nappes. The focal depths show that most of the earthquakes occurred down to 22 km depth within the crust. However, the recovered seismicity in our 3-D Vp and Vp/Vs models suggests three earthquakes at the base of the crust in the Giudicarie-Lessini region (depth slices of 30 and 45 km in Figure 8a).

Upper Crustal Structure
The low Vp anomaly in the MoB (Figure 8a) seems to extend to the south down to a depth of 5 km and is attributed to the marine and freshwater Molasse buried beneath the Helvetic and Penninic units Pomella et al., 2015;Zerlauth et al., 2014). The high Vp values of anomaly A in the FV region (Figures 8a, 9c, and 9d) confirm the observation of Diehl et al. (2009) andBrückl et al. (2007) and may be related to dense Mesozoic limestones and dolomitic rocks (Bressan et al., 2012). Velocity values between 6 and 6.8 km/s are expected for the limestones and dolomites (Castagna et al., 1985). Moreover, a prominent cluster of seismicity down to a depth of ∼17 km seems to be associated with this high Vp anomaly (Figures 9c and 9d), which is part of the eastern Southern Alpine deformation front. Sadeghi-Bagherabadi et al. (2021) also observed a similar feature in their Vs model and showed that a seismic cluster is confined to the S-wave velocity of ∼3.1 and ∼3.6 km/s down to about 20 km.
The Vp/Vs model in the FV region and external Dinarides (Figure 8b) indicates predominantly high values (>1.85) down to 10 km depth, which fits the values reported for dolomites and limestones (∼1.75 and ∼2.0, respectively; Castagna et al., 1985). However, the amplitude of our Vp/Vs anomaly becomes smaller at 10 km depth. The Vp/Vs values in the range of 1.79-1.87 confirm observations of Bressan et al. (2012) and are related to Cretaceous and Jurassic limestones.
The high P-wave velocity body east of NGF (anomaly B) coincides with Permian magmatic rocks in the eastern Southern Alps (purple line at depth slices of 10 and 20 km of Figure 8a; Tadiello & Braitenberg, 2021). This higher velocity anomaly may indicate the continuation of Permain magmatic rocks down to the base of the upper crust. A similar feature was observed in the study by Sadeghi-Bagherabadi et al. (2021) with higher Vs values.
The higher P velocities in the region of the TW (anomaly C) can be explained by the European basement nappes (Penninic units) that are exhumed and exposed at the surface (e.g., Rosenberg et al., 2018;Scharf et al., 2013;Schmid et al., 2013). This high-velocity anomaly is underlain by a lower velocity layer in some parts (anomaly D), although the Vp rises again at depth of around 15 km. However, resolution of this thin layer is limited by our model parameterization (vertical node spacing of 5 km.

Lower Crustal Structure
Anomaly E has seismic velocities typical of the lower crust (i.e., 6.8-7.25 km/s; e.g., Holbrook et al., 1992). The boundary between lower and intermediate crust is taken to be the 6.8 km/s iso-velocity line which is obtained by averaging Vp measurements of typical intermediate and lower crustal rocks from the Ivrea crustal cross-section exposed in the westernmost Southern Alps (Burke & Fountain, 1990). The lower crust defined this way is anomalously thick in the central part of our study area, a feature which we refer to as the lower crustal "bulge." Figure 10 focuses on this lower crustal bulge and its bearing on previously proposed structural models along the TRANSALP transect.
Two models, the so-called "Crocodile" and "Lateral Extrusion" models, based on active-source reflection seismology (Castellarin et al., 2006;Gebrande et al., 2006;Lüschen et al., 2004Lüschen et al., , 2006 as well as line drawings from Vibroseis and explosive seismic reflections of Gebrande et al. (2006) are shown in Figures 10a and 10b, respectively. Both models feature a wedge of Adriatic lower crust within the orogenic crust. However, the wedge is more pronounced in the "Lateral Extrusion" model. This model derives its name from the proposed connection at depth of conjugate, oblique-slip faults (the IN, SEMP, and PAF faults; see Figure 8a) that together accommodated eastward motion of the exhumed orogenic crust in the Tauern Window in the latest Oligocene to Miocene time (e.g., Favaro et al., 2017;Scharf et al., 2013).
Our tomographic images in the TRANSALP section indicate that the lower crustal bulge is located mainly beneath and south of the PAF and south of a prominent S-dipping thrust zone, the "Sub-Tauern Ramp (STR)." As interpreted in Figures 10a and 10b, the STR reaches from the Inntal Fault southward some 80 km . This large thrust zone is thought to be associated with the main collisional phase of N-directed thrusting their Figure 6). The age and origin of the lower crustal bulge are unknown; it could be made up of either European or Adriatic lower crust depending on how the faults beneath the Tauern Window are drawn. In Figure 10c, we interpret the tip of the bulge to truncate and therefore postdate the STR. In this scenario, the lower crust detached from its Adriatic upper mantle and protruded to the north. It was emplaced directly onto the downgoing European lower crust. The lower crustal bulge would thus form the northern limit of the Adriatic Indenter at depth. It is tempting to relate this protrusion or wedge to N-S-directed shortening, exhumation, and lateral extrusion of the Penninic units in the TW in Miocene time. Our model contrasts with both models of TRANSALP, which envision emplacement of Adriatic lower crust at shallower depth onto European intermediate crust (e.g., Castellarin et al., 2006;Gebrande et al., 2006). We emphasize that our interpretation does not represent a unique solution. Other interpretations involving both European and Adriatic lower crust are currently being tested by comparing the volume of the lower crustal bulge with estimates of shortening in the Eastern and Southern Alps.
In the Western and Central Alps along the NFP-20E transect, Kissling et al. (2006) imaged a much thinner lower crustal wedge than the bulge depicted in Figure 10c. It has been proposed that this wedge comprises the Adriatic lower crust that was inserted northwards during Oligo-Miocene collision beneath S-directed back-folded and back-thrusted upper crust (Schmid et al., 1996(Schmid et al., , 1997.

Crust-Mantle Discontinuity
The crust-mantle velocity discontinuity (Moho) is approximated by an iso-velocity contour at Vp = 7.25 km/s, similar to the Moho proxy adopted by Diehl et al. (2009). This yields a Moho depth that is consistent with that obtained from other methods (e.g., controlled-source seismology, receiver functions, local earthquake tomography (LET), and S-wave velocity; Diehl et al., 2009;Hetényi, Plomerová, et al., 2018;Kummerow et al., 2004;Sadeghi-Bagherabadi et al., 2021;Spada et al., 2013). Our Moho proxy (dashed black line in the cross-sections of Figure 9 and iso-velocity contour of Figure 11) changes depth south of the PAF, indicating a significantly shallower Adriatic Moho than European Moho. This is consistent with the southward subduction of the European lithosphere.
Other methods reveal an upward step or a gap in the Moho going from Europe in the north to Adria in the south.   Gebrande et al. (2006), respectively. The solid red lines in (a) and (b) indicate the evolutionary models of "Crocodile" and "Lateral Extrusion" (Castellarin et al., 2006;Gebrande et al., 2006;Lüschen et al., 2004Lüschen et al., , 2006, respectively. The solid red lines in (c) illustrate our alternative interpretation based on our LET model and the iso-velocity contour of 6.8 km/s as the upper boundary of the lower crust (see text), and also using the line drawings of Gebrande et al. (2006) , 2013). However, discrepancies are seen beneath the TW and PAF. Due to vertical and horizontal velocity gradients in LET (instead of a velocity contrast, e.g., in receiver functions), it is difficult to clearly differentiate the Adriatic and European Moho or recognize their current position at depth (also indicated with a question mark in Figure 10c). Therefore, other findings, for example, from seismic reflection and receiver functions are needed to provide further insight into the Adria-Europe collision.

Conclusion
The high-quality arrival-time data set from dense seismic networks in the Eastern and eastern Southern Alps and the well-constrained hypocenters  allowed us to calculate high-resolution 3-D models of Vp and Vp/Vs down to depths of 50 and 10 km, respectively. In the Vp model, the upper crust is clearly resolved. Moreover, the lower crust defining the front of the Adriatic Indenter in the Eastern Alps is well observed in the tomographic images. The reliability of the derived crust-mantle discontinuity was studied using a synthetic test with a close-to-realistic 3-D velocity model. The tested Moho topography was recovered without strong artifacts or smearing in the well-resolved region. It is indicated by an increased velocity gradient instead of a discontinuity.
A high-velocity anomaly (6.8-7.25 km/s), found at a depth of about 30-50 km from the NGF toward the eastern end of the TW and imaged in this detail for the first time, indicates thickening of the lower crust beneath and south of the TW and PAF. Our tomographic Vp model supports several models of indentation, including one that we propose here in which the lower crustal bulge represents the tip of the Adriatic Indenter and comprises mostly Adriatic lower crust. The lower crust detached from its mantle and upper crust and protruded to the north. We think that thickening and northward indentation of the lower crust may have triggered exhumation and eastward lateral escape of the previously accreted European nappes in the TW. Tectonic balancing of late Oligocene to recent shortening may yield further insight into the make-up of the lower crustal bulge (European, Adriatic, or both?) beneath the Eastern and eastern Southern Alps.
The Moho topography as indicated by an iso-velocity contour of Vp = 7.25 km/s (Moho-proxy), derived from our LET, confirms previous findings of a shallower Adriatic Moho (32-40 km) than the European Moho (>45 km). The change in the Moho depth between the two plates generally coincides with a Moho step or gap as imaged with other seismological methods. Our new model is consistent with south-directed subduction of the European lithosphere beneath the Adriatic lithosphere in the Eastern Alps. Figure 11. Moho depth based on the iso-velocity contour of Vp = 7.25 km/s. The unresolved region is masked with white color and the poorly resolved areas (with spread value smaller than the threshold of 2.7) are faded. Figure A2. Assessment of the solution quality of the Vp (a) and Vp/Vs (b) models using classical checkerboard test with alternating anomalies over four nodes of the fine parameterization. This test resembles the test in Figure 5, however, the altering anomalies of Vp and Vp/Vs were placed at depths of 5, 20, and 45 km, while leaving the 0, 10, and 30 km layers unperturbed. Please see the caption of Figure 5 for further features of the figure.