Resolving Seismic Anisotropy of the Lithosphere–Asthenosphere in the Central/Eastern Alps Beneath the SWATH-D Network

The Alpine orogeny is characterized by tectonic sequences of subduction and collision accompanied by break-off events and possibly preceded by a flip of subduction polarity. The tectonic evolution of the transition to the Eastern Alps has thus been under debate. The dense SWATH-D seismic network as a complementary experiment to the AlpArray seismic network provides unprecedented lateral resolution to address this ongoing discussion. We analyze the shear-wave splitting of this data set including stations of the AlpArray backbone in the region to obtain new insights into the deformation at depth from seismic anisotropy. Previous studies indicate two-layer anisotropy in the Eastern Alps. This is supported by the azimuthal pattern of the measured fast axis direction across all analyzed stations. However, the temporary character of the deployment requires a joint analysis of multiple stations to increase the number of events adding complementary information of the anisotropic properties of the mantle. We, therefore, perform a cluster analysis based on a correlation of energy tensors between all stations. The energy tensors are assembled from the remaining transverse energy after the trial correction of the splitting effect from two consecutive anisotropic layers. This leads to two main groups of different two-layer properties, separated approximately at 13°E. We identify a layer with a constant fast axis direction (measured clockwise with respect to north) of about 60° over the whole area, with a possible dip from west to east. The lower layer in the west shows N–S fast direction and the upper layer in the east shows a fast axis of about 115°. We propose two likely scenarios, both accompanied by a slab break-off in the eastern part. The continuous layer can either be interpreted as frozen-in anisotropy with a lithospheric origin or as an asthenospheric flow evading the retreat of the European slab that would precede the break-off event. In both scenarios, the upper layer in the east is a result of a flow through the gap formed in the slab break-off. The N–S direction can be interpreted as an asthenospheric flow driven by the retreating European slab but might also result from a deep-reaching fault-related anisotropy.


INTRODUCTION
The development of seismic anisotropy is significantly affected by tectonic processes and the dynamics of the Earth's interior (Savage 1999). The lithosphere and the asthenosphere are the dominant source regions of seismic anisotropy in the crust and upper mantle. In the crust, anisotropy results from the stressinduced alignment of cracks (Christensen, 1966;Nur and Simmons, 1969;Nur, 1971;Crampin 1987;Yousef and Angus 2016) in alternating layers of different seismic velocities like folded sediments (Backus 1962;Savage 1999) or (in the lower crust) an alignment of intrinsically anisotropic minerals. The latter also occurs in the upper mantle lithosphere as frozen-in anisotropy. In the asthenosphere, seismic anisotropy is dominated by the lattice-preferred orientation of olivine crystals that align along the strain direction (Silver 1996;Karato et al., 2008), in response to mantle flow processes (Long and Becker 2010). The most robust observable for measuring seismic anisotropy is the splitting of shear waves, which is commonly inferred from core-mantle converted phases (Silver and Chan 1991;Savage 1999). Shear-wave splitting characterizes anisotropy by a fast axis, which is parallel to the preferred orientation of the mantle minerals in a simple mantle flow-dominated anisotropic model, and a lag time, which scales with the strength of anisotropy and the extent of the anisotropic volume. However, the diverse tectonic history and complex present-day dynamics influence the anisotropy of the upper mantle and lithosphere. Also, large or local shear faults may produce an extended crustal and lithospheric anisotropy above an asthenospheric mantle flow (e.g., Savage 1999;Reiss et al., 2016;Reiss et al., 2018). Lithospheric fragments of subducted slabs can carry frozen-in anisotropy. They result in a pronounced layered anisotropy when they sink (after possible break-off) into the mantle, due to a replacement flow of asthenospheric mantle material (Qorbani et al., 2015a).
The Alpine orogeny is a result of an ensemble of subduction, extension, and collision sequences driven by the convergence of the African and Eurasian plates (Dewey et al., 1973;Dewey et al., 1989;Handy et al., 2010). Prior to the collision of the Adriatic-African and European plates at around 35 Ma, Alpine Tethys and the distal European margin subducted during the convergence initiated at 84 Ma with Adria as the upper plate (Schmid et al., 2004;Handy et al., 2010). With the collision, the Adriatic plate was pushed by the African plate deep into Europe, which led to the formation of a complicated subduction geometry (Handy et al., 2010;Király et al., 2018). Concurrently, a slab break-off was initiated after the completed subduction of the Tethyan units and separated the crustal European units from the subducted lithosphere below the Eastern Alps (von Blanckenburg and Davies, 1995;Lippitsch, 2003;Harangi et al., 2006;Kissling et al., 2006). This was possibly accompanied or preceded by a subduction polarity switch (Schmid et al., 2004;Vignaroli et al., 2008;Handy et al., 2010;Molli and Malavieille 2011;Handy et al., 2015). Especially the subduction polarity below the Eastern Alps is under debate (Kästle et al., 2020), as tomographic studies indicate a northward dip of a possibly Adriatic subduction below Europe (Schmid et al., 2004;Ustaszewski et al., 2008;Schmid et al., 2013), which might be a result of a polarity switch accompanied by a break-off of the previously subducted European slab (Handy et al., 2010;Schmid et al., 2013). For further detailed descriptions of possible models, we refer to the study by Kästle et al. (2020). The present-day tectonic structure can be identified as accreted European units in the Central and Eastern Alps, partially showing metamorphism from deep subduction (Froitzheim et al., 1996;Handy et al., 2010, see also Figure 1). This is separated from the Southern Alps, consisting of Adriatic units by the Tonale Fault (TF) and the Periadriatic Fault (PAF). The main irregularity of the orogen parallel E-W trending faults arises by the Giudicarie Fault (GF), which formed during the indentation of the Adriatic microplate (Handy et al., 2015).
To address the question of subduction polarity and tectonic evolution in the transition of the Central to the Eastern Alps, the recently established seismic network SWATH-D provides a unique opportunity due to its unprecedented dense station coverage (AlpArray Seismic Network, 2015;Heit et al., 2017;Hetényi et al., 2018a). The 154 stations cover the critical transition between the TRANSALP (Lüschen et al., 2004) and the EASI experiment (Hetényi et al., 2018b) along the Periadriatic Fault (Handy et al., 2010). While the Moho geometry in the TRANSALP experiment was still unclear (Lüschen et al., 2004), the EASI experiment seems to image the Adriatic Moho below the European plate, which indicates an Adriatic subduction (Hetényi et al., 2018b;Kästle et al., 2020). While tomographic studies lack evidence for the subduction polarity in the Eastern Alps (Kästle et al., 2020 and references therein), they clearly suggest a European subduction in the Western and in the Central Alps (Lippitsch, 2003;Zhao et al., 2016). The result from EASI therefore indicates a polarity switch, possibly in the area covered by the SWATH-D complementary network. The present study addresses open questions regarding the slab break-off, slab geometry, and tectonic evolution in this key area of the Alpine system by analyzing the splitting of core-mantle phases due to seismic anisotropy beneath the network. Based on the findings from previous studies and our individual splitting measurements, we specifically account for two-layer anisotropic models which can account for both frozen-in lithospheric and asthenospheric mantle flows, or a present-day two-layer flow in the asthenosphere.

