Shear wave splitting in the Alpine region

To constrain seismic anisotropy under and around the Alps in Europe, we study SKS shear wave splitting from the region densely covered by the AlpArray seismic network. We apply a technique based on measuring the splitting intensity, constraining well both the fast orientation and the splitting delay. Four years of teleseismic earthquake data were processed, from 723 temporary and permanent broad-band stations of the AlpArray deployment including ocean-bottom seismometers, providing a spatial coverage that is unprecedented. The technique is applied automatically (without human intervention), and it thus provides a reproducible image of anisotropic structure in and around the Alpine region. As in earlier studies, we observe a coherent rotation of fast axes in the western part of the Alpine chain, and a region of homogeneous fast orientation in the Central Alps. The spatial variation of splitting delay times is particularly interesting though. On one hand, there is a clear positive correlation with Alpine topography, suggesting that part of the seismic anisotropy (deformation) is caused by the Alpine orogeny. On the other hand, anisotropic strength around the mountain chain shows a distinct contrast between the Western and Eastern Alps. This difference is best explained by the more active mantle ﬂow around the Western Alps. The new observational constraints, especially the splitting delay, provide new information on Alpine geodynamics.


I N T RO D U C T I O N
The Alps are formed by the collision of the Adriatic and European plates, (e.g. Trümpy 1960;Stampfli et al. 2001;Schmid et al. 2004;Handy et al. 2010). While the surface expression of the mountain chain is well known, there remain many open questions that relate to its deeper portions, to the forces that drive the deformation, and to the deformation itself. Mantle flow must play an important role in mountain building, and in the formation of the wider Mediterranean tectonic system in general, for example Faccenna et al. (2014). Its role needs to be better understood. Luckily, the deformation that occurs during flow in the subsurface can result in preferred orientation of minerals, and becomes visible to seismic waves through the phenomenon of seismic anisotropy. Seismic anisotropy can thus be a tool for constraining the geometry of deformation (Silver 1996).
Several approaches exist for measuring seismic anisotropy. The most common one uses the phenomenon of 'shear wave splitting' that occurs for seismic body-wave phases, for example SKS phases. A seismic shear wave that travels through an anisotropic medium splits into a fast and slow wave. The splitting can be described by the two 'splitting parameters' φ and δt which stand for the azimuth of the faster S-wave velocity and the delay time between the fast and slow S-waves, respectively (Silver & Chan 1991).
We have processed waveform data recorded in the framework of the AlpArray project (Hetényi et al. 2018) including its marine component. That data set is more complete than any earlier data set for the region in terms of both data coverage and length of records at the seismic stations (Kummerow & Kind 2006;Barruol et al. 2011;Bokelmann et al. 2013;Salimbeni et al. 2013;Qorbani et al. 2015Qorbani et al. , 2016Salimbeni et al. 2018 andPetrescu et al. 2020a, among others). We have processed the data with an algorithm that is based on the splitting intensity method by Chevrot (2000). That technique resolves well both splitting parameters, and we apply it in an automatic mode. Processing the data automatically, without human interference, renders the results reproducible. This also allowed us to process a very large data set. We propose uniformly estimated splitting parameters for the whole Alpine region including the Ligurian Sea. Thus, our results are well-suited to discuss earlier propositions of seismic anisotropy beneath and around the whole mountain chain.
This paper presents the essence of the technique (Chevrot 2000) and its automatic application. Additional technical details are given

D E F O R M AT I O N A N D S E I S M I C A N I S O T RO P Y U N D E R T H E A L P S
Formation of mountain chains is deeply linked to plate tectonics. Collision of the Adriatic and European plates has caused the mountain belt to rise (Dewey et al. 1989;Stampfli et al. 2001;Schmid et al. 2004;Handy et al. 2010). The plate convergence is accompanied (or even driven) by asthenospheric flow in the upper mantle (Savage 1999;Lebedev et al. 2006;Long & Becker 2010). Understanding that flow associated with plate convergence is important therefore. The in situ deformation of rocks can give clues about the mantle flow. Seismic anisotropy can witness large deformation at upper-mantle depths, for example Bormann et al. (1993). Anisotropy preserves especially the most recent episode of deformation that becomes imprinted in the medium in the form of a preferential alignment of minerals parallel to the belt axis (Margheriti et al. 1996). The collision of the African and Eurasian plates and the resulting subduction of slabs are the key features that may have led to the formation of seismic anisotropy in and around the Alpine region.
In case of the Western Alps, Barruol et al. (2011) find a simple mountain-chain-parallel pattern of fast orientations, due to the mantle flow beneath the halted continental collision. Their observation complemented the earlier splitting measurements along the TRANSALP profile by Kummerow & Kind (2006) in the western part of the Eastern Alps. Bokelmann et al. (2013) have extended this pattern to the rest of the Eastern Alps. They have found an intriguing progressive rotation of fast orientations along the whole Alpine chain. A more extended picture of mantle deformation under the Alps was delivered by Qorbani et al. (2015), who measured the splitting parameters for more stations in the Eastern Alps. They have suggested the presence of two layers of anisotropy in the Eastern Alps, one of which could be associated with a detached slab of the subducted plate. Salimbeni et al. (2013) have analysed data from the Southern Alps and Appenines, and have tried to correlate fast seismic orientations with possible fossil and asthenospheric deformation that may be expected from the geodynamic situation of the region. Löberich & Bokelmann (2020b) have analysed data from stations in the Central Alpine region using a new technique developed by Löberich & Bokelmann (2020a), to infer the orientation of the flow plane directly from shear wave splitting observations. They suggested that the seismic anisotropy under Switzerland is most likely caused by asthenospheric flow around the Alps ('horizontal shear plane'). Petrescu et al. (2020a) have analysed AlpArray data from 113 stations across the Central Alps, and have suggested that fast orientations are mostly mountain-belt parallel. They suggested that spatial variations of fast orientation agree with some of the geological divisions around the Alps. Link & Rümpker (2020) have analysed differential Ps -XKS splitting from the AlpArray data set, to simultaneously infer parameters for the crust and the upper mantle. They find a spatial pattern with high degree of spatial coherence in fast orientations. The anisotropy of the transition region from the Alps to the Carpathians and the Pannonian Basin was studied by Qorbani et al. (2016). They showed, how the fast orientation pattern continues to the east of the Alps. The Carpathian-Pannonian region was then also investigated by Petrescu et al. (2020b), who again extended our knowledge of the anisotropy pattern further east. Some constraints on seismic anisotropy have also been provided using P-waves, for example Hua et al. (2017), mostly agreeing with the earlier suggestions.
What stands out in these studies is that although the surface has undergone complex geological deformation, the pattern of fast orientations in the Alpine region appears to be relatively simple. The splitting delay times has generally received less attention than the directional component, since the delay is usually subject to larger errors. However, the strength and depth of the anisotropic perturbation can be related to the size of the splitting delay between the two phases (fast and slow), and the splitting delay would clearly add important new constraints. Our current study follows the effort of the above-mentioned papers and presents a broad view of anisotropy beneath the Alps and neighbouring regions.