Individual Shear-Wave Splitting Measurements
When a shear wave propagates through an anisotropic medium, the wave is split into a fast and a slow phase polarized perpendicular to each other (Savage 1999). The polarization directions coincide with the horizontal projection of the fast and slow directions of the anisotropic medium. Usually, core-mantle converted phases (e.g., SKS, SKKS, PKS, and PKKS) are used for the analysis, as their initial polarization is determined by the event-station geometry, when they leave the Earth's core. Assuming a radially symmetric isotropic Earth, no splitting occurs, and the initial polarization in the radial direction is preserved in the recording at the receiver. In contrast, after the propagation through an anisotropic medium, seismic energy may be present on the transverse component due to shear-wave splitting. This anisotropic effect can be investigated by applying an inverse splitting operator with the aim to remove the splitting from the recorded data (Silver and Chan 1991). Here, we use the SplitRacer toolbox (Reiss and Rümpker 2017) to analyze the splitting of the core-mantle converted phases, which also includes standard procedures for requesting and processing large amounts of data from (virtual) seismic networks based on FDSN web services.

Assembling Data and Processing
The SWATH-D (Heit et al., 2017) network provides 2 years of continuous seismic data at 154 stations between 2017 and 2019. This data set is accessible from the GEOFON archive. We include data from permanent stations in the area for longitudes between 10 and 14.5°E and for latitudes between 45.5 and 47.5°N. The data are assembled from the networks BW, CH, IV, NI, OE, OX, SI, and SL (see Figure 1 and Supplementary Table S1 in the supplements), which are also part of the AlpArray backbone (AlpArray Seismic Network, 2015;Hetényi et al., 2018a). Teleseismic events with magnitudes above 5.8 and within the distance range from 89 to 140°are selected for download. Data with gaps are discarded in the first quality check. The traces are cut into 100-s windows centered around the arrival of an expected core-mantle converted phase. A signal-to-noise ratio is calculated based on the normalized signal energy in a 25-s window subsequent to the expected arrival time of the phase and the normalized noise energy in a 20-s window preceding this arrival time. Phases with a signal-to-noise ratio below 2.5 are discarded. The traces are band-pass filtered with corner frequencies of 0.02 and 0.25 Hz for all stations. The windows containing the phases for later analysis are then selected by visual inspection. The long period signal of the selected phases is used to determine their initial polarization (Rümpker and Silver 1998), which allows us to identify a possible station misorientation. We correct for the mean misorientation and consider temporal changes at affected stations (see also Supplementary Table S2). The phases are individually analyzed for splitting in a grid search for a pair of splitting time and fast axis direction, minimizing the remaining transverse energy after applying the inverse splitting operator (Silver and Chan 1991). The fast axis direction is measured clockwise with respect to north in this study. The results are categorized after manual inspection as good, average, and null according to the quality of the measurement. Poor measurements are discarded and not considered for further analysis. The quality is assessed based on the ellipticity of the waveform, noise contamination, energy reduction, and error of the results. Phases with linear particle motion and low noise contamination are considered as null measurements. If phases show an elliptic particle motion and a considerable improvement in energy reduction on the transverse component after applying the inverse splitting operator, the measurements are considered as good for low noise contamination and average for moderate noise contamination. If the error exceeds realistic measures, the phases are discarded as poor measurements. Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 679887

Analysis for Two-Layer Anisotropy at Station Clusters
Previous studies have found evidence for a two-layer anisotropy in the eastern part of the region covered by the SWATH-D network (Qorbani et al., 2015a). We therefore search for two anisotropic layers in our analysis. Instead of fitting calculated effective (or apparent) splitting parameters to the observed results obtained from an event-by-event splitting analysis, we choose to perform the more stable two-layer joint-splitting approach (Silver and Savage 1994;Rümpker and Silver 1998;Homuth et al., 2016). This also allows us to include null measurements and serves to improve the statistical significance of the results. In the two-layer joint-splitting approach, an inverse splitting operator is calculated for a set of four parameters (Silver and Savage 1994;Rümpker and Silver 1998), a splitting time, δt, and a fast axis, ϕ, characterizing the anisotropic effect of each layer. By applying the inverse splitting operator to an XKS waveform, the splitting of the two-layer model is removed from this phase. How well the existing splitting is removed is expressed by the remaining transverse energy. The inverse operator is applied for a range of splitting parameters, the delay time, and fast axis direction, for each of the two layers. The remaining energy is stored, which results in a four-dimensional matrix of the remaining transverse energy, which is directly related to the misfit function of the two-layer anisotropic model. For a single phase, a search for the minimum will result in highly ambiguous solutions. We therefore include (or join) all available phases by summing up the four-dimensional misfit function to search for the overall best fitting solution, which ideally minimizes the transverse energy of all phases simultaneously. This technique has produced more stable results in previous studies (Homuth et al., 2016;Reiss and Rümpker 2017) and should produce equivalent results to those of the waveform-based technique of Menke and Levin (2003). The joint-splitting technique is much less computationally expensive and less dependent on model assumptions than the modeling technique of Yuan et al. (2008), who, in return, were able to search for a plunge of the fast axis. However, a detailed comparison of the different techniques available is beyond the scope of this article. For our data set, 2 years of data are not sufficient for reliably determining two-layer splitting parameters at a single station. Therefore, we introduce a modification of the two-layer jointsplitting approach provided in the SplitRacer software. In view of the dense station spacing of the SWATH-D network and by considering the extent of typical XKS Fresnel zones in the upper mantle, we can assume that core-mantle converted phases largely sample the same anisotropic regime for neighboring stations. Therefore, we aim for the simultaneous analysis of neighboring stations (groups), with similar anisotropic patterns. As a measure of similarity between two data sets, we use their correlation coefficients (Everitt et al., 2011). We found that the four-dimensional distribution of transverse energy allows sufficient characterization of anisotropic features and waveform variability beneath a station. Therefore, we choose this as the basis for our measurement of correlation. The following describes our stepwise analysis procedure: • We first apply the regular joint-splitting analysis for two anisotropic layers to every station j and store the individual contribution, T ij (x), of every phase i to the overall misfit function, TE j (x), separately. The variable x depicts the parameter vector (δt 1 , ϕ 1 , δt 2 , ϕ 2 ), where the indices 1 and 2 depict the lower and upper layers, respectively. • In a second step, we cross-correlate the four-dimensional energy tensor (or misfit function) of each station to produce a correlation matrix, c jk Here, we only consider stations that have at least one phase categorized as good (see also Figures 2A,B). This provides a measure of similarity of the energy tensor from each station and, therefore, similar anisotropic structures sampled at the stations. • Based on this assumption, a cluster analysis using a complete linkage is performed (Everitt et al., 2011). The hierarchical clustering is a stepwise approach. Data points are merged into one group, based on the distance to each other (measured using the correlation coefficient). The new distances of all data points to the merged set are then based on the maximum distance (or the maximum value of the dissimilarity matrix) from them to one of the points in the set. In the hierarchical clustering, this procedure can be done from all data points treated as single individuals to one group, which contains all data points (for more details, see Everitt et al., 2011). The ideal number of groups has to be chosen by the analyst and depends on the data set and the purpose of the investigation. The correlation matrix can then be sorted for the groups, to visualize the similarity of the clusters, which produces patches of high values around the diagonal elements of the matrix for similar groups (see also Figure 2C). • For each group found in the cluster analysis, a joint-splitting analysis of all phases measured at the dedicated stations is performed. Thus, four splitting parameters are identified, which minimize the overall transverse energy of all phases in this group. We base the joint splitting on a bootstrap approach with 1,000 assembled subsets (Efron 1979), to test the stability of the result. Each subset is sampled randomly from the set of phases in each group with replacement preserving the size of each set. The joint analysis for the subset is found for the stack of all individual energy tensors previously stored in the joint analysis of the stations. We discard results where both layers show a similar fast axis direction (ϕ ± 20°) or a phase shift of 90°[(ϕ + 90°) ± 20°], as this would be an indication for a single-layer anisotropy. If no two-layer anisotropic model can be found, a single-layer joint-splitting analysis is calculated. We determine the final result from the mean and the standard deviation of a Gaussian curve fitting the distribution from the bootstrapping.