Data
AlpArray was a large multinational project in Europe, in which 55 institutions from 17 countries gathered to jointly study the Alpine region (Hetényi et al. 2018). The main infrastructure part of the project was the broad-band seismic network which was set up by the European institutions to instrument the wider Alpine region with almost 300 mobile seismic stations to fill the gaps in the permanent networks coverage. The AlpArray network covered the area with an evenly spaced station grid, with station separations of not more than 52 km. The main portion of the seismic experiment was carried out between 2016 and 2018, but some stations have been in the field before and after. Details about installation and performance of part of the AlpArray network can be found in Fuchs et al. (2015Fuchs et al. ( , 2016. The AlpArray project had also a marine component: 30 ocean bottom seismometers (OBS) were deployed in the Ligurian Sea between June 2017 and March 2018. The AlpArray region was defined as 250 km distance from a smoothed 800 m contour line around the Alps, see fig. 2 in Hetényi et al. (2018). For the purpose of our current study, we generalized the region by circumscribing an ellipse-like shape around the AlpArray region. Altogether, our list comprised 765 broad-band seismic stations in the area of interest taking advantage of the neighbouring permanent networks as well as of the ocean bottom seismometers.
We have gathered earthquake information from the monthly Lamont-Doherty Earth Observatory catalogue. 1 We selected events of magnitude greater than M w = 6 that have occurred between January 2015 and March 2020 in epicentral distance between 85 • and 135 • from the centre of the AlpArray. That list contains 304 events. In total, we collected 151 000 records. For each event, the data was available roughly at 500 stations on average. In other words, we processed around 200 records from each station on average. After the initial pre-processing and data quality checks, the list of exploited events shrank down to 250 earthquakes and the list of used stations to 723. Depending on the selection criteria used, successful measurements were obtained on average for 18-26 events per station for the most strict and most relaxed criteria (see below), respectively. It means, that only 9-13 per cent of all events recorded at every particular station (around 200 events) yielded useful splitting measurements. However, in general, different events were recorded at different stations, and hence the proportion of usable records in the entire data set is higher than 13 per cent.

Shear wave splitting
Several approaches exist for determining shear wave splitting parameters. Bowman & Ando (1987) proposed a technique based on cross-correlation of horizontal components. The difference between the lag times of maximum and minimum cross-correlation represented the splitting time. The rotation angle with the highest crosscorrelation represented the anisotropy orientation. Silver & Chan (1988) developed the well-established transverse energy minimization technique. In a variant of this technique, so-called 'eigenvalue technique', the second eigenvalue of the covariance matrix can be minimized (Silver & Chan 1991). That technique forms the core of SplitLab, a powerful toolkit written in Matlab R by Wüstefeld et al. (2008). These techniques operate on a single event-station pair and produce an estimate of the two splitting parameters individually from each record.
A rather different approach has been suggested by Chevrot (2000). He determined only a 'splitting intensity' for each eventstation pair (rather than splitting delay and fast orientation). The splitting intensity method utilizes a set of SKS measurements made at a station, each one providing an individual splitting intensity. The backazimuthal dependence of splitting intensities from different events are then used to infer the splitting parameters δt and φ. A sinusoid is fitted to the individual splitting intensities as a function of backazimuth. Maximum of that sinusoid represents the splitting delay time, and zero-crossing of that sinusoid represents the fast axis azimuth.
The splitting intensity measurement makes use of the approximate relation between the two horizontal components of a vertically incident S-wave travelling through a homogeneous anisotropic medium. The relation has been presented by Silver & Chan (1991) and Monteiller & Chevrot (2010): the linearly polarized wavelet w(t) splits into two S waves with a delay time δt between them (fast and slow). If that delay is much smaller than the period, the radial and transverse components, R(t) and T(t), can be approximately written as The transverse component can be seen as the scaled derivative of the radial component R (t): The splitting intensity is s = δt sin 2β, and the splitting delay time δt is the maximum it reaches. The phase of that function then relates to the fast axis orientation φ as β = φ BAZ − φ. β represents the angle between the earthquake backazimuth φ BAZ and the fast axis orientation φ. More details on the derivations and on computational details can be found in the Appendix and in Chevrot (2000) and Monteiller & Chevrot (2010). Methods operating on single station-event pairs extract much information from the waveform data, and the backazimuthal variation of splitting delays can be analysed; yet they tend to be sensitive to noise in the measurements (Vecsey et al. 2008) and the error bars of the splitting delay are often too large to be used in a meaningful way. The Chevrot-technique is promising, in that it extracts a simple linear measure from the data of all earthquakes recorded at the station. It has the potential to constrain the splitting delay better, which is helpful for studying the lateral and vertical variations of seismic anisotropy. Furthermore, the splitting intensity measurements can subsequently be used for anisotropic tomography studies (Chevrot 2006). It can yield information on anisotropy beyond single-layer anisotropy, even though that has been assumed in the derivation of the splitting intensity. Splitting intensity method can be applied at stations which provide sufficient number of suitable earthquake records, meaning at stations with long enough deployment. AlpArray project supplied hundreds of earthquakes of adequate magnitudes in the required distance range. Together with the high density of stations, the splitting intensity method is especially suitable to be applied in the Alpine region.

Noise bias and Wiener filter
The use of multiple events and the simple set up of the splitting intensity method provide some computational stability, which is an important feature for automatic processing. However, we think there is a bias in the presence of noise, and thus a need to deal with it. The issue emerges due to the inherent assumption of a noise-free radial component (see the Appendix). Synthetic tests under various noise conditions confirmed this bias. In the Appendix, we suggest two ways of alleviating this problem. The Wiener filter can suppress this issue sufficiently, by eliminating the assumed additive random noise (Brown 1983). This study thus makes use of the Wiener filtering, as has been also done in Monteiller & Chevrot (2010). The Wiener filter is applied to the bandpass-filtered data (10-100 s), see below.

Automatic data selection and processing
Our computational procedure uses the Python TM library called Split-Wave which was implemented by Walpole et al. (2014). It utilizes the above-mentioned single-event-station-pair methods, which we enhanced with our implementation of the splitting intensity method of Chevrot (2000), to process the data in a completely automatic mode.
We first corrected the records for the instrument response. The data were linearly detrended and downsampled to 10 Hz. Afterwards, the waveforms were bandpassed between 10 and 100 s with a zero-mean Butterworth filter. From each record, we selected a time window of 5 min centred around the SKS theoretical arrival times. These arrival times were calculated by the tau-p method (Buland & Chapman 1983;Crotwell et al. 1999) using the ak135 velocity model (Kennett et al. 1995). Splitting measurements were performed on shorter time windows that were selected by the z-trigger detector according to Swindell & Snell (1977) (an alternative to STA/LTA algorithm).
On the input of the trigger, we had squared amplitudes of the horizontal motion in that pre-selected 5-min window (sum of the squares of the horizontal components). The trigger characteristic function was calculated as a sum of the input amplitudes over a moving window of 100 samples length (10 s). From that function, its mean was removed and it was divided by its standard deviation to normalize it. The thresholds for triggering were set to 2.2 for trigger on and to 0.8 for trigger off. From the resulting triggered windows, we selected the one which was the closest to the theoretical SKS arrival and we only used it, if it was within 30 s from that theoretical arrival. After that, the Wiener filter was applied to the rotated horizontal components in that selected shorter time window. Fig. 1 shows two examples of 5-min windows of horizontal components (black and blue in panels (a) and (d)) with predicted SKS arrival (red bar) and with the triggered time window (magenta bars). Panels (b) and (e) show details of the triggered windows with Wiener-filtered horizontals and with the derivative of the R component. To test the stability of the time window selection and of the retrieved parameters, the measurement was performed 10 times on each waveform. We were increasing the length of the used time window by adding increments of 0.5 s length to each end of the original time window. This was repeated up to adding 4.5 s to each end of the original window. The longest window was hence 9 s longer than the original triggered window. All 10 values of splitting intensity measured in all these time windows were subsequently used in estimating the splitting parameters, if they passed the selection criteria (see below). Fig. 1, panels (b) and (e), show example of time windows extended by 4.5 s to each side (in both examples) of the originally triggered window-see the filtered traces exceeding the magenta vertical bars. Applied to all 151 000 records, this process led to a total of 90 000-154 000 automatically processed time windows with successful splitting intensity measurement depending on the selection criteria used. On average, we used 8.9-9.2 time windows per event per station.

Calculating the splitting parameters
After the splitting intensities are gathered, plotting these values against event azimuth shows considerable scatter, see red and black symbols in Fig. 2, for two stations in panels (a) and (c) as an example. This scatter is strongly reduced by using selection criteria. These include an upper limit of |s| < 2 s for the splitting intensity, and a cosine similarity criterion with the scalar product R · T. The cosine similarity C was set at a threshold of greater than 0.7, to select only SKS phases which show a reasonable correlation between the transverse component and the derivative of the radial component (Foreman et al. 2014). The cosine similarity was estimated for every time window independently. Hence, in general, some of the 10 time windows could fall below the cosine similarity threshold of 0.7 or its splitting intensity can exceed the upper limit of 2 s and could not be used then. The final splitting parameters are then estimated by fitting a sinusoidal curve through the splitting intensity values that satisfy the selection criteria (see the red symbols in Fig. 2). Another requirement for accepting a station is a minimum coverage of at least four of the non-overlapping 15 • azimuthal windows. The dominance of the 180 • periodicity in Fig. 2 suggests that the scatter in the black symbols is usually due to random noise, rather than due to unmodelled structural effects. Errors of splitting parameters (fast axis orientation and splitting delay) are estimated from fitting the sinusoidal splitting function to the measurements of splitting intensities for all events and all lengths of time windows at a given station. The two examples in Fig. 2(a) and (c) are given with their respective errors for both parameters. Errors for all the stations which yielded splitting parameters are given and discussed later.

Azimuthal harmonic order
The analysis is based on the simple assumption of a single nondipping horizontal layer for which the splitting parameters are inverted. Our results are hence yielding an effective anisotropy estimated over the possible vertically variable anisotropic layers. To test the strength of this assumption one can consult the Fouriertransformed splitting intensities with respect to azimuthal order k, see Fig. 2(b) and (d). For the ideal non-dipping single layer case, the second-order term will dominate (180 • periodicity). Synthetic tests by Chevrot (2000) have shown that a dipping layer of more than 30 • would cause a shift to higher amplitudes of the first-order term (360 • periodicity). The higher-order terms (k > 2) suggest either more complicated anisotropic structure or the effect of noise. However, larger values at higher orders can easily arise as an artefact of the uneven backazimuthal coverage. Higher-order terms need to be handled carefully therefore, and they are not very useful as practical diagnostic. For further information one may refer to eqs (8) and (9) in Chevrot (2000).

General view
The distribution of automatically determined splitting parameters in Fig. 3 depicts fast orientations and splitting delay times for the Alpine region and also around it. The map shows results from 512 stations, which passed our above-mentioned criteria for data quality and the three strict splitting intensity measurement requirements. The generally good spatial coherence of the fast axes orientations is striking, especially considering that the data processing has been entirely independent from station to station. The most significant observation is the large-scale pattern of mountain-belt-parallel fast orientations, which also continues in the outer Alpine regions. Such a pattern had already been described by previous studies, see Barruol et al. (2011), Bokelmann et al. (2013), Qorbani et al. (2015), and partially also by Salimbeni et al. (2018) and Petrescu et al. (2020a). The pattern shown in Fig. 3 is not just confined to the Alps; it extends off the mountain range especially to the north. The distribution of splitting delay times, given by colours and lengths of the bars in Fig. 3, also shows good spatial correlation. There appears to be some correlation with topography, with the strongest splitting occurring under or adjacent to the mountain chain, and the weakest splitting occurring in areas remote from the mountain chain and in the Po Plain. There are other regions outside the Alps, which exhibit interesting features: high splitting delaytime anomalies for the Swiss Prealps, which extends northward to the Rhine Graben; the northern Apennines display high delay times for the back arc and in the foreland; the very consistent pattern in fast orientation and (weak) anisotropy in large parts of the Bohemian Massif; and a strong and coherent splitting at the southeastern edge of the Alps. The spatial coherence of both the fast orientations and of the splitting delay times for neighbouring stations gives reason to believe that the measurement reflects real larger-scale (compared to station distance) structural properties. Our study offers new and more complete spatial coverage of both splitting parameters than earlier studies. Fig. 4 shows four measures explaining the quality of the splitting parameters presented in Fig. 3 for each station. Fig. 4(a) presents the R 2 -value. It is defined as 1 − (S res /S tot ), where S res is sum of squares of the residuals (the difference between the data and the prediction of the sinusoidal model) and S tot is total sum of squares    . Panel (a) shows splitting intensities (SI) observed at the permanent station BNALP as a function of event azimuth. The size of the dots is scaled by the signal-to-noise ratio. Larger dots indicate more reliable measurements. The red dots show measurements which satisfy the similarity criterion and at the same time also the upper limit criterion (and have then been used for fitting a sinusoidal curve-the splitting intensity function). The resulting splitting parameters φ and splitting delay δt are indicated at the bottom of the plot together with their uncertainties. Dots at the same azimuth are measurements from the same event sampled at multiple time windows. Panel (b) shows the amplitude of the harmonic spectrum as a function of azimuthal order in black for all measurements and in red for the selected satisfactory measurements only. The second term dominates in both cases, but it is much more prominent after the selection criteria are applied. Panels (c) and (d) show the same for the temporary station A141A. of the difference between data and their mean. That ratio of the two sums of squares is equal to the ratio of unexplained variance to data variance. Possible values of R 2 range between 0 and 1, with 1 representing the case of perfect data fit. We see, that there are clear spatial clusters of higher and lower R 2 -values. Panel (b) shows the number of time windows with successful measurement yielded at each station. We again see that these numbers are higher in the centre of the area of interest and lower at the edges. Panel (c) of Fig. 4 presents errors of splitting delays and panel (d) presents errors of fast axis orientations. Lower errors are associated to a higher number of events. The errors shown in panels (c) and (d) represent the uncertainty of a fit of the sinusoidal curve into the splitting intensity measurements for a given selection criteria: upper limit of |s| < 2 s, cosine similarity C ≥ 0.7 and four 15 • azimuthal windows covered. They do not reflect other possible scenarios of selection criteria.

Variation along the Alps
First, we describe the variation of splitting parameters along the Alpine chain, and especially the variation of the orientations of fast axes. A progressive rotation of fast axes along the Western Alps had been identified by Barruol et al. (2011), and such a progressive rotation also appears in the new measurements (Figs 3 and 5a). This is noteworthy, since we are comparing the manually retrieved set of splitting parameters from the earlier study with the automatically retrieved set of parameters of the current study. The similarity confirms the earlier suggestion of progressive rotation of fast axes in the Western Alps, and it also suggests that the automatic method provides stable results.
We use the Alpine axis defined in Bokelmann et al. (2013), and split our 512 stations shown in Fig. 3 into those up to the distance of 125 km from the Alpine axis (red dots in all three panels of Fig. 5; 227 stations) and those further from the axis (black dots in all three panels of Fig. 5; 285 stations). The closer stations roughly cover the region investigated by Bokelmann et al. (2013) and by other studies referenced therein. Panel (b) shows clear similarity of the resulted orientations between our selection of closer stations and the results from Bokelmann et al. (2013), see the green solid line taken from fig. 3 of the latter paper. In both data sets, the progressive rotation comes to a halt at around 600 km along the profile of the Alpine axis (see Fig. 5(a) and (b)).
The variation of our fast orientations is slightly less similar to Bokelmann et al. (2013) in the easternmost part, starting from around 1100 km, see Fig. 5(b), where the new results (red dots) start to have lower azimuths of orientations than the general trend from previous result (green line). For the easternmost Alps, Qorbani et al. (2015) suggested the presence of two anisotropic layers. We will discuss that more complex anisotropy later. It is interesting that the rotation of fast orientations can not only be found in the region close to the mountain belt (red dots in Fig. 5), but also to the north of the Alps (dots shown by black). Comparison with the orientation of the Alpine axis itself confirms that the fast orientations are not strictly parallel to the Alps; there is a constant offset of around 25 • with respect to the orientation of the Alpine axis, until the distance of about 1000 km (Fig. 5(b)).   fig. 3 in Bokelmann et al. (2013) suggests that the error bars of the fast axis orientations are smaller in our current study than in the previous measurements. Splitting delay times have been considered to be relatively badly constrained in the earlier studies, see the comparison with Barruol et al. (2011) below. In our study, the uncertainties of the splitting delay times obtained by the new splitting intensity technique are also much smaller than in the previous works (see also the tables with the splitting delay errors in Bokelmann et al. (2013), Qorbani et al. (2015) and Barruol et al. (2011)). Indeed, the average error among the 512 stations presented in Figs 3 and 5 is ±2 • for the fast axis orientations and ±0.1 s for the splitting delay times.
The lower uncertainties of our current results allow us to evaluate also the variation of the splitting delay times along the Alpine chain . Fig. 5(c) again shows the set of 512 stations split into 227 close stations (red) and 285 more distant stations (black). As the splitting delays have larger variation and do not exhibit a clear trend, we added running averages of the splitting delays for the closer stations (dark red solid line) and for the more distant stations (grey solid line). In general, the closer stations have slightly higher splitting delays than the more distant stations. In the closer-station set, we see local splitting delay maxima between the distances of 600-700 km and between 1000-1100 km. These correspond to high-delay regions on the eastern side of the Swiss Prealps, and in South Tyrol (Dolomites) seen in Fig. 3. After the distance of 1100 km (easternmost Alps), the splitting delays tend to slightly diminish.

Comparison with Barruol et al. (2011)
Now, we will focus on the variations of measured anisotropy in the western part of the Alps. We compare our splitting parameters with those from Barruol et al. (2011), in the following abbreviated as B11. The red bars in the map in Fig. 6(a) indicate fast axis orientations from B11. Black bars are measurements from our current study. The data sets overlap at 35 stations, if we consider as 'the same' station also three pairs of stations, which were just renamed or shifted by only hundreds of meters since the B11 study. We can use those 35 stations, labelled by their current names in Fig. 6, to analyse the differences in fast orientations and splitting delays between the two studies. From Fig. 6, panels (a) and (b), we see that the differences peak near zero; the results are thus in general consistent between the two techniques. The differences form normal distributions; just station BRANT near the Swiss-French border appears to be an outlier. Errors in our current study are smaller than those in B11, see the bar charts in panel (c) for the distribution of splitting delay errors and in panel (d) for the distribution of fast orientation errors. Measurements in B11 (red) have larger mean error. For the splitting delay times, errors in B11 peak at 0.2 s while our current study has the peak (black histogram) at 0.03 s. Similarly, fast orientation errors in B11 peak around 10 • while our splitting intensity method measurement has most of the errors concentrated between 1 • and 2 • for the given set of 35 overlapping stations. The distributions of differences around zero in Fig. 6(b) correspond to the error distribution of B11. One reason why the errors differ so much between the splitting intensity method used in our study and the technique used in B11 is that the former is determined from multi-event measurements at every station, while the latter refer to single event measurements. They can therefore not be directly related.

D I S C U S S I O N
We have obtained a spatially coherent pattern of fast orientations in and around the Alps, which complements earlier splitting measurements. Our study provides a complete distribution of splitting parameters in a reproducible fashion. Both the fast axis orientations and the pattern of the splitting delays show good spatial coherence.
We enhance in particular the retrieved information on the splitting delay, which is systematically and uniformly constrained. That information is crucial when it comes to deriving the strength of anisotropy in the upper mantle.

Splitting intensity, other techniques and selection criteria
In this study we have primarily used the splitting intensity technique that had been proposed by Chevrot (2000) and Monteiller & Chevrot (2010), and we have used it in a fully automatic implementation. More commonly used methods, such as the transverse minimization technique by Silver & Chan (1988) and the cross-correlation technique by Bowman & Ando (1987), are more influenced by noise in the data, as was already clear in earlier attempts establishing automatic procedures (Wüstefeld et al. (2008), Teanby et al. (2004)). The splitting intensity technique instead has turned out to be relatively insensitive with respect to the choice of the time window, and therefore more suitable for automatic procedures.
Comparison with the earlier studies by Barruol et al. (2011) (Fig. 6) in the Western Alps and Bokelmann et al. (2013) (Fig. 5) along the Alpine arc have generally shown good agreement, and our results are consistent also with previous results on the Eastern Alps when one anisotropic layer only has been considered ( fig. 2 in Qorbani et al. 2015).
Earlier studies have concluded that multiple layers might be present when the uncertainty in the splitting parameters becomes large (e.g. Lamarque & Piana Agostinetti 2020; Petrescu et al. 2020a), and simplistic models might be too rough to explain the predominant anisotropy then. The splitting intensity technique allows inspecting the azimuthal variation spectrum (as in Fig. 2), and in theory, we can define the areas where a single anisotropic layer satisfies the observed data on that basis, versus the areas where a more complex layered anisotropic structure is required by the data. We have found however that the higher orders of that spectrum are easily perturbed, if the azimuthal distribution of splitting intensities is not dense. We have therefore chosen to look at the goodness of fit represented by the R 2 -value in Fig. 4 instead; in other words, we inspect the fit of the predicted sinusoidal curve to the data. For the eastern part of the Alps, where a two-layer anisotropy has been suggested by Qorbani et al. (2015), we indeed find lower R 2 . This may in principle be due to unmodelled complex anisotropy, or due to higher noise in the data (worse measurement conditions). It is notable in Fig. 4(a) that the central part of the study region is associated with better data fit than around it, not only to the east, but also to the south and the west. The poorer data fit in the south and west could be due to more complex anisotropy as in the east, however, there are no published suggestions of complex anisotropy for these regions. The poorer data fit could possibly have also another reason. Fig. 4(c) and (d) indicate that the errors in splitting delay and fast orientation are lower in that central region as well. Fig. 4(b) shows the number of measurements that were used for the analysis, and it becomes apparent that the central region is also associated with more measurements that could have been recovered. This suggests that the primary effect is that measurement conditions are better in eastern Switzerland, western Austria and surroundings. This can be due to the type of station installation (permanent or temporary), and due to the local settings like the geological subsurface (on sedimentary basin or not), etc. We believe that all four measures that are shown are affected similarly by measurement conditions (and subsurface complexity)-hence the correlation between the four measures. In the regions of lower R 2 -value and higher uncertainties, more complex anisotropy is possible, especially in the east, but that is not required from these observations. Fortunately, the AlpArray experiment has instrumented the area for several years, and a fair number of strong earthquakes from various azimuths has been recorded; yet insufficient backazimuthal coverage is a limiting factor, which is reflected in the formal uncertainty (Fig. 4). Enforcing a minimum backazimuthal coverage is one of the three selection criteria, together with an upper limit of 2 s for the splitting intensity, and a minimum value for the cosine similarity criterion (see above). In the following, we discuss the dependence of the splitting map on whether we enforce these criteria or not, to test the stability of the outcome. We show the splitting maps for weaker sets of selection criteria in Fig. 7, with imposing the coverage of only at least three of the 15 • azimuthal windows (compare it with four azimuthal windows required in Fig. 3). First, we kept the cosine similarity the same as in Fig. 3 to be greater than 0.7, see Fig. 7(a). In this case, we retrieved the splitting parameters from 576 stations-compared with 512 stations in Fig. 3. Releasing the azimuthal coverage criteria and keeping the cosine similarity threshold means that the 512 measurements from Fig. 3 are obtained the same in Fig. 7(a), and 64 additional measurements (stations) are retrieved. We see, that the additional measurements have variable splitting delays; the additional measurements with high splitting delay do not match the overall pattern of fast orientations (outlying long bars with yellowish colours), while the additional measurements with low splitting delay match the orientation pattern well (short bars with dark colours). Fig. 7(b) again requires the azimuthal coverage of three azimuthal windows for the onland stations and there is no requirement for the OBS stations. In addition, we lowered the cosine similarity threshold to 0.5 for both onland and OBS stations. We see many more additional stations providing splitting parameters (650 in total). Generally, the whole image is shifted to lower splitting delays. This effect is described in the Appendix: a lower cosine similarity threshold leads to accepting more data, which in turn corresponds to a higher noise level in the measurements-and a stronger bias to smaller splitting delays. The decrease of splitting delays can be larger than the measurement error for each of the individual cosine similarity scenario, suggesting that the formal errors are optimistic. Note that in the case of no azimuthal requirement for the OBS stations in the Ligurian Sea, the records from 12 of them delivered splitting parameters, even though the OBS stations were deployed only for a limited time period of 9 months. Compare with Fig. 3, where there was no OBS measurement retrieved, and with Fig. 7(a), where only one permanent OBS station yielded the splitting parameters. The fast orientations of the 12 OBS measurements show larger scatter, but on average they connect well with the more stable estimates onland in Italy, Southern France and Corsica. This is an interesting result, suggesting a common origin of the seismic anisotropy under the Ligurian Sea, Southeastern France and Corsica, indicating a likely asthenospheric origin, as suggested by Barruol et al. (2004) and by Lucente et al. (2006).
Note that we have not removed any outliers (with respect to coherence with neighbouring stations) from the two splitting maps in Figs 3 and 7(a), to accentuate that the splitting maps are determined using a well-defined set of criteria (reproducible). We have removed six outliers with δt > 2.6 s from Fig. 7(b) (three onland and three OBS). All other δt are smaller than 2.1 s. The most important observation is, that even at the last map of splitting parameters (Fig. 7(b)) retrieved with very loose selection criteria, the spatial coherence of both the fast orientations and splitting delays remains preserved. In fact, the weaker criteria in Fig. 7(b) create only few outliers in the fast orientations (only three removed onland), less than in Fig. 7(a), where there are around 12 obviously outlying stations onland.
Apparently, the larger number of measurements (154 400 in panel (b) compared to 93 400 in panel (a)) has a stabilizing effect on the fast orientations, while its stronger noise creates a bias in the splitting delays to lower values. The splitting delays in Figs 3 and 7(a) are probably more representative of the true strength of the anisotropy effect though. The good continuity of fast orientations within the Alpine region also confirms that sensor misorientations are provided properly in the station metadata.
To see what is the effect of the selection criteria on the results, we calculated splitting parameters for many different scenarios perturbing two of the selection criteria. We selected 21 cosine similarities (see eq. 4) in the range 0.4-0.8 with a step of 0.02 and for each of these cosine similarities, we used 4 different azimuthal windows requirements (at least 1, 2, 3, and 4 of the 15 • windows covered). Hence, we set up 84 different scenarios. For each station, we calculated the standard deviation of all time delays and fast orientations obtained over these scenarios. These standard deviations are, on average over all stations, 0.17 s for the time delay and 8.4 • for the fast orientations. As the averages are biased towards higher values by the outlying stations with very large spread of values, we also calculated the median values for each splitting parameter over all stations. For the delay times, the median uncertainty is 0.11 s and for the fast orientations, it is 3.4 • . These values show that different selection criteria do not have a significant effect on the results. Still, lower cosine similarity systematically shifts the splitting delays towards lower values but allows more data to be used. Higher cosine similarity threshold lowers the number of measurements, but allows for more precise splitting delay determination.

Interpretation and geodynamic inferences
The fast orientations that we observe in this study agree with those that have been already seen in earlier studies (Barruol et al. 2011;Bokelmann et al. 2013;Qorbani et al. 2015Qorbani et al. , 2016Petrescu et al. 2020a). Our results confirm and extend the picture that we already have about the Alps. The new map of fast orientations shows an excellent spatial coherence, which is remarkable, since stations have been processed entirely independently from each other. Each station provides a smoothed image of subsurface anisotropy, according to the lateral extent of the Fresnel zone. For a dominant period of 15 s, the Fresnel zone would have a diameter of ca. 150 km at a depth of 150 km (e.g. Alsina & Snieder 1995). Since this averaging width is smaller than the arc of the Western Alps, it is clear that the rotation of fast orientations in the Western Alps is a real feature that requires interpretation. In and close to the Alps, fast orientations tend to be roughly parallel to the mountain chain, except for an intriguing nearly constant angular offset in the Western and Central Alps (see Fig. 5) that will be addressed below. This study also adds information about the fast orientations under the Ligurian Sea, which align with the spatial patterns seen in the areas around (i.e. Italy, southern France and Corsica as seen in Salimbeni et al. (2008Salimbeni et al. ( , 2013Salimbeni et al. ( , 2018, Barruol et al. (2011), Plomerová et al. (2006 and Barruol et al. (2004)).
We have interpolated the splitting parameters onto a regular grid across the Alpine region for the measurements retrieved for the sets of selection criteria that are shown in Figs 3 and 7(b). The interpolated versions of these splitting measurements are shown in Fig. 8. The combination of azimuthal coverage and cosine similarity requirements used is marked on the respective maps. Both maps are important: Fig. 8(a) provides a better representation of splitting delay, while Fig. 8(b) gives a better representation of fast orientations, due to the smaller number of outliers in the data set that it is derived from. Beside the bias in splitting delay, the interpolated patterns are quite similar in essence.
What is particularly new in this study with respect to previous studies is the stable pattern of the strength of seismic anisotropy (splitting delay) for the region, which provides new aspects for interpreting the Alpine geodynamics. At large scale, we note the prevalence of strong splitting delays clustered in/near the mountain chain, while weak splitting is retrieved at larger distances from the mountain axis, and towards the west and east away from the Alps (i.e. the Massif Central, the Bohemian Massif, the Pannonian Basin and the Dinarides). This correlation between splitting delay and topography (Figs 3, 7 and 8) indicates that some of the deformation causing seismic anisotropy is due to the Alpine orogeny, rather than fossil anisotropy derived from earlier deformation times, like for example during the Variscan orogeny.
Another striking pattern at large scale is the asymmetry in splitting orientations between the western and eastern terminations of the Alps. While splitting delays are large and fast orientations are rotating rapidly around the entire western end of the mountain chain, the eastern end is characterized by weak splitting and slow spatial rotation of the fast orientations. The fast orientations and splitting delay times which we obtain for the stations in the Western Alps, in northwestern Italy, and now also for the Ligurian Sea (see Figs 6 and 7), agree and support the hypothesis of toroidal flow around the Western Alps that has been discussed by Barruol et al. (2011), which is driven by the rollback of the Adriatic plate towards E, and possibly aided by the rollback of the European plate slab in the Western Alps towards NW. Fig. 9 compares the splitting map of Fig. 3 with two tomographic depth slices from the model of Kästle et al. (2018). Panel (a) corresponds to 100 km depth and panel (b) corresponds to 150 km depth. Several of the low-velocity zones in the western half of our study region are associated with strong anisotropy, for example to the west of the Apennines, the Ligurian Sea, and to the west of the Western Alps. These are regions where one would expect the 'toroidal' flow around the Western Alps (Barruol et al. 2011) and toward the Apenninic slab to be particularly strong, and this flow would indeed be expected to be associated with high temperature (low seismic velocity) and strong seismic anisotropy. Evidence from Löberich & Bokelmann (2020b), based on resolving the flow plane orientation from seismic data, suggests anisotropy related to asthenospheric flow.
Within the Alps, there are two regions which stand out by their consistently strong splitting delays (up to about 1.7 s in Fig. 3). These are located at the eastern side of the Swiss Prealps, and in South Tyrol (Dolomites; Figs 3, 7 and 8, see also the two local maxima of the dark red line average in Fig. 5(c)). The splitting orientations in both these anomalies are consistent and are SW-NE oriented. Comparison with the tomographic depth slices in Fig. 9 (Kästle et al. 2018) shows less correlation though. The western splitting maximum (marked by green ellipsis and labelled as 'Swiss Prealps' in Fig. 9) is located within or at the northern edge of a positive velocity anomaly. This may suggest that the splitting anomaly is either associated with the subducting lithosphere in the upper mantle, or it is connected with the flow around the Alps. The eastern splitting maximum (marked by green ellipsis and labelled as 'South Tyrol Dolomites' in Fig. 9) is located to the SE of the anomaly. There is a nearly 100 km gap of lower splitting delay (decrease of around 40 per cent) between these two (best visible in Fig. 8). This band of weaker splitting is oriented SW-NE and has been noted before (e.g. Löberich & Bokelmann 2020b) using earlier data that are not included in this study.
Above, we have alluded to the systematic offset between fast orientations and the orientation of the mountain chain. This offset suggests that the anisotropy is not primarily related to the Alpine lithosphere. Asthenospheric flow seems more likely. Such an asthenospheric toroidal flow around the westernmost Alps has been discussed in earlier studies, that is Barruol et al. (2011), Bokelmann et al. (2013 and Salimbeni et al. (2018). It has been interpreted as an effect caused by the eastward retreat of the Apennines slab (Jolivet & Faccenna 2000;Barruol et al. 2004;Lucente et al. 2006). The counter-clockwise rotation of the Adriatic plate, around a pole located between the Po Plain and the Western Alps (Calais et al. 2003), has led to the opening of the Ligurian Sea and to the counter-clockwise rotation of the Corsica-Sardinia block. The resulting asthenospheric flow is recorded in the Apenninic backarc.
Tracking flow lines from the seismic anisotropy measurements (especially the interpolated map in Fig. 8b), and assuming that fast orientations represent relative motion orientation, shows a more interesting pattern than simply being parallel to the outline of the Alps. The asthenospheric flow is modulated by the presence of the Alpine lithosphere, and that bounding lithosphere is moving itself. The spatial relation with tomographic anomalies is intriguing, but not simple. Reconstructions of Adriatic plate motion (e.g. Handy et al. 2015) show motions in the Western Alps that are perpendicular to the fast orientations, and this is probably not a coincidence. It can be expected that those relative motions (between Adriatic and European plate) have a relation with mantle flow below. This is especially the case for rollback, and such a scenario has recently been proposed for the Western Alps (Schlunegger & Kissling 2015). That explains a variety of observations about the Western Alps that were previously difficult to reconcile. In (oceanic) subduction zones it is quite common to find trench-parallel mantle flow under the slab (Long & Silver 2009), and this may be also the case here. Indeed, it appears that this kind of flow would embed well into the general toroidal-flow pattern around the Western Alps. Probably the seismic anisotropy pattern reported here will be useful for constraining geodynamic flow models, and arriving at a self-consistent view of flow under that larger region.
In the Eastern Alps, we detect E-W splitting orientations rotating towards east to ESE-WNW directions. The splitting delays are lower here compared with the Central and Western Alps. Our fast orientations are in agreement with the average splitting parameters reported by Qorbani et al. (2015). Their average splitting delays were larger though, and they attributed the effect to the presence of two layers of anisotropic material, the deeper one of lithospheric origin in the detached slab, and the upper one of asthenospheric origin. The splitting parameters in this study do not allow validating this model in a simple way. They can be used however in a subsequent anisotropic tomographic inversion, to elaborate on the presence of two anisotropic layers. Such a model is of interest, since it can help to better understand the depth extent of the deformation associated with the escape of the tectonic units towards the east, driven by the opening of the Pannonian basin (Gutdeutsch & Aric 1988;Ratschbacher et al. 1991).
Several models of subduction under the European Alps have been proposed, considering either south-directed subduction for both Western and Eastern Alps (Lüschen et al. 2004), a subduction polarity flip (Lippitsch et al. 2003;Kissling et al. 2006;Zhao et al. 2016;Hua et al. 2017), a 'shorter' subduction beneath the Eastern Alps, reaching maximum 200-250 km depth (Lippitsch et al. 2003;Kästle et al. 2018) or a subduction extending down to the mantle transition zone (e.g. Koulakov et al. 2009;Dando et al. 2011;Mitterbauer et al. 2011). Even though there is no clear consensus on the origin, shape and polarity of the subduction in the Eastern Alps yet, all cited previous works observe that subduction and they agree on the location of the divide and on the lack of a shallow slab (observing the fast velocity anomaly for depths larger than 250 km), at longitudes east of 14 • (Lippitsch et al. 2003;Mitterbauer et al. 2011;Zhao et al. 2016;Kästle et al. 2018). Our observations are coherent with the divide observed by the tomographic models. When looking at the area as a whole, and at the traces of the tectonic lines ( Fig. 8 in particular), we see resemblance between these orientations notably in the divide between the two strongly anisotropic regions, in the vicinity of the Giudicarie line. This raises questions on the nature of surface-mantle coupling, and on how deep processes drive deformation that we see at the surface, especially for this particular major fault zone. Seismic anisotropy and especially SKS splitting may become particularly important for resolving such questions.

C O N C L U S I O N S
We have presented new SKS splitting parameters for the Alpine region that have been robustly determined. The spatial resolution is unprecedented, and it includes in particular information on splitting delays that has been relatively uncertain so far. The high station density of the AlpArray refines the picture of seismic anisotropy in the Alpine region, and it gives us clues about the deformation in the subsurface and its causative effects. While the constraints agree well with previous results, for example the continuous rotation of fast orientations along the Western Alps and the lack of rotation in the Central Alps, there are new features that emerge, pointing at the role that the Alpine orogeny plays in the mantle deformation; yet it is just one of several geodynamic processes that are present in the area. Our study includes results obtained from the OBS measurements in the Ligurian Sea. Even though these results have higher uncertainties due to the shorter OBS deployment, they cover previously unexplored area and give new perspective connecting the anisotropy measurements onland. Within the Alps, we have identified two regions of particularly strong splitting, which are at the northern and the southern edge of the subducting lithosphere in the Central-Western Alps. The two zones are separated by a belt of much weaker anisotropy that has an intriguing relation with surface tectonics. This lays the focus on the nature of surface-mantle coupling, and on how deep processes drive deformation that we see at the surface. Questions like these indicate that seismic anisotropy and especially SKS splitting is a primary observational tool for understanding geodynamic processes in the upper mantle.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online.
In the online-only supporting information (OnlineOnlySupplement.zip), we provide text files with the obtained splitting parameters shown in Figs 3 and 7(b). These files contain station coordinates and names, and splitting delay times and fast axis orientations with their errors for 512 and 650 stations, respectively. We also add text files with similar information shown in Fig. 8(a) and (b), which are the interpolated versions of Figs 3 and 7(b). These files contain coordinates of the fine grid (roughly 5 km spacing) and values of splitting delay times and fast axis orientations interpolated onto this grid. Note, that in Fig. 8, the fast orientation bars are given on sparser grid to keep the figure legible. In addition, we also provide a file with coordinates used for cropping the region in Fig. 8.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

A P P E N D I X : S P L I T T I N G I N T H E P R E S E N C E O F N O I S E
We will investigate here first, what is the relation between the apparent and true splitting intensity in the presence of noise. Then we will show a numerical example and explain, how Wiener filter helps to correct for the bias.

A1 Bias in the splitting intensity method
The splitting intensity technique of Chevrot (2000) assumes the waveform model where T is a matrix which contains m transverse component records with n data points, r is the derivative of the radial component and s is the splitting vector (a set of individual splitting intensities), see eq. (A1) in the appendix of Chevrot (2000). In this model, noise (matrix N) affects only the transverse component. To see how eq. (A1) performs under more general noise conditions, we introduce 'noisy' radial and transverse components: where R 0 and T 0 denote the initial noise-free radial and transverse components with the noise N, N 1 and N 2 (independently constructed noise vectors). The Chevrot technique is used to estimate the splitting intensity from real measurements as where s * in our notation represents the splitting intensity biased by noise. Chevrot (2000) and Monteiller & Chevrot (2010) implicitly assume that this is equal to see, eq. (4) in Monteiller & Chevrot (2010). If we insert the noiseprone radial and transverse components into the initial eq. (A5), we get The mixed terms and 2R 0 (t)N (t) dt can be neglected, if we assume that N and N 2 are uncorrelated with T 0 and R 0 . Then we get The noise level can be defined by its mean μ and variance σ 2 . For N(t) we assume μ = 0. The last equation then becomes with the number of samples n. The relation between the noise-biased and unbiased splitting intensity then becomes what can be simplified by substituting to Eq. (A14) shows that the splitting intensity is biased in the presence of noise, and the apparent splitting intensity s * is expected to be smaller than the true value s, as x is always positive. The bias tends to increase with the noise level characterized by σ 2 . It affects only the splitting intensity, and through it, the splitting delay is affected as well. The nature of the bias in the splitting intensity is the effect on the radial component, which also needs to be considered. The other splitting parameter, the fast orientation, is not biased by the presence of random noise.

A2 Noise bias from a simple numerical example
We have examined the bias caused by noise in a synthetic test. The waveforms were calculated for a set of given input parameters (see Fig. A1, top-left panels) using Ricker wavelets (Ryan 1994). Since we wanted to investigate the behaviour of the calculated splitting intensity values under different noise conditions, band-limited additive white Gaussian noise was added to the synthetic north and east components (Fig. A1, bottom-left panels). These components were rotated into radial and transverse components and the splitting intensities were calculated using the initial eq. (A5). Grey dots in the top-right panel of Fig. A1 show splitting intensity measurements taken under different noise conditions, as a function of the noise variance σ 2 . Splitting intensity is not a unique function of σ 2 . It rather depends on the realization of the noise waveform. We try many different noise waveforms, and plot the splitting intensity as a function of σ 2 that corresponds to them. Some noise waveforms perform better (by chance) than others. We see that a highly likely outcome of adding random uncorrelated noise to the R component is a reduction in splitting intensity value. An increasing bias of splitting intensities to smaller values is clearly visible. The red cross in the upper plot and the red line in the histogram shows the reference splitting intensity, calculated from the noise-free case. It clarifies that the apparent splitting intensity values become biased with increasing noise level. The blue line describes the theoretical decay of the splitting intensities from eq. (A12), when the noise-free splitting intensity (s 0 ), the number of samples and the integral of squared noise-free radial derivative is inserted. Figure A1. Top-left panels show noise-free synthetic horizontal waveforms (Ricker wavelets) and their rotation into radial and transverse components. Noise was then added to these synthetics, see bottom-left panels. Repeating this for various noise conditions shows a decrease of splitting intensity with increasing noise level (top-right panel). Blue line indicates the expected decay of splitting intensity with respect to noise level. Bottom-right plot shows a histogram of biased splitting intensities; the true value is given by the red line. Figure A2. Histograms of splitting intensity in the biased (grey) and corrected (green) cases under various noise conditions. The energy correction shifts the splitting intensities back to the expected value.
If we now move on from a fixed-parameter scenario to a range of possible outcomes, we can look at the distributions from the biased and unbiased splitting intensities in Fig. A2. The variation of the input splitting delay time between 0-2.5 s leads to a respective shift in the expected splitting intensity values s 0 . We see that various noise conditions lead to an underestimation of the calculated splitting intensities. The effect is larger for larger input values. Figure A3. Effect of filtering: left-hand panel shows how average cosine similarities improve (rise) after Wiener filtering is applied (green as opposed to red). This leads to a positive shift in absolute values of average splitting intensities in the middle panel and to on average larger estimates of the splitting delay in the right-hand panel.

A3 Correcting the bias, Wiener filter
Eq. (A14) suggests a way to correct the bias in the splitting intensity measurements by estimating the noise level σ 2 for x. This is illustrated in Fig. A2. We note that the correction removes the bias on average. The maxima of the distributions are close to the expected value, while it increases the spread of the distribution.
Another way of addressing the issue of the noise bias is by attempting to eliminate it from the traces as good as possible. This is the approach chosen by Monteiller & Chevrot (2010), which we have implemented in our study as well. A Wiener filter (Robinson & Treitel 1967) solves this problem in a least-squares sense. The mean-squared error between the filtered input and the desired output is minimized, to find the best filter operator. This replaces the waveform on the radial component by a reference waveform, thus making traces as similar as possible. With an expected dominant period between 5-20 s for an SKS-wave (Monteiller & Chevrot 2010), we have convolved our initial waveforms with a Ricker wavelet (Ryan 1994) with the main period at around 15 s. This made the technique more robust against incoherent noise. Monteiller & Chevrot (2010) also suggest averaging over small azimuthal windows. This seems to be a reasonable option, when data quality is only moderate.
The splitting parameters in this study were determined with the Wiener-filtered splitting-intensity method. In applications where the noise is stronger, one may also remove the bias explicitly. The filtering increases the signal/noise ratio in the measured data. This gives rise to, on average, higher cosine similarity between radial derivative and transverse, see the left plot in Fig. A3. This also results in a shift towards on average greater absolute splitting intensities (middle panel in Fig. A3). This has two effects: (1) with respect to the final selection criterion of a minimum cosine similarity of 0.7, more SKS phases could be considered as 'split' and were used to constrain the splitting function, after the Wiener filter was applied, (2) the larger number of good measurements, as well as the on average greater absolute splitting intensities, leads to a more accurate fit and a general shift towards higher splitting delays. This positive shift in the splitting delays was expected from the synthetic tests we performed earlier. The distribution maximum of splitting delays is found now at values around 1.4 s (Fig. A3, right-hand panel).