Individual Measurements
The splitting analysis of individual events provides, in total, 3,123 pairs of splitting parameters. 792 of them are categorized as good, 1,467 as average, and 863 as null measurement. The splitting time for good and average measurements varies mainly between 1.3 ± 0.5 s and shows no major lateral variation. The dominant fast axis direction is between 50°and 65°, although there is a trend toward larger angles of 90°and up to 145°, which appear mainly with growing longitude. At first sight, the individual results can be interpreted as gentle clockwise rotation of anisotropic fast axis orientation from west to east (see Figure 3). However, the permanent stations in the east show a large variation of the individual results, which can also be seen to be less pronounced in the west. We further investigate the azimuthal distribution of the good measurements with respect to their station position in longitude (see Figure 4). The splitting measurements change laterally, mainly with longitude. The lateral and azimuthal change of the splitting parameters can be visualized in one representation, if we consider the longitude position of the station, where the measurements are recorded, as the radius in a polar plot, where the polar angle represents the back azimuth of the incoming event (see Figure 4A). While the measurements have the same lateral pattern for events coming from the opposite direction (see 45°-90°and 225°-270°, respectively), their fast axis orientation changes for events coming from different angles. This is clear evidence for the events sampling the same anisotropic origin, whereas the azimuthal variation points toward a possible anisotropic layering. Stations with a longitude greater than 13°E show the same fast axis direction for events coming from 180°to 225°and for events coming from 270°to 315°, indicating a 90°periodicity. This favors a multilayer anisotropic origin in the east. In a more classical representation of the splitting time and fast axis direction (see Figure 4B), the azimuthal variation for the stations located east of 13°E becomes evident. We estimate the two-layer characteristics of the single measurements by fitting synthetic apparent splitting parameters for a two-layer anisotropic model to our measured parameters. The synthetic splitting time and fast axis direction are calculated using the expressions introduced by Silver and Savage (1994). The apparent parameters are calculated for a two-layer model, where the splitting time of each layer is varied for 0.2 s between 0 and 4 s. The fast axis direction of each layer is varied in 10°steps between 0 and 180°. One synthetic pair of parameters for each measured pair is calculated considering the initial polarization defined by the back azimuth of the event. This allows us to define a misfit between the measured and synthetic parameters and leads to a favored model with minimum misfit. We find the best fitting model, which is defined by two equally strong anisotropic layers, an upper layer with a splitting time of 1 s and a fast axis direction of 120°and a lower layer defined by a splitting time of 1.6 s and a fast axis direction of 70°. This agrees well with the azimuthal pattern observed by Qorbani et al. (2015a) and their anisotropic model. However, while the central part of the network between 12°E and 13°E shows no significant change in the fast axis direction, a similar periodicity of weak azimuthal variation can be seen for stations west of 12°E. We investigate this further by analyzing the azimuthal variation for stations in the longitude range between 11 and 12°E. Here, we observe again a clear azimuthal variation in the splitting-parameter plots (see Figure 4C). The corresponding best fitting two-layer anisotropic model shows one dominating upper layer with a splitting time of 1.4 s and a fast axis direction of 50°. The lower layer is weaker but still considerably strong with a splitting time of 0.6 s. The fast axis direction of the deeper layer points toward 170°. This indicates a layered anisotropy for the western part of the network also, which will be further tested in a joint analysis of multiple events for two anisotropic layers in the following section. Based on the single measurements, we would expect a much smaller two-layer effect with one dominating layer (possible locally confined) in the west, compared to the east with the very pronounced variation in the anisotropic parameters.

Joint Analysis of Multiple Events for Two-Layer Anisotropy
In the following section, we test the entire data set for layered anisotropy, by searching for splitting parameters of two anisotropic layers which best fit the data. We include all stations, as this analysis also allows us to test for one dominating layer with the other layer showing only minor splitting. Therefore, no previous selection on stations to be used for this analysis is necessary. A cluster analysis is performed, where we search for stations with similarities in their four-dimensional matrix of the remaining transverse energy derived from a two-layer joint-splitting analysis, as described above. We discard stations without any good measurements for this analysis to avoid a distortion caused by noise-contaminated stations. From the complete linkage using the correlation matrix of station-wise correlated fourdimensional misfit functions, we can visually identify two main groups of similar energy tensors (see Figures 2A,C), which are separated at 12.5°E-13°E. However, we allow for additional variation within these two groups by introducing six different groups in total to consider minor lateral variation in the (one-or) two-layer anisotropic pattern (compare Figure 2A). The results from the joint splitting of the six groups all show robust splitting results ( Figure 5). Three groups (1-3, see Figures 5A-C) are characterized by an upper layer with a fast axis between 55 and 70°and a lower layer with a fast axis between 1 and 37°. The splitting time of these groups (1-3) show a dominance of the upper layer with a value of 0.8-1.2 s. Only group 1 (see Figure 5A) shows a clear, pronounced second layer with a splitting time of 0.5 s and a small error for splitting time and fast axis direction in both layers. Group 2 ( Figure 5B) is characterized by a considerable splitting of 0.4 s, but the error from the bootstrapping statistics of ±0.3 s shows the minor confidence in this parameter. Therefore, a one-layer solution might explain the data equally well. The splitting time for group 3 ( Figure 5C) itself is very weak with only 0.2 s and is therefore also explained equally well using a one-layer model. The larger error in the fast axis direction of the lower layer is also indicative of the minor relevance of the two-layer anisotropy. For the remaining three groups (4-6, see Figures 5D-F), the lower layer shows strong splitting with a splitting time of 0.8-1 s. The fast axis direction for the lower layer shows values between 65°and 73°. While group 4 ( Figure 5D) exhibits an upper layer with a splitting time of only 0.4 s, the other two groups (5 and 6, Figures 5E,F) show an equally pronounced splitting as the lower layer with 0.9 s. All three groups show a similar direction for the fast axis of the upper layer (106°-114°). Group 4 ( Figure 5D) also shows a larger error for the fast axis in the weaker upper layer, which is indicative of a smaller significance of the two-layer anisotropy, similar to group 2 ( Figure 5B). The two main groups (1-3 and 4-6) clearly separate into a western and an eastern subdomain at a longitude between 12.5°E and 13°E ( Figure 6). This coincides with the location of the EASI profile (AlpArray Seismic Network, 2014) and also with the location between 12 and 12.5°E, where Qorbani et al. (2015a) found a transition of anisotropic effects in this region. The western part is (at least locally) characterized by a lower layer of almost N-S aligned fast axis direction (unit D in Figure 6C) and an upper layer of approximately 60°fast axis direction consistently measured for the whole western area (unit C in Figure 6C). The eastern part is characterized by a lower layer with the same direction as the previous upper layer, a fast axis direction of 60°(unit B in Figure 6C). The fast axis in the upper layer of the eastern region is aligned roughly along 110°(unit A in Figure 6C). The larger error of the fast axis direction associated with group 4 may indicate that the error is introduced by the transition from the eastern to the western area, as parts of the rays coming from the west might intersect a larger portion of the eastern regime. This lateral change is not considered in the rather simplified two-layer anisotropy, where only vertical variation is assumed. The confidence in the two-layer anisotropy for group 4 ( Figure 5D) is supported by the small error in the splitting time of only 0.1 s.

DISCUSSION
We find a lateral and a vertical variation in the anisotropic properties of the Eastern Alps, similar to Qorbani et al. (2015a). However, the statistical approach based on bootstrapping in combination with the joint-splitting analysis also allows us to identify a two-layer anisotropy for the area west of 13°E, which was previously thought to be characterized by only a single anisotropic layer (Barruol et al., 2011;Bokelmann et al., 2013;Qorbani et al., 2015a;Petrescu et al., 2020a). The dense station spacing also allows a high lateral resolution of the anisotropic pattern, which supports the abrupt change between the eastern and western parts with a distinct boundary at 12.5°E-13°E. The two-layer joint splitting allows us to separate the study area into a western and an eastern subdomain. Both domains show a layer with a fast axis of ∼60°, defining the upper layer in the west and the lower layer in the east. While the lower layer in the western subdomain shows only weak anisotropy and N-S-direction, the upper layer in the eastern subdomain is equally pronounced as the lower layer and shows a fast axis of 100°-115°. When observing a possible vertical variation of seismic anisotropy, the crust is likely contributing at least part of this effect. Qorbani et al. (2015b) found wide similarities of single anisotropic directions to kinematic data for lower crustal material in the Tauern Window and alignment with fault structures. This is interpreted as vertical coherent deformation of the crust and the lithospheric mantle, mainly controlled by the indentation of the Adriatic microplate. Our study area is further extended to the south and covers a larger longitudinal range, which includes the main Alpine suture marked by the Tonale Fault (TF) (see also Figure 1 and Figure 7) and the Periadriatic Fault (PAF). The Giudicarie Fault (GF) differs from the mainly E-W aligned fault systems and forms a displacement of the two major faults. To investigate a possible alignment with deepreaching surface geological features, we back-projected the individual measurements into a depth of 120 km (see Figure 7), which is the maximum base of the lithosphere for the area of the Tauern Window (Jones et al., 2010;Bianchi and Bokelmann, 2014;Qorbani et al., 2015b). While Qorbani et al. (2015b) found a very good agreement of the fast axis direction and kinematic data with the faults in the small area around the Tauern Window [in particular with the Mölltal Fault (MF)], we see such an alignment only for a limited number or only part of the faults. At the Mölltal Fault (MF), only measurements that are very close to the fault appear to align with the surface structure. They extend to the north to a region of more variable fast axis direction. A large depth extent of the fault and therefore a considerable effect is not to be expected from this very localized apparent alignment. The measurements around the Lavanttal Fault (LF) system might also be interpreted as an aligned pattern of the fast axis with the fault. However, the measurements appear to have a rotational pattern around the station location, which coincides well with the pattern of the fast axis direction for the neighboring station in the S-W and the north. This favors a deeper origin (e.g., multilayer anisotropy) as for these other stations, no alignment with surface features is visible and the rotational pattern indicates a 90°periodicity. It is remarkable that the Periadriatic Fault (PAF) and the Tonale Fault (TF), as major sutures for the European-Adriatic collision, show no significant correlation with the fast axis directions. Only the Giudicarie Fault (GF) with the northern and southern parts seems to produce a more significant imprint on the fast axis direction. The single measurements align with the fault and also show (to its west) an irregular NNE-SSW alignment along the fault, which is not visible for the measurements in the vicinity. In summary, we only see the Giudicarie Fault (GF) as a major surface feature, which seems to connect to larger depths, producing a possible anisotropic effect in the shear-wave splitting. The remaining faults seem to produce either only a very locally confined effect or produce, surprisingly, no effect at all, for example, the Periadriatic Fault (PAF) and the Tonale Fault (TF). Seismic tomography can provide insight into slab geometries and lithospheric thinning as well as asthenospheric upwelling in terms of positive and negative velocity perturbations. However, so far, studies based on body waves show only a limited agreement of velocity anomalies with flow patterns derived from shear-wave splitting (see Qorbani et al., 2015a and references therein;Lippitsch, 2003;Koulakov et al., 2009). A recent surface wave tomography by Kästle et al. (2018), Kästle et al. (2020) provides significantly improved resolution for relatively shallow structures (between 50 and 200 km). The latter shows a distinct weakening in the positive velocity perturbation of the European slab at a depth of ∼120 km at the transition to the Eastern Alps. The start of the weakening of the positive anomaly is located at a longitude FIGURE 7 | Single splitting measurements with quality good and average projected into a depth of 120 km (gray bars). For better visualization, the measurements are averaged on a regular grid, which are shown as black bars. Same tectonic map as in Figure 1.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 679887 of ∼13°E and therefore coincides well with the changing anisotropic pattern found in our study. This weakening can also be identified in images from body wave tomography (Zhao et al., 2016 ;Kästle et al., 2020). The weakening of the positive velocity anomaly has been interpreted as the discontinuous European slab and marks a break-off at depths of about 120-200 km. If we assume 4% anisotropy and a mantle shear velocity of 4.5 km/s, the required thickness of the anisotropic volume is about 90-100 km to satisfy the splitting times found for the upper layer in the east. The depth range between the lithospheric base and the slab fragment is found to be around 80 km and agrees with the required volume (within errors).
The dominant anisotropy (of about 60°) in the western part coincides well with the Alpine orogeny and is aligned parallel-tosubparallel to the strike of the mountain belt. This is evidence for a vertical coherent deformation with coupled anisotropy from the lithosphere to the upper asthenosphere as proposed by Silver (1996). The fast axis is not perfectly aligned with the mountain axis but shows a constant trade-off of about 10-20°. While this FIGURE 8 | Schematic model of the lithospheric/asthenospheric interaction in the study region. A to D correspond to the anisotropic units shown in Figure 6 (A-B) European plate subducts below the Adriatic plate. During the subduction process, the European plate is pushed back and produces a weak N-S directed plate motion-induced flow at the bottom of the slab (D). At around 35 Ma, a slab break-off initiated, which opened a gap for asthenospheric flow at a depth between ∼120 and ∼200 km. (A) Asthenosphere produces a toroidal flow around the retreating European slab (B and C). The opening of the break-off causes a reorientation of the minerals parallel to the orogen C into the new flow direction A of 115°N. (B) As a result of the compressive regime, anisotropy develops parallel to the orogeny (B and C) within the crust and lithosphere due to vertical coherent deformation (VCD). This anisotropy is stored in the subducting European plate. The opening of the slab break-off results in asthenospheric flow of 115°N (A). While the compressional regime in the west remains constant with the VCD causing the measured upper-layer anisotropy and plate motion-induced flow causing the lower-layer anisotropy, the change of deformation in the eastern part is perturbed by the break-off event. Here, the upper-layer anisotropy is caused by the asthenospheric flow above the detached slab, which preserves the VCD of the former collision causing the lower-layer (frozen-in) anisotropy B.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 679887 disagrees with an ideal coherent deformation of the crust and the mantle lithosphere, the overall regional changes of anisotropy in the Alpine area clearly rejects the simple asthenospheric flow model. In that case, the fast polarizations are expected to align more-or-less parallel to the absolute plate motion. The high variability within the plate therefore favors a shallower origin of the anisotropy. The deviation of the orogen-parallel alignment may be indicative of fossil anisotropy that is mainly preserved in the deeper crust in combination with the mantle lithosphere related to a previous stage of the European-Adriatic collision. The fast axis direction would then point to the alignment of the pre-collisional front during the subduction of the Piemont-Valais oceanic units (Handy et al., 2015). An alternative scenario for the origin of the observed pattern in the western part of the study area is a slab-controlled flow that is mainly driven by the Adriatic intender causing a steep European subduction. The blocked mantle flow is then forced to flow around the retreating European slab, causing the alignment parallel to the orogeny (see Figure 8A, and Király et al., 2018;Petrescu et al., 2020a). Both scenarios seem equally likely for the western region, where the dominant upper layer exhibits fast axes oriented subparallel to the orogenic belt. However, our observations for the eastern subdomain show discrepancies with a slab-controlled orogen-parallel flow behind the Alpine Front; as part of the asthenospheric domain behind the slab is reoriented in the direction of the upper layer (representing the flow through the gap produced by the breakoff), we would expect to observe a weakening in the slabcontrolled splitting. As a consequence, the volume pushed back by the retreating slab is shortened, leading to a weaker anisotropic effect parallel to the orogeny. However, the measured direction and the strength of the anisotropy are preserved in the lower layer of the eastern subdomain and resemble the anisotropy of the upper layer in the west. Therefore, our two-layer splitting results support a lithospheric origin and the vertical coherent deformation model with orogen-parallel anisotropy or fossil anisotropy parallel to the pre-collisional subduction front (see Figure 8B). The upper layer in the east can be interpreted as an asthenospheric flow that is controlled by the break-off event, which initiated 35 Ma ago (Handy et al., 2015) when the subduction of the Tethys oceanic plate was completed. The contact of the Adriatic plate with the more rigid European continental plate favored the break-off and changed the dynamic condition in the collision. This led to a perturbation of the coherent deformation with a cancellation of the anisotropic effect in the upper lithosphere. The orogeny-parallel direction remains preserved in the detached slab fragment, which we measure as lower-layer anisotropy in the east. It may be noted that splitting observations further to the east in the Pannonian Basin and the Carpathians exhibit a similar fast axis direction to the upper layer in the eastern subdomain (Petrescu et al., 2020b and references therein). This is evidence for a larger-scale feature that is connected to the upper-layer anisotropy found in the eastern subdomain of our study area.
In contrast, the N-S directed anisotropy, which we resolve in the western subdomain, is quite weak and is only locally a stable feature, while the fast axis directions of all groups in the area point in a similar direction. The generally small delay times indicate less significant anisotropy or a non-horizontal fast axis, which may also be present across a wider region, but is suppressed, due to one or more dominating anisotropic layers. If we interpret the N-S directed anisotropy in the east as an asthenospheric flow below the lithosphere, such large-scale distribution would be necessary. The asthenospheric flow would then represent the return flow of mantle material following the dip of the subducting plate and submerge below the retreating European slab that is pushed by the Adriatic intender (see Figure 8 and compare with Funiciello et al., 2006 for laboratory experiments). The close vicinity of the group with clear two-layer results in the west to the Giudicarie Fault might favor this geologic feature as the origin of an apparent two-layer effect. As the back-projected single measurements also show an alignment with the fault (see Figure 7), this provides additional evidence for a deep-reaching fault, which affects the splitting. While the order of the layers is a counterargument, as the N-S layer (parallel to the strike of the fault) is found to be the lower layer in contrast to the shallow surface feature, the conditions for the two-layer analysis are not ideal. The number of measurements and the azimuthal coverage are limited. It seems logical that the laterally confined feature of the fault leads to an effect on events of only a limited azimuthal range. That, in combination with the limited azimuthal event coverage, can lead to an apparent two-layer model. Still, the strong effect on the splitting measurements requires a deepreaching fault likely extending into the mantle lithosphere. This is not supported by Moho contours of the suture, where a displacement, as seen at the surface, is not visible (Spada et al., 2013;Handy et al., 2015). Whether the measured two-layer effect in the western subdomain originates in the mantle or is influenced by the crust could be resolved with a longer measurement period of a similarly dense network. This would allow the application of techniques based on receiver functions or ambient noise, which may be able to resolve shallower features in more detail.

CONCLUSION
This study addresses the open questions regarding a possible slab break-off accompanied by a subduction polarity flip in the transition of the Central to the Eastern Alps. We approach this question based on tectonic implications given by shearwave splitting measurements at the dense SWATH-D complementary network and stations of the AlpArray backbone in the same area. As previous studies imply a twolayer anisotropy in the Eastern Alps (Qorbani et al., 2015a), we investigate this further and search for a possible two-layer anisotropic model. A depth-dependent anisotropic structure is also supported by the observed azimuthal variation of the single splitting measurements. The temporary character of the SWATH-D experiment and the associated limited amount of data do not allow for a conventional single-station two-layer analysis in this area. Therefore, we combine measurements of multiple stations based on a cluster analysis to resolve the splitting parameters of two layers. While allowing for further variation in up to six groups, the cluster analysis yields two main groups of similar splitting effects separated approximately at a longitude of 13°E. The western part is characterized by a lower layer of N-S direction and an upper layer with a fast direction of roughly 60°. The eastern part shows a lower layer with a fast direction of 60°and an upper layer with a fast direction of 110°, which is in agreement with previous observations (Qorbani et al., 2015a).
As the upper layer of the western part shows the same direction as the lower layer in the east, a common origin for the anisotropy seems likely. We interpret this as lithospheric anisotropy frozen in the European subducting slab with vertical coupling into the asthenosphere. Weak N-S aligned anisotropy in the west indicates either a return flow at the bottom of the retreating European slab pushed by the Adriatic plate or, alternatively, a deep-reaching Giudicarie Fault, which may produce an apparent two-layer effect. The abrupt change of the two-layer anisotropic pattern with an equally pronounced splitting of the upper layer in the east is evidence for a break-off event initiated 35 Ma ago (Handy et al., 2015). While the slab fragment is sinking into the mantle (a connection to the European slab is likely to remain), the gap provides an opening for asthenospheric flow from the Alpine Front to the Pannonian Basin.
While this interpretation provides evidence for the subduction of the European plate below the Adriatic intender with a break-off event in the eastern part of the Alps, the splitting measurements cannot constrain a possible subsequent subduction of the Adriatic plate below Europe, filling the gap (Lippitsch, 2003;Handy et al., 2015;Kästle et al., 2020). Therefore, further investigation is required to resolve the subduction polarity in the Eastern Alps.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data can be found here:

AUTHOR CONTRIBUTIONS
FL developed and performed the data analysis and wrote the first draft of the manuscript. FL and GR jointly discussed the results and their interpretation and revised the manuscript.