Arrival angles of teleseismic fundamental mode Rayleigh waves across the AlpArray

The dense AlpArray network allows studying seismic wave propagation with high spatial resolution. Here we introduce an array approach to measure arrival angles of teleseismic Rayleigh waves. The approach combines the advantages of phase correlation as in the two-station method with array beamforming to obtain the phase-velocity vector. 20 earthquakes from the first two years of the AlpArray project are selected, and spatial patterns of arrival-angle deviations across the AlpArray are shown in maps, depending on period and earthquake location. The cause of these intriguing spatial patterns is discussed. A simple wave-propagation modelling example using an isolated anomaly and a Gaussian beam solution suggests that much of the complexity can be explained as a result of wave interference after passing a structural anomaly along the wave paths. This indicates that arrival-angle information constitutes useful additional information on the Earth structure, beyond what is currently used in inversions.

It has been observed frequently that seismic phases do not always arrive from the direction that is expected if the Earth had a simple 1-D structure. Arrival-angle deviations are particularly pronounced for surface waves; they are usually found to be up to ±15 • (Levshin et al. 1994;Laske et al. 1999) or ±10 • with extremes up to ±30 • (Laske 1995;Cotte et al. 2000;Maupin 2011). These deviations are not mere curiosities, but they can be important for constraining the 3-D structure of the Earth. The assumption about great-circle-parallel propagation biases surface wave tomography (Lay & Kanamori 1985). For that reason many studies have corrected the two-station measurements of phase velocities for the true arrival angles (Baumont et al. 2002;Pedersen 2006;Kolínský et al. 2011;Foster et al. 2014b;Chen et al. 2018). Some (e.g. Meier et al. 2004;Kästle et al. 2016;Soomro et al. 2016) have assumed that the arrival-angle deviations average down to zero after considering many measurements along the same station pair. Array-based methods (beamforming) are often used to overcome the problems of the two-station method. Phase velocity is measured simultaneously with the true arrival angles, meaning the slowness vector is obtained * www.alparray.ethz.ch (Cotte et al. 2000(Cotte et al. , 2002Pedersen et al. 2003;Maupin 2011;Foster et al. 2014a;Kolínský et al. 2014;Pedersen et al. 2015;Chen et al. 2018). It has also been suggested, based on a linearization of ray-theoretical equations (e.g. Woodhouse & Wong 1986), that the deviation from the great-circle azimuth is proportional to the transverse derivative of phase velocity (and that amplitude anomalies are sensitive to the second lateral derivative). This approach was used by Laske (1995) and Yoshizawa et al. (1999) for inversion of polarization measurements into phase maps. The additional information can improve global phase velocity maps (Laske & Masters 1996). The approach was generalized by Larson et al. (1998) for anisotropic structures.

Spatial variations of arrival-angle deviations
Over the last years, data from large arrays have become available to study how arrival-angle deviations vary in space: Pollitz (2008) has studied the wavefield in western US using data from large USArray project observing band-like patterns of amplitude anomalies and concluding that most of the pattern is part of the incident wavefield. A similar conclusion was made by Lin et al. (2012). Foster et al. (2014a) have again shown evidence of intricate spatial patterns of arrival-angle deviations across the USArray and they tried to model the observed distribution, but only with mixed success. They suggested that arrival angles are very sensitive to the source location. Pedersen et al. (2015) studied the LAPNET/POLENET network, and they also discussed the results with respect to varying earthquake location. Chen et al. (2018) found spatial patterns of arrival-angle deviations that resemble the stripe-like pattern observed by Foster et al. (2014a). They also discuss the distinct arrival-angle deviations measured at the same place for different events. With reference to Pedersen et al. (2015), they attribute this observation to multipathing due to complicated structure.

Complex wave propagation
Surface wave propagation and scattering is a complex phenomenon. This complexity has been approached by a variety of techniques, each of which involves specific approximations. Yet, attempts to explain the actual cause of arrival-angle (phase and amplitude) deviations are rare. Usually, the deviated arrival angles are attributed to simple off-great-circle propagation, meaning that the waves are avoiding the low-velocity regions, over broad range of periods (Lay & Kanamori 1985;Alsina et al. 1993;Levshin et al. 1994), to multipathing when late arrivals of waves are split from the main wavegroup and propagate along longer path (Maupin 2011;Chen et al. 2018) or to diffraction and scattering (Friedrich 1998;Maupin 2001;Bodin & Maupin 2008). Pedersen et al. (2015) suggest to explain the different arrival angles at the same array from sources located close to each other by scattering and diffraction near to the source. A similar observation is discussed also by Maupin (2011) and Foster et al. (2014a).
The wavefield gathers its complexity along the path before coming to the network. To account for this, superposition of many plane waves (Stange & Friedrich 1993;Friedrich et al. 1994;Friedrich & Wieland 1995) was used to fit jointly phase and amplitude at stations in Southern Germany showing that most of the large anomalies of the wavefield were imported from outside of the network, as had been predicted by Wielandt (1993). Friedrich (1998) came to a frustrating conclusion that any data observed by a network can be explained without allowing any structure under the network at all attributing all the anomalies only to the incoming wavefield. Similarly, the shape of non-planar incoming wavefield was inverted simultaneously with local phase velocity field by Bruneton et al. (2002), Pollitz (2008) and Salaün et al. (2012). The bias resulting from structures outside of the tomographic grid was removed using the jackknife inversion by Chevrot & Zhao (2007). Two-plane wave approach (Li et al. 2002;Forsyth & Li 2005) was used to account for effects of multipathing interference that occur between the source and the array and applied in northeastern United States and southeastern Canada (Li et al. 2003).
Single scattering theory for complex surface wave propagation was developed by Snieder (1986) pointing out that Rayleigh waves are mainly forward-scattered. A multiple scattering approach was proposed and numerically tested by Friedrich et al. (1993). Pollitz (1994) introduced a calculation scheme for phase perturbations caused by small (few wavelengths) scatterers. Spetzler et al. (2001Spetzler et al. ( , 2002 discussed the limits of the ray-theory-based studies showing that scattering is needed to obtain higher resolution of the Earth's structure and that increasing the data coverage is not enough. However, Sieminski et al. (2004) argued that for regionalscale tomography, we can counterbalance the shortcomings of the ray theory in considering the finite-frequency effects by denser ray coverage.
A complex geometry of heterogeneities was studied by Maupin (2001) including mode-coupling and anisotropy. Finite-frequency effects were incorporated by Ritzwoller et al. (2002) by introducing sensitivity kernels to global 'diffraction' tomography taking into account the scattering over the first Fresnel zone. To explain the effects of off-great-circle propagation, finite-frequency effects were considered by Yoshizawa & Kennett (2004) considering the 'influence zone' (1/3 of the first Fresnel zone, Yoshizawa & Kennett 2002). The same technique was applied by Isse et al. (2006) to the Philippine Sea region. Zhou et al. (2004) introduced 3-D sensitivity kernels for phase, amplitude, and arrival angles showing a strong influence of offgreat-circle structure to all these three observables. Global upper mantle structure was determined using that approach by Zhou et al. (2006) resolving small-scale structures taking into account wavefront healing and radial anisotropy. Chevrot & Zhao (2007) used 3-D sensitivity kernels to invert for the structure of the Kaapvaal craton. Amplitude and phase data were inverted using 2-D sensitivity kernels and representing the incoming wavefield by two-planewave approach of Yang & Forsyth (2006a) for Southern California by Yang & Forsyth (2006b) and by Yang & Ritzwoller (2008) for the western United States (USArray). Finite-frequency sensitivity kernels were calculated by de Vos et al. (2013) for the two-station method. As the above referenced studies, our work also focuses on a particular (yet different) characteristic of complex surface wavefield propagation.

AlpArray
Here we propose to observe in depth the arrival-angle deviations from earthquakes distributed in different directions, and to find a general explanation for the arrival-angle deviation patterns mapped over the large region covered by the AlpArray network. Fig. 1 shows broad-band stations available for this study. These are temporary stations (full green triangles) deployed within the AlpArray project inside the region defined as 250 km distance from a smoothed 800 m contour line around the Alps (yellow line, Hetényi et al. 2018). In addition, we use also permanent stations (blue triangles) both inside and outside of the Alpine region across the European continent in the area limited by 39 • N-51 • N and 3 • W-25 • E. The AlpArray project encompasses also ocean bottom seismometers (empty green triangles) which were recording during 2017 but the data were not available at the moment our study was performed. Details about installation and performance of part of the AlpArray network can be found in Fuchs et al. (2015Fuchs et al. ( , 2016. The mean station distance is 30-40 km, see Hetényi et al. (2018), which represents a finer spatial sampling than in the case of USArray (Foster et al. 2014a).

Our study
In this paper, we further develop the beamforming approach of Kolínský et al. (2014) who used single-array phase-velocity measurements to obtain a local 1-D velocity model beneath the array. We use fewer earthquakes than Pedersen et al. (2015); however, we study the patterns of arrival-angle deviations over a large area covered by the AlpArray network, using the floating subarray (miniarray) concept. Any subset of stations of a large network can turn into a local array. This allows us to map the arrival-angle deviations in space at the receiver side and draw conclusions about their origin along the propagation from the source. All the abovementioned kinds of arrival-angle deviations, previously explained as off-great circle, multipathing, diffraction, scattering and closesource effects, appear in our observations as well. Our paper focuses  on arrival angles exclusively, for the sake of clarity. Phase velocities will be studied later, once the wave-propagation phenomena are well-understood. We compare the observed arrival-angle deviations with modelling, in an approach based on the model of Nolet & Dahlen (2000). The above-mentioned studies introducing phase and/or arrival-angle sensitivity kernels (Ritzwoller et al. 2002;Zhou et al. 2004;Yang & Forsyth 2006a) or predicting phases and/or arrival angles directly (Pollitz 1994;Maupin 2001) could in principle be used for our purposes as well. The Nolet & Dahlen (2000) approach is attractive for its simplicity (as well as computational efficiency). This will allow us to propose a different and actually simpler explanation for the observed effects than the latter papers.

M E T H O D 1 : A R R AY M E A S U R E M E N T
We provide an array measurement of surface wave phase velocities, adopting the method proposed by Kolínský et al. (2014). It takes advantage of two traditional approaches. The first one is the classical two-station method of measuring phase velocities between pairs of stations in the time domain as described by Kolínský et al. (2011). The second approach is array beamforming for determining the phase velocity and the true arrival angles of propagation for a given array. Selection of the fundamental mode by filtering in the frequency domain, tapering in the time domain, and estimation of the time delay between a pair of stations is similar to the two-station method. We however overcome the main problem of the two-station method, by using the beamforming afterwards. So, we do not need to assume anything about the propagation paths. We do not need to know the origin time of the event, and neither the location of the event. Any array design has a limited range of wavenumbers that can be resolved, and so beamforming suffers from potential biases if more signals are present in the record that we cannot resolve. We alleviate this by beamforming a signal which had already been separated before by filtering and tapering. This also allows us to measure the phase velocities at very low wavenumbers which are below the threshold resolvable by a given array when applying the beamforming in a traditional way. Our method consists of following steps.

Selecting the fundamental mode
We start by processing each of the records separately station by station. We follow the procedure described by Kolínský & Brokešová (2007); using the classical method of Fourier transform-based multiple filtering to analyse the dispersive records (Dziewonski et al. 1969). We use non-constant relative resolution filtering applying around 80-100 narrowband Gaussian filters in the frequency domain. The filters are wider towards higher frequencies. The width of the filters does not only depend on frequency, but it is further adjusted according to the properties of the given signal so that optimal resolution is achieved both in frequency and time domains (Kolínský 2004). The filters are hence different from event to event. However, they are kept the same for all the records from a single event. Each of the filtered spectra is then transformed back to the time domain. As a result, we obtain a set of quasi-monochromatic signals, see Fig. 2, orange lines. The frequency of these signals is not constant for two reasons: first, applying the symmetric Gaussian filters to asymmetric spectra results in a systematic shift of the dominant frequency. Second, the filters have finite width and all the frequencies of the filtered spectra are contained in the resulting time-domain signals. To avoid any systematic error, we do not use the central frequencies of the filters in further analysis. Instead, we calculate instantaneous frequencies for every time sample using the analytical signal corresponding to every quasi-monochromatic signal (Levshin et al. 1989). Envelopes of these signals are calculated as a modulus of the analytical signals.
The group-velocity dispersion curve is then determined. We pick all local maxima of all envelopes. The fundamental mode is wellrecognized at long periods above 60 s. Usually, there is only one envelope maximum at this period range in the given time window (1.5-2.0 hours). We start the picking at this maximum moving to shorter waves and picking always the maximum closest in time to the previous maximum at longer wave signal. We select the groupvelocity dispersion curve, which does not depend on amplitude-for shorter waves, our picked maxima are not necessarily the highest ones at the given filters. For this group velocity dispersion curve, epicentre coordinates and origin time are used. The group velocity dispersion is affected by measurement errors. However, we do not need the dispersion curve itself for further analysis. Instead, we taper the quasi-monochromatic signals in the time domain at a predefined length with respect to this dispersion curve. We use a window of 4 periods generally (see Fig. 3). We keep 0.8x the period around the maximum at the same amplitude. 1.6x the period to lower and higher times is smoothed to zero by a cosine window. Fig. 3 shows examples of quasi-monochromatic signals at 50 and 100 s. For the 100 s wave, there is only one wavegroup, which is tapered and smoothed. For the 50 s wave, we see more wavegroups. The one corresponding to the fundamental mode is not of the highest amplitude. However, using the method of continuous ridge picking, we can find the proper group and we taper it as well. We follow that procedure for all the filtered signals, see Fig. 2, blue lines. As well as filter widths in the frequency domain, also the taper lengths in the time domain are kept the same for all records of a given event. Relating the tapering window to the number of periods that results in longer time-domain signals for long periods, can be understood also as a transformation of the filters in the frequency domain, which are narrower for low frequencies. The width of the filters in the frequency domain is set so that the instantaneous frequency of the tapered wavegroup is nearly constant. Because in the next steps we use only these four-period tapered signals for cross-correlation, the error of the group velocity estimation does not influence the result. The same applies for the possible error of origin time and source location-for further calculation, only the tapered signals are used and any information about the event is not needed, meaning, any potential source parameter error does not influence the result.

Subarray measurement
The measurement of time differences between quasimonochromatic signals is carried out according to the method described by Kolínský et al. (2011Kolínský et al. ( , 2014. Filtered and tapered quasi-monochromatic signals are shifted in time sample by sample and their correlation is computed for each time shift. One of the stations is set as a reference, and time shifts corresponding to the highest correlation are determined consecutively for the other stations. At this stage, we use only the tapered fundamental mode without any assumption about the propagation direction. We do not struggle with the unknown cycle number for shorter waves as it was the case in the study by Kolínský et al. (2014). Our arrays are so small (half-aperture = 80 km) that even the shortest periods of 25 s are long enough (∼90-100 km) so that we observe only one correlation maximum (time delay) in the preset time window of correlation shift. We set this time window according to the half-aperture of our array so that in case the wave propagates right along the longest distance between two stations, the time window is long enough to cover a delay of a priori assumed (conservatively low) phase velocity. Usually, we use a time window of ±35 s (lowest possible phase velocity set to 2.3 km/s for the half-aperture of 80 km). In case the wave travels in any other direction, the resulting time delay is always smaller than 35 s.
The only exception is the regional event from Greece. For that close earthquake, we extended the range of periods down to 10 s. For such short waves (λ ∼ 35-40 km), two correlation maxima are sometimes obtained for the given correlation time window of ±35 s. We again use the same procedure as in the case of the group velocity dispersion curve above: we follow the continuous dispersion curve (time delays) starting at long periods towards shorter periods and the proper time delay is selected also for these short waves.
Then, we determine the true velocities and directions of surface wave propagation. Assuming the wave is planar, knowing the array geometry (interstation distance vector r i for i-th station with respect to the central one) and having the measured time differences t i , we have a set of equations for two unknown components of the slowness vector s for each period k: see also Kolínský et al. (2014). Fig. 4 shows the procedure. Plot (a) is a map view of a subarray centred at station A008A and plotted in Cartesian coordinates. For example, we use a wave of 70 s period coming from the Kamchatka earthquake. In such a case of an array of 15 stations, 14 time differences t i=1..14 are determined with respect to the selected reference station A008A and hence 14 equations for two unknowns are available. Linear regression is used to determine the slowness vector for each period. We fit a plane to all the time measurements. Coloured lines in Fig. 4(a) show the continuation of the wavefront across the subarray, with colour indicating time (as for the times at the stations). Time differences range from −15.6 s for NNE station A084A to +20.0 s for station A018A in SSW. Station A005A has a significant residual being off the plane and we removed it for the final calculation. The slope of the intersection of the regression plane with the time-X plane gives the x component of the slowness vector (s x = −0.08 s/km, slope measured with respect to X axis). The slope of the intersection of the regression plane with the time-Y plane gives us the y component of the slowness vector (s y = −0.24 s/km, slope measured with respect to Y axis). The intersection of the regression plane with the X-Y plane gives the wavefront direction. As a result, we have both the absolute value of the phase velocity as well as its direction, depending on period ( Fig. 4 shows one period only). The direction is measured clockwise between the north and the ray, which is perpendicular to the wavefront (green line with arrow). For the given example, the phase velocity is c = 1/ s 2 x + s 2 y = 3.99 km/s  and the true arrival angle is 18.5 • (after A005A has been removed).
Calculating the above procedure for the whole period range, we obtain a local phase-velocity dispersion curve, where the velocity itself is already corrected for the true arrival angles, which vary with period. The procedure makes use of the assumption of plane wave propagation. That assumption is well-satisfied by all event distances and period in this study at the aperture of our subarrays.

Subarray design
To image the wave propagation across wider region, we used a floating subarray design. Any station of the network is considered as a centre of a subarray with uniformly defined properties. After testing many subarray sizes, we chose an aperture (diameter) of 160 km. This size is small enough to map local variations of true arrival angles and wide enough to resolve long-period waves up to 170 s, which is the limit we are able to measure consistently using the AlpArray network taking into account various sensor types. Hence, for each station, we look for neighbouring stations in a distance less than 80 km around. To avoid the measurement of very small time delays, we exclude stations closer than 20 km. Station pairs with interstation distances between 20 and 80 km are then used to calculate the time delays with respect to the central station. Array beamforming is provided, local dispersion curves as well as true arrival angles are computed and the whole procedure is repeated for all stations in the network. To provide reliable measurements, we only consider the station to be a centre of a subarray if at least five other stations are found in the aperture around (minimum subarray size = 6). If there are less of them, the beamforming is skipped. However, such a station can still be used as a neighbouring station when some other of its neighbours turns into a subarray centre.
As an example, Fig. 5 summarizes the results for three events and one subarray centred at station A008A. We represent the results as a palette of vectors, where colour represents the period, length the phase velocity and direction the true arrival angle. Stations in the 82 km radius aperture are shown by white triangles. The selection differs from event to event. We have 15 stations for Kamchatka (as in Fig. 4) and 16 stations for Sumatra and Ascension earthquakes. Stations A002A and A005A were used for the Kamchatka and Sumatra events (A005A was later removed from the Kamchatka event as described earlier), and stations A002B and A005B for the Ascension earthquake. The colour scale of periods is the same for all events; however, the ranges for individual earthquakes are different (42-160 s for Kamchatka, 70-160 s for Sumatra and 24-150 s for Ascension). We see velocities ranging from 3.7 to 4.4 km/s (local dispersion curve). The spread of true arrival angles is around 10 • peak-to-peak for all the events. For comparison, the geometrical great-circle paths are plotted by a solid white line, see also the white line in Fig. 4(a) for the Kamchatka example. For Kamchatka and Sumatra events, the true arrival angles are systematically deviated to lower values. For the Ascension event, the true arrival-angle range is distributed around the geometrical value. Other stations of the network are shown by green triangles. Obviously, the number of stations of each subarray differs because of the network geometry. The geometry of a given subarray may also differ event by event because of variable availability of records, as shown in the example. However, the subarray design criteria (distance range of 20-80 and minimum size of 6) are kept the same for all the events used in our study.
The floating subarray approach allows mapping the wave propagation in different ways. The absolute value of the phase velocity can be assigned to the location of each subarray centre. Keeping the location and moving over the period range, we have a local dispersion curve. Keeping the period constant and moving over the  whole network produces a phase-velocity map. Evaluating these two latter observables will be a topic of next studies. Similarly, we can check the period dependence of the true arrival angle for each subarray having 'azimuthal dispersion'. Following the same period over the network, we can plot the true angles of propagation as a map of vectors pointing in the direction the wave is coming from (as in Fig. 5). Turning these vectors (e.g. Fig. 5) by 90 • , we can trace the wavefronts, as the wavefront is perpendicular to the vector of phase velocity. To keep the approach universal, instead of the absolute arrival angles measured from north, we will rather inspect the deviations of the true propagation angles with respect to the great-circle paths, both given as a clockwise angle from the north. Fig. 4 shows the meaning of the arrival-angle deviation: it is the difference between the measured arrival angle of the ray (18.5 • in the example) and the geometrical great circle backazimuth (21.8 • in the example) yielding the arrival angle deviation of −3.3 • for the 70 s wave propagating from the Kamchatka event. Maps of these deviations are the matter of our study. We will study how these deviations evolve in space, over the period range and how they differ event to event.

Time residuals
To evaluate the fit of a plane wave to the data from a subarray, we calculate a cascade of time residuals. A time residual is the difference (L1 norm) between the measured time delay at a given station with respect to the reference one and the time delay predicted by the plane-wave fitting. We obtain these residuals for each event-subarray-station-period combination. We can then average these time residuals over the neighbouring stations to get one mean residual characterizing the fit of the plane wave of the given period for the whole subarray. We can as well average these time residuals over the period range for evaluating each single neighbouring station (individual station residual). Last, we can average over both the periods and neighbours to obtain a single mean time residual (total residual) characterizing the whole subarray. An example of a distribution of residuals is given in Fig. 4(b), where a histogram of mean (neighbours-averaged) residuals of all 429 subarrays at period of 70 s for the Kamchatka earthquake is given. They peak between 0.1 and 0.3 s, which is also the case for the example subarray A008A in plot (a), where the residual is 0.127 s. The use of these residuals is described later in the text, and in detail in the Appendix.

M E T H O D 2 : S I N G L E -S TAT I O N M E A S U R E M E N T
To complement our observation of wavefronts measured by the floating subarray technique, we can also trace the wavefronts by single-station measurement. For this purpose, we pick the zerocrossings of the ground velocity amplitude, see Fig. 3, red lines. Times of these zero-crossings cannot be used to obtain the phase velocity in any sense -phase velocity cannot be measured by a single station not knowing the phase at the source. However, for our purposes, one can pick the same phase at all the stations in the network. We are interested in the evolution of the wavefronts across the network and we do not need the absolute values of these time picks. Usually, we pick three zero-crossings of the selected fundamental mode wavegroup under the tapered envelope. Ideally, they are exactly a period apart (if picking only negative-to-positive crossings as in Fig. 3). Assuming continuity of wavefronts across the whole network, we correct the measurements for the cycle skips selecting any of the three picks so that the times are continuously evolving in space. To visualize this single-station measurement over the network, we will later plot these zero-crossing time values as a contour map with lines for every 10 s of wave propagation for each period. These contour lines can be directly compared with the wavefront vectors obtained from the subarray measurement described above. This single station method is used only to qualitatively benchmark the results of our subarray measurements. All the results of this paper are based solely on arrival-angle deviations calculated using subarrays. Subarray allows us to directly measure the arrival-angle deviations together with estimating the quality of the measurement (residuals). This is not possible using the single-station phase picking since when producing the contours (wavefronts), some spatial smoothing needs to be applied and estimating the quality of the measurement is difficult.

DATA
We selected 20 earthquakes spanning the two first years of the AlpArray project (2016 and 2017). The selection consisted of the following steps: we listed all the M > 7.0 events for that time period from all over the world. We added M > 6.5 events from missing regions (Indonesia, Tajikistan), and added M > 5.5 events for closer region around Europe (Greece). This resulted in 58 events in total. We removed events with big aftershocks (M > 5.0) following closer than 1 hour after the main shock so that we did not have overlaps of the records. We removed events from the same location except for four pairs, which we kept for comparison. We checked the records removing events with poor surface wave amplitudes and we made sure our list covered roughly the whole azimuthal range including events from rare azimuths even when having poorer SNR (Bouvet Island). Fig. 6 shows the final set of 20 earthquakes including greatcircle lines pointing to the centre of the AlpArray region (47 • N, 10 • E). Table 1 lists the events including location, date, origin time, magnitude and depth. Fig. 6 shows also distances to the station A037A which is the closest station to the 47 • N and 10 • E. Letters A, B, C and D label four pairs of events at similar location. The distances (green), AlpArray region of interest (yellow) and detailed map boundary (magenta) used for next figures are shown in Fig. 6 as well. The biggest event in the given time period was the Mexico 2017-09-08 M = 8.1 earthquake. The basic idea of the AlpArray project was to have simultaneous recording across the Alpine region. However, some stations were being added still during the 2016 up to mid 2017. At the beginning, the region was therefore not covered uniformly and the western part was missing. For each earthquake, we used all available temporary AlpArray stations as well as all permanent stations available for the European continent within the range of 39 • N-51 • N and 3 • W-25 • E.  Our station list of all potentially available broad-band stations contained 948 items. Considering only the stations which were available at least for one of the events, we dealt with 890 unique stations after the download. All these 890 stations are plotted in Fig. 1. However, much less stations were downloaded for each particular event, see Table 1. The range spans from 654 to 804 records and is 749 per event on average. The events in the Table 1 are listed by date. We see increasing record numbers from the first event in January 2016 (Alaska) peaking in June 2017 (Guatemala). This corresponds to temporary stations being added to the AlpArray backbone. Since July 2017 (Komandorskiy Island), the number decreases again. This is because not all the temporary stations are online and so for the most recent events, data from some of them were not available at the moment since they are collected in periodical campaigns in the field occasionally. Fig. 1 shows that in addition to the AlpArray region depicted by the yellow line, we have three other localities where station coverage allows for subarray measurement: the Apennine Peninsula, the Vogtland region and the Pyrenees. We excluded the Pyrenees when plotting the arrival angle deviation maps because the subarrays are sparse and do not map the region continuously. The Apennine Peninsula and the Vogtland region are kept in our array analysis. The station coverage allows for the overlapping subarray measurement complementing the results in the Alpine region. Ocean bottom records were not available at the time when our study was conducted. For the actual measurement, however, many records needed to be removed because of the poor data quality. Procedures for assuring good data quality are discussed extensively in the Appendix.

General remarks
In the next paragraphs, we describe our observations both for particular events as well as general features which are common for all or several events together. Looking at Table 1, we see that the average size of the subarrays is fairly uniform, while the area covered by all subarrays altogether increases with the continued AlpArray deployment during 2016 and 2017. The highest station density is in the Apennine Peninsula thanks to the Italian permanent networks; the minimum subarray size was set to 6 for all events. Period ranges are also given in Table 1. They differ from event to event due to the  (nghbrs), stations used solely for the single station measurement (sngls) and number of stations not used at all (not) are given for each earthquake. Average (avrg) and maximum (max) sizes of the subarrays as well as the period range are shown. Eight earthquakes form event pairs, labelled as A-D (see Fig. 6 and the text).
The last column refer to the figure number where the period-dependent arrival-angle deviations are shown ( Fig. 11 and Supporting Information Figs S1-S19). different magnitudes, different distances to the network and the conditions along the ray paths. During the processing, we usually started with a broader range. After preliminary results were obtained, we narrowed the range for periods, where robust results were obtained. This means that even for events of comparable size and comparable distance, we may generally use different period range for the analysis. Fig. 7 shows the most comprehensive way of presenting our results from array measurements. At the top, we see the phase-velocity vectors for 550 subarrays for the broad period range of 30-170 s obtained from the records of the Mexico event. Plotting the vectors follows the same scheme as in Fig. 5. Besides the subarrays, we plotted also the stations used as neighbouring stations for subarray calculation but not used as a centres of any subarray (31 stations, blue triangles) as well as stations used exclusively for the single-station phase time picking (75 stations, light blue triangles). Phase time contours are not shown in the figure since the figure encompasses the whole period range and the phase time contour map can only be presented for single selected period (see Fig. 8). In addition to this, we also plotted the stations for which the data were downloaded for this particular event, however, due to the poor quality we could not use the records (137 stations, red triangles). Since we plot only the subarray results not smoothing any quantity over the map, we kept also the Pyrenees region (four subarrays) in the map. The map at the bottom shows a detail of the central AlpArray region marked by the magenta rectangle in the top plot. We added the station names followed by the number of stations used for each particular subarray.

Full period range overview
In both plots, results with mean residuals < 2.0 s are plotted in the map. It means, not all the periods are plotted for all the stationssee, for example, the region at the Slovenia-Croatia border, where the shortest waves are missing. Looking at the colour distribution of all the vectors, we see the shorter periods (30-40 s, dark blue vectors) are generally coming from lower backazimuths than the mid-range periods (50-80 s, light blue vectors). For the northern part of the network, the longer periods (above 100 s, yellow and red vectors) are deviated further to higher backazimuths. In France and Germany, the spread of the arrival angles with period is rather narrow. In the Apennine peninsula and Slovenia, in contrary, the range of the arrival angles is rather significant. Following a single colour (particular period), we see smoothly changing patterns across the network -subarrays close to each other give similar results. Such a figure gives us an overview of the measurement, however, due to its complexity, it is difficult to emphasize the details. It is also not easy to compare the measured arrival angles with the geometrical ones. In the following figures, we separate the plots by period and rather plot directly the deviations, meaning the difference between the measured arrival angles and geometrical backazimuths (magenta lines in the bottom plot of Fig. 7) for each subarray (see also Fig. 4 for the definition of arrival-angle deviation). We show the results of array measurement (Method 1) in comparison with single-station phase-time contours (Method 2). We plot only a selected period for two events: 100 s for Sumatra (two maps on the left) and 60 s for Mexico (two maps on the right, compare  it with Fig. 7). For the maps on the top, we again plot the absolute arrival angles; however, we now rotated the vectors from the array measurement by 90 • to emphasize the wavefronts (red lines) instead of the ray paths. Such a representation allows for direct comparison with the phase time contours from the single-station measurement (yellow lines plotted every 10 s of propagation). To compare both measurements with the geometrical great-circle propagation, we plot the great-circle wavefronts by magenta lines in the background. We see a good match of the two observations for both events. Both the array wavefronts and the phase time contours from single-station measurement deviate from the great circles in the same direction and approximately by the same amount. For the Sumatra event, the true arrival angles are smaller in the northern part (waves propagate north of the great circles) and higher in the southern part (Apennine Peninsula, waves propagate south of the great-circles) of the network being almost parallel in the centre and towards the west. One might say that the rays are converging toward the AlpArray region from both sides. For the Mexico events, we see exactly the opposite feature: the rays are diverging having lower true arrival angles in the northern part (waves propagating south of the great circles) and higher in the southern part (waves propagating north of the great circles). Two maps in the bottom show the same array observation as the upper two maps but plotted in terms of arrival-angle deviations. We plot the difference between the geometrical ray paths and the measured phase-velocity vectors (see the detailed map at the bottom of Fig. 7, for the explanation of the deviations, see Fig. 4) or, the same, the difference between the measured wavefronts and the great circle wavefronts (upper two maps of Fig. 8). The colour scales of the two lower maps are the same and they are also kept the same for all the following figures in the present paper. We also keep the positive values corresponding to clockwise deviations as in Maupin (2011);Foster et al. (2014a) and Pedersen et al. (2015); see the difference between green and white lines in Fig. 4 as an example. Both maps for the Sumatra and Mexico events have negative deviations in the north and positive in the south; however, the waves propagate from almost the opposite directions, see the arrows and white dotted lines showing the direction of wave propagation and geometrical ray paths (spaced approximately a wavelength apart) and compare it with Fig. 6. Representing the subarray measurements by plotting the arrival-angle deviations allows for discussing the details of the wave propagation for each particular event as well as a mutual comparison between the events. It also allows seeing directly the size of the deviations, which, in both cases in Fig. 8, reach ±15 • . In the next figures, only the results of subarray measurements are shown and discussed. Single-station phase picking gives us a useful overview of the general behaviour; however, we cannot obtain precise quantitative estimate of the arrival-angle deviations. Fig. 9 shows arrival-angle deviations measured at the period of 50 s for six events. The upper pair of results is for earthquakes located very close to each other (Aleutian and Komandorskiy Island, see Fig. 6). The areas covered by the two measurements are slightly different but we note a strong similarity between the resulting maps. The deviations appear as north-south elongated stripes with changing polarity having the same position and comparable magnitude for both events. The two maps in the middle show the deviations for earthquakes located in the same direction, however, with different epicentral distance being more than doubled for Indonesia compared to Tajikistan (see Table 1 and Fig. 6). For the Tajikistan event, we see very narrow stripes of deviations switching between positive and negative values spanning the whole range of ±20 • . For the Indonesia event, the situation is different. All the deviations are negative, reaching over −30 • in some places. That is why we have chosen to show that map corrected by −12 • , meaning, we added 12 • to all the deviations. After this modification, the stripes are elongated roughly along the geometrical backazimuths. The bottom pair of maps shows the results for other locations. The Ascension Island earthquake produces consistent stripes somewhat tilted from the propagation direction spanning a broad range of ±12 • .The Bouvet Island earthquake gives us one dominant stripe of positive deviation followed by narrow negative and again smaller but broader positive deviation in the western part. For all the events, the deviation stripes are mutually fairly parallel and they are elongated almost along the great circle direction for the two upper events and with particular angle for the four lower events. White dotted lines showing the great circles are approximately a wavelength apart. Fig. 10 follows the same layout as Fig. 9 for a period of 90 s and four events. Again, the upper two maps represent results for measurement of deviations for two events at the same location (Ecuador 2016-04-16 M = 7.8 and Ecuador 2016-05-19 M = 6.9). Again, both results are similar. The deviations for the bigger event (left plot) are smoother and exhibit continuous stripes while those for the smaller event (right plot) are more scattered. However, the overall picture is the same. Lower panels show deviations for another two events. For the Mariana earthquake, we see mainly negative deviations across the whole region with stripes being only slightly pronounced. For the Chile event, we obtained deviations elongated almost along the great circles with pattern changing smoothly from positive over negative back to positive deviations across the whole region from NW to SE.

Period dependence
The period-dependence of the arrival-angle deviations is shown in Fig. 11 for the South Atlantic 2016-08-19 M = 7.4 earthquake. We plot eight maps covering the period range from 50 to 125 s. Note that for each of the maps, the number of subarrays is generally different as we keep a constant threshold of the mean residual < 2.0s for plotting the deviation maps. The stations used for respective maps (periods) are shown by magenta triangles. Subarrays with higher mean residuals for given period are not used for respective maps. Looking at each of the maps separately, we again see smooth stripelike patterns of deviations changing between positive and negative values elongated in NNE-SSW direction. Comparing the maps, we see that the patterns are similar (with an exception of the shortest 50 s period), but they are moving from east to west with increasing period. White dotted lines are kept at the same position for all eight maps (spaced approximately a 50 s wavelength apart) to allow the comparison. Starting at the 60 s map, there is a positive deviation stripe located approximately along a line connecting 15 • E in the southern and 18 • E in the northern part of the map. Following that stripe for longer waves, it is already shifted to between 13 • E (in the south) and 17 • (in the north) for 80 s and further to the west being shifted to 10 • E (in the south) and 11 • (in the north) for the 125 s map. The stripe is not only shifting to the west, but it is also changing its angle with respect to the geometrical backazimuths. For the shorter waves, it is more aligned to the direction of the geometrical ray paths, however, for longer waves, it is more deflected off the ray paths being almost N-S oriented. The same applies also for the other stripes of negative and other positive deviations. As the pattern shifts to the west, some stripes are moving out of the AlpArray, see for example the negative deviation being almost in the middle of the area at 70 s and only touching the western margin of the measurement for 125 s. Another stripe of positive deviation is, on the other hand, moving into the picture from east. At 90 s, the east margin shows still a negative arrival-angle deviation. At 100 s, it starts to turn into warm colours and for 125 s, a consistent stripe of positive deviation clearly shows up in the east.
Period-dependent deviation maps for the other 19 events we investigated are given as Supporting Information following the same scheme as Fig. 11 and labelled as Figs S1-S19, see also the last column of Table 1.
Looking at the deviation maps of different periods, we see that not only the stripes are moving laterally over the region, but also the distance between the stripes (or their width) is varying. Generally, the stripes are getting broader for longer waves.  Fig. S6) or the negative deviation at the west of the region measured for the Alaska event, Supporting Information Fig. S1, which stays at the same width from 40 to 82 s period. Looking at Figs 9 and 10 where the deviations for 50 and 90 s are plotted for different events, we see that the width of the stripes is not proportional to the period (wavelength). The Tajikistan and Indonesia events show much narrower stripes than the other four events for 50 s in Fig. 9 and as well, the two Ecuador events in Fig. 10 again show much narrower stripes then the other two events for 90 s period.

Comparison of event pairs
The observational results include arrival-angle deviation maps from four event pairs which are from almost the same location. Pair A (see Table 1), the Aleutian and Komandorskiy events are compared at 50 s period, see Fig. 9. For the whole period range results, see Supporting Information Figs S15 and S17. The distance differs only by 80 km (less than 1%) and the epicentres are 140 km apart. The magnitudes differ a lot (6.8 and 7.7), however, we do not see any differences in the pattern which could be assigned to the difference in magnitude of the events. Although the number of subarrays differs a bit, the deviation pattern is fairly the same. Pair B consists of the two Ecuador events (2016-04-16 and 2016-05-18), see Fig. 10 for 90 s period and Supporting Information Figs S4 and S5 for the whole period range results. In this case, the epicentres are even closer, differing only by 30 km (0.3%) in distance and by 40 km in location. The difference in magnitudes is similar to the previous pair (6.9 and 7.8). The bigger event shows more continuous and smooth deviation stripes while the smaller one suffers from more uneven results and the deviations are of higher values. This can be contributed to a lower signal-to-noise ratio and hence less reliable measurement of the smaller event although the residual distribution is very similar. The signal quality already plays a role in case of such a distant event while it did not affect the Aleutian event compared to Komandorskiy, since these latter were thousand kilometres closer. However, despite that, we still see the same pattern of deviations also for both Ecuador earthquakes. Pair C, the Guatemala (Supporting Information Fig.  S16) and Mexico (Supporting Information Fig. S18) earthquakes differ in distance by 120 km (1.2%). The difference in location is 180 km and the magnitudes differ even more than in the previous two cases (6.9 and 8.1). We see the same patterns for the period range of 40-110 s for both events again. Pair D consists of the two South Atlantic events (Figs 11 and Supporting Information Fig.  S6) of almost the same distance (difference of 47 km, 0.4%) but quite a bit different location (330 km). Magnitudes are comparable (7.2 and 7.4). Again, the patterns can be considered as being the same showing, however, a significant west-east shift. While the M = 7.4 event located more to the west resulted in the pattern shifted towards the east of the imaged region, the M = 7.2 event located east of the first one in South Atlantic resulted in similar pattern shifted to the west of the AlpArray region. The similarity of both patterns suggests that waves coming from both earthquakes pass the same structural feature. The flipped west-east shift of the pattern provoke the question, whether there may be a small-scale scatterer somewhere along the path which acts like an obstacle casting a shadow over the AlpArray network. When moving the source to the east, the shadow moves to the west.
Our study includes one pair of events located in almost the same direction but with significantly different distances, see Table 1 and Fig. 6 for the Tajikistan and Indonesia events (Supporting Information Figs S11 and S14). The results for 50 s period are compared in Fig. 9. Beside the strong static shift for the Indonesia event (−12 • ), both events show pronounced narrow stripe-like patterns of arrivalangle deviations. While for the Tajikistan event the stripes have a certain angle with respect to the geometrical ray paths, for the Indonesia event, they are almost parallel to it.
Looking at Figs 8-11, we see that the positive and negative deviations not only have a stripe-like pattern. They also repeat systematically. At each plot, we see at least one negative and one positive stripe-like deviation of the same width. For some events, we obtained even more stripes [Tajikistan, 50 s, Fig. 9 (Supporting  Figure 11. Arrival-angle deviations for eight selected periods for the South Atlantic 2016-08-19 M = 7.4 event. Note that the number of subarrays for each period generally differs, since we only use subarrays with mean residual lower than 2.0 s for the given period. White dotted lines are the great circles. They are the same at all eight maps to allow for comparison among the different periods. Downloaded from https://academic.oup.com/gji/article-abstract/218/1/115/5315763 by Biblio Planets user on 09 September 2019 Information Fig S11); Ecuador events, 90 s, Fig. 10, (Supporting Information Figs S4 and S5) to name the best examples]. If each of the stripes (negative/positive) was caused by its own structural anomaly (low/high velocity anomaly), it would require a very unlikely coincidence of having a structural feature with regularly changing properties (velocity) from positive to negative values repeatedly several times over a scale of hundreds to thousands of kilometres. The repetition of the stripes raises the question whether the deviations are perhaps caused by waveform complexities arising from the interaction with simpler structural features.
A full range of periods is shown for the South Atlantic M = 7.4 2016-08-19 earthquake in Fig. 11 and for all the other events in Supporting Information Figs S1-S19. We observe a systematic shift of stripe-like pattern with period in Fig. 11, from west to east with increasing period. Perhaps also this behaviour might be explained by a very particular type of heterogeneity, with different structural anomalies located in different depths, causing a different effect at different wavelengths. This pattern shift with period is actually observed for all the 20 earthquakes regardless of their distance or azimuth, which again raises the question whether such patterns can be explained by a simpler shape of anomaly, but allowing for wavepropagation complexities given by the finite-frequency interaction with such a simple anomaly.

M O D E L L I N G O F A R R I VA L -A N G L E D E V I AT I O N S
A relatively straight-forward way of modelling the effect of a heterogeneity along the ray paths has been proposed by Nolet & Dahlen (2000), who gave a simple equation predicting the wavefield of surface waves perturbed by a velocity anomaly in 2-D. They defined Q, a function of distance x along the ray, distance R perpendicular to the ray, and angular frequency ω, as a perturbation to a unit plane wave travelling with a phase velocity c, which has the following form: where λ = 2π c/ω is the wavelength, L is the half-width of the anomaly and τ max is the maximum time delay of the initial waveform at x = 0, referring to the point where the ray leaves the anomaly. The anomaly is a simple box-car anomaly placed in a homogeneous space. Its strength and geometry is controlled by τ max and L. The phase delay of the perturbed wave with respect to the original wave is then given as the phase of the complex perturbed wave where τ (x,R,ω) is the resulting phase time delay. In the paper by Nolet & Dahlen (2000), all results are given in terms of four dimensionless variables: x/L representing the distance along the ray, R/L representing the distance along the wavefronts, L/λ controlling the width of the anomaly and ωτ max its strength. Such an approach allows to investigate the wave propagation without considering particular frequencies and velocities. To demonstrate this, in Fig. 12 (top three plots), we reproduce the results from Nolet & Dahlen (2000) for the surface wave fractional phase delay τ /τ max similarly as it is shown in Fig. 6 of their paper. The x-axis represents the distance perpendicular to the ray R scaled by the half-width of the heterogeneity L. The y-axis is the actual time delay τ scaled by the maximum time delay τ max which is set constant for all three plots as τ max /T = 0.188 where T is the period of the wave. In the paper by Nolet & Dahlen (2000), the value of τ max /T = 0.1 was used. We plot twenty lines in each plot representing the fractional time delays τ /τ max for distances starting with x/L = 0 (green line) to x/L = 9.5 (red line) being equally spaced by 0.5x/L. Three plots are given for three different geometries L/λ. If we consider the wavelength λ as constant, the plots differ in the half-width of the heterogeneity L which is the widest in the left plot and the narrowest in the right one. If we consider the heterogeneity to have the same width L, the three plots give us the dependence on wavelength λ, which is the shortest in the left plot and the longest in the right one. The top three plots in Fig. 12 depend only weakly on τ max , since we plot the ratio τ /τ max , and τ scales almost with τ max (as good as the exponential has a slope of 1 close to zero). This can be seen when comparing our plots in Fig. 12 with Fig. 6 of Nolet & Dahlen (2000). Even we use τ max /T = 0.188 instead of their original τ max /T = 0.1, the difference is negligible. For our purposes, however, we are interested in the arrival-angle deviations which can be calculated from the phase time delays τ (x,R,ω). First, we recalculate the time delay to wavefront position w(x,R,ω) = cτ (x,R,ω). The position w(x,R,ω) represents the distance the perturbed wave travelled further or closer with respect to the plane wave. Then, we calculate a derivative of the wavefronts w with respect to lateral distance R. The arrival-angle deviation A is given as arctangent of this derivative and hence it is enough to calculate the lateral derivative of the phase time delay. Results are presented in the bottom three plots of Fig. 12 again for the same three half-width to wavelength ratios L/λ as above. Green lines represent the arrival-angle deviations (again) corresponding to the closest time delay, meaning the deviations right after passing the anomaly. Red lines are the deviations in a distance of x/L = 9.5 corresponding to the red time delay above. It is noteworthy that all the six plots in Fig. 12 depend neither on period nor on wavelength nor on velocity if τ max /T is constant. The arrivalangle deviations are given as a spatial derivative of the time delay. For lower L/λ, this derivative is larger since the time delay is spread over shorter distance. Even the initial time delay (green curve in the top plots) is the same for all L/λ, the space derivative is doubled, if the L/λ is halved, which is clearly to be seen moving from leftmost bottom plot to the right. Even both the time delay τ as well as the arrival-angle deviations A obviously do depend on the velocity, scaled plots like in Fig. 12 do not, because the change of velocity produces also a change of the wavelength λ and as we keep constant L/λ, the L scales with the velocity as well. Increasing the velocity without changing anything else is represented by decreasing L/λ, meaning, we move again from left to right in the plots of Fig. 12. The same applies for the geometry of the anomaly given by L. If we decrease L or increase λ (keeping other ratios constant including the ratio τ max /T), the arrival-angle deviations increase (moving left to right in Fig. 12). The last (third) parameter, which makes the deviations higher is obviously the initial time delay τ max /T itself (not shown in Fig. 12, as all the six plots are drawn for the same τ max /T).

Figure 12.
Top plots: Evolution of fractional phase-time delay with propagation distance, for the three values of half-width-to-wavelength ratio L/λ, similarly as in fig. 6 of Nolet & Dahlen (2000). The horizontal axis is the scaled cross-path distance R/L. 20 curves give the effect at scaled propagation distances ranging from x/L = 0 (green curve) to x/L = 9.5 (red curve). Bottom plots: corresponding arrival-angle deviations in degrees (see the text). All six plots are given for τ max /T = 0.188.
To demonstrate the results in real dimensions, we calculated arrival-angle deviations for three hypothetical earthquakes on a global scale, see Fig. 13. For the top three models, we keep the same three L/λ = 2.0, 1.0 and 0.5 as in Fig. 12. The period (T = 100 s) and the velocity (c = 4.0 km/s) are set so that the wavelength λ = 400 km. The three half-widths of our anomalies are now L = 800, 400 and 200 km. We use τ max /T = 0.188, as in Fig. 12, which, in case of T = 100 s, means τ max = 18.8 s. Such a time delay may be given, for example, by a heterogeneity of 1000 km length and 7% strength. We consider low-velocity anomalies in Fig. 13 (red rectangles). The ratio τ max /T = 0.188 was chosen so that the highest arrival-angle deviation reached for L/λ = 0.5 equals ±20 • , see the right bottom plot in Fig. 12, green line. The magenta rays in Fig. 13 correspond to magenta vertical lines in Fig. 12 representing R/L = ±1 to ±5. Green fractional time delays from the top plots of Fig. 12 are now recalculated to distorted wavefronts in Fig. 13 using the velocity c = 4.0 km/s. This applies as well for all the other time delays (wavefronts) with a step of 0.5x/L ending with a red wavefront at the distance of x/L = 9.5. Colour scale of arrivalangle deviations is also the same in both Figs 12 and 13 reaching ±20 • (what equals to all the other colour scales throughout this paper).
Bottom three models in Fig. 13 show the results for the same geometrical settings with period decreased to 50 s. The three anomalies are kept exactly the same, meaning, their lengths and strengths as well as the velocity of the homogeneous medium (c = 4.0 km/s) around the anomalies are kept the same. The absolute value of τ max is also kept at 18.8 s, which means that the ratio τ max /T = 0.376 is doubled for the bottom three models with respect to the top ones, where it was τ max /T = 0.188. The wavelength is halved to λ = 200 km and hence the ratios L/λ are now doubled to 4.0, 2.0 and 1.0. Even the doubled L/λ means the arrival-angle deviations are halved, the double increase of τ max /T with respect to top three models now produces more than doubled arrival-angle deviations since the dependence of arrival-angle deviations on τ max is exponential. As a result, the arrival-angle deviations for T = 50 s at the bottom three models are, at their maximum, higher than those for T = 100 s in the top models. Also note that as opposed to the top three models for T = 100 s, the highest arrival-angle deviations are now found not right after passing the anomaly but, for example, at the distance of x/L = 7.5 which corresponds to 3000 km from the L = 400 km anomaly, see the middle bottom model.
The arrival-angle deviations are calculated to a distance of 10 000 km from the anomaly for all six cases and for a width that slightly exceeds the range of R/L = ±5. We see the patterns observed and discussed above. Smaller anomaly produces narrower lobes of deviations with higher amplitude. The lobes are also tilted with respect to the geometrical ray paths. A broader anomaly produces wider lobes with lower amplitude. They persist for longer  Fig. 12 and are doubled in the bottom models. To allow comparison, a rectangular representation (oblique Mercator projection) has been chosen for the Earth's surface. Magenta lines correspond to the magenta lines in Fig. 12 showing the scaled cross-path distance R/L. distance, being more aligned with the geometrical ray paths than those from a narrower anomaly. We see that such an anomalies can produce deviation patterns consisting of lobes of changing sign and repeating in space. Considering the size of our observation area covered by the AlpArray reaching approx. 10 • north-south and 20 • east-west (one rectangle of the meridian-parallel grid in Fig. 13), we see that we can indeed expect to observe two to four lobes of deviations depending on the anomaly size and its distance to the station network.

Main results
We have observed intriguing spatial patterns of arrival-angle deviations for essentially every earthquake. These patterns generally have a stripe-like form elongated in the direction of geometrical propagation or with a certain angle to it. The patterns change much less with distance from the event than in the direction transverse to the geometrical ray paths. The same was observed also by Foster et al. (2014a), and called 'banded appearance' therein. On the other hand, we do not see any similarity in the patterns for events from different directions, which indicates that the patterns are neither due to the subsurface under the AlpArray network nor the stations themselves. If the patterns were caused by structural anomalies within the network, it would be also very unlikely to see them always aligned with the azimuth of the event. Clearly, the observed deviations are caused by structural features outside the region. The wavefronts propagating across AlpArray had already been distorted before they reached the region.

Confirmation of the stripe features
The fact that earthquakes at the same location (pairs A-D) produce similar spatial patterns of arrival-angle deviations can be regarded as a confirmation of the robustness of the measurement. That the patterns are real propagation features is also confirmed by the comparison of the two different techniques in Fig. 8, where it was shown that the array measurements produce rather similar results to the direct determination of wavefronts from single-station phase measurements. The latter can serve only for comparison though, since the array measurements of phase (velocity and) arrival direction are based on the entire waveform rather than a single phase measurement in time. As an output of the beamforming, we obtain directly the values of the arrival angles and hence the wavefront distortion. From the single-station measurement, we can only interpolate the contours to the map what involves smoothing and hence the results do not yield clear details as the array method does. The array analysis also allows more rigorous data quality procedures (see the Appendix, part 'Stage 5') which is missing for the singlestation phase picking. Several maps, for example, in Fig. 8, show small point-like features though, which probably represent errors in individual measurements (remember that the maps have not undergone any smoothing). Comparison of array and single-station wavefront measurement was done by Foster et al. (2014a) in a similar way also concluding that array approach allows to characterize the wavefronts quantitatively.
Magnitudes of the teleseismic earthquakes we used span from 6.6 to 8.1 (not accounting for the regional M = 5.5 Greece event), see Table 1. Finite size of the rupture, in general, affects surface wave propagation. The question whether the finite rupture influences also the arrival-angle deviation observation can be assessed by looking again at the pairs of collocated events. Excluding the two events in South Atlantic Ocean (pair D, Table 1) which have similar magnitudes (7.4 and 7.2), we see that the other three pairs (A-C) differ in magnitude a lot (6.8 vs 7.7 for the Aleutian and Komandorskiy events, pair A; 6.9 vs 7.8 for the Ecuador events, pair B and 6.9 vs 8.1 for Guatemala and Mexico, pair C). Although we can note differences at the longest periods where the signal is worse for the smaller events, the observed patterns are so similar that the effects of different source sizes can be considered as negligible. Compare, for example, 50 s deviation maps of Guatemala and Mexico events in Supporting Information Figs S16 and S18. They look almost the same both in shape and size even the rupture length of the Mexico event was about 250 km while for the Guatemala event it was first tens of kilometres only.

Average properties
We will investigate statistical properties of the arrival-angle deviations from the array beamforming. Fig. 14 shows histograms of arrival-angle deviations (with 1 • bins), plotted for selected periods (right upper corner of each plot) and compiled using the results for all twenty earthquakes. The number of subarrays is given for each plot as well. Since our histogram distributions are not perfectly normal, we use the 68.2 nd percentile instead of the standard deviation to characterize the width of the distributions. This value is shown by light green lines as negative and positive range with respect to the mean depicted by dark green line and given by number in upper left corner of each plot. If the distribution was normal, that percentile would correspond to standard deviation. However, standard deviation is sensitive to outliers while the percentile is not. Using the mean values and 68.2 nd percentile as a replacement of the standard deviation, we calculated the normal distributions plotted by blue lines in Fig. 14. This shows to which extent our measurements follow the normal distribution. Around 60 and 70 s they differ, however, for other periods the histograms copy the normal distribution fairly well.
The resulting distribution of arrival-angle deviations is summarized in Fig. 15 depending on period. The width of histograms from Fig. 14 becomes smaller for longer waves (magenta line in Fig. 15, left vertical axis), however, the change between 40 and 150 s is only 2 • . The distributions are significantly wider only for waves shorter than 40 s sampling mostly the crust. However, in that shorter-period range, the number of measurements is also much lower, as shown by the thick grey line (the total number of subarrays for the given period, right vertical axis in Fig. 15). The thin grey line shows the number of subarrays with deviations falling within ±30 • range, which is the interval shown in Fig. 14. The difference between the two latter curves represent the insignificant number of outliers falling over ±30 • . The span of arrival-angle deviations therefore does not appear to be significantly decreasing with period, although we see a slightly lowering trend. This is somewhat in contrast with a recent study by Chen et al. (2018), who reported significantly decreasing deviations for longer periods. Their study was limited to a period range 30-60 s. However, it even emphasizes the difference to our findings, since we see constant widths of arrival-angle deviation distributions from 40 to 75 s. The deviations slowly decrease only for waves longer than 75 s and increases for waves shorter than 40 s. A period range closer to ours (25-100 s) was used by Foster et al. (2014a). They also reported the deviations being smaller for longer waves, however, they also say that for all the periods, the deviations span the same range. They see a greater number of deviations of higher size for shorter waves. Neither of the latter two studies, however, performed a statistical assessment of the deviation distribution and the statements therein are based on visual inspection of the respective deviations maps only. Also note that our statistics is based on measurements with uneven distribution of subarrays. If there is, for example, a positive deviation in a region with higher subarray density, we get more counts in our histograms than if there was the same deviation in a region with sparser subarray locations.
The green line in Fig. 15 shows the mean absolute deviation calculated the same way as in Pedersen et al. (2015), meaning, it is a mean of absolute values of all the deviations. It has the same shape as the magenta line showing the width of distributions. Comparing our results with Fig. 3 in the latter study, we see that our mean absolute deviations are higher by about 1 • than those presented by Pedersen et al. (2015) for waves shorter than 40 s, matching, however, perfectly at 40 s and then at 90 and 100 s period. In the range between 40 and 90 s, our mean absolute deviations are higher by about 2 • , especially for 60 and 70 s period. This is not surprising, since the array size used by Pedersen et al. (2015) was 275 × 460 km while our subarrays have an aperture of only 160 km. If there were stronger deviations forming sharp stripes (Aleutain 40-65 s, Supporting Information Fig. S15; Komandorskiy, 40-50 s, Supporting Information Fig. S17; Ecuador M = 7.8 and M = 6.9, 50 s, Supporting Information Figs S4 and S5 to name some), they would be smoothed down by an array of the LAPNET size. Our array aperture (160 km) matches roughly one third of the wavelength of the longest period (125s, ∼500 km) and around the whole wavelength of the shorter periods (40 s, ∼150 km). Even using such small subarrays, we also cannot exclude that some narrow stripes of deviations were smoothed down. However, looking at Fig. 13 and noting that the distance between the positive and negative stripes is approximately a wavelength right after the anomaly, where the stripes are the narrowest, we see that an array size of one wavelength is sufficient to resolve even those narrow stripes well.
The red line in Fig. 15 (left vertical axis) represents the mean value of deviation distributions (also given by numbers in Fig. 14). We see a smoothly changing curve with negative values most of the time and slowly reaching zero for the longest period. The negative values are apparently associated with a few events at large distance, especially the Indonesia (Supporting Information Fig. S14) and the Mariana (Supporting Information Fig. S7) events, which show a static negative shift of arrival-angle deviations. That static shift decreases with period. However, we may still say that according to our findings, the mean deviation is not zero in the Al-pArray region as expected by Kästle et al. (2016); Soomro et al. (2016) and Meier et al. (2004), who assumed that the deviations Figure 15. Mean arrival-angle deviation depending on period (red line, see also dark green lines in Fig. 14), width of the distribution corresponding to the 68.2 nd percentile (magenta), total number of measurements (bold grey) and number of measurement in the ±30 • range (thin grey line). Green line shows mean deviation calculated from the absolute values of all deviations. For this figure, also the events from the Supporting Information have been used. cancel out after averaging the ray paths over many events. On the other hand, the number of events in our study is limited, however, still representing broad range of azimuths. This raises the question whether the observed non-zero mean could possibly explain the discrepancy between the earthquake-based two station phase-velocity measurement and ambient noise results as pointed out by Kästle et al. (2016).

Scale and shape of the deviations
A clear observation from the colour maps is that the patterns are elongated along or close to the direction of geometrical propagation, with along-great-circle lengths that are mostly much larger (spanning the whole region of observation ∼1200 km) than the size of the subarray (160 km). Perpendicular to the great-circles, there are half-widths (distance between ridges and valleys) that can be as small as 200 km, but can go up to 600 km and more. This agrees with results by Foster et al. (2014a) who reported the width of the stripes (bands) being from 2 • to 5 • (∼220-550 km). The strikes are mostly parallel to the great-circles, but occasionally deviate up to 20 • , and very rarely up to 30 • (for a few short-period measurements). This is again similar to Foster et al. (2014a) who observed the stripes both parallel or deviated from the great circle paths. In their study, in addition, some of the stripes are also curved while others are linear. We do not see any curved stripes probably due to the AlpArray being much smaller than USArray. We also concluded that the patterns originate from a structure exclusively outside of the AlpArray region, while Foster et al. (2014a) hypothesized seeing the effects of structure also within the USArray since some of the bands of deviations evolved across the North America. Nevertheless, they attributed much of the deviation complexity to originate outside of the USArray region as well.

Comparison with the modelling and potential inferences
We have seen in Section 6 that the observed deviation patterns can in principle be predicted by finite-frequency effects from simple spatially-isolated anomalies. We will not attempt to match specific patterns here (we defer this to a future paper), but we address the question whether the general features can be explained. In our simple model, we have three parameters, which control the azimuthal deviations. These are the distance of the anomaly from the point of observation, the width of the anomaly and its strength given by the initial time delay τ max , which is a product of the anomaly length and velocity perturbation. From an observational point-of-view, the effects fall into four categories: (i) The amplitude of arrival-angle deviations: We observe arrival-angle deviations commonly up to ±20 • and we model them in the same range. The stronger the anomaly, the higher are the deviations. This is not shown in the figures, but it is implied directly from the eqs (2) and (4). We also see that the narrower the anomaly is (in the range tested in Fig. 13), the more pronounced are the deviations. The deviations decrease with distance from the anomaly for Downloaded from https://academic.oup.com/gji/article-abstract/218/1/115/5315763 by Biblio Planets user on 09 September 2019 the two main (most central) stripes, however, they increase for the other side stripes. Stronger anomaly of the same width or narrower anomaly of constant strength can be considered as having higher lateral gradient (first derivative) of the phase velocity. According to ray-theoretical predictions by Woodhouse & Wong (1986), this would lead to higher arrival-angle deviations as well.
(ii) Width of the stripes: From Fig. 13 we see that wider anomalies cause broader stripes. The stripes also become broader when moving further from the anomaly even though this effect is not that large as it might be seen in Fig. 13 since the map projection distorts the distances. Stronger anomalies, on the other hand, make the stripes narrower.
(iii) Alignment of the stripes: Rather than 'V' shape (constant misalignment), the stripes are shaped like 'U' (misalignment changes with distance), meaning that the stripes become more aligned with the event-station great-circle paths for larger distances from the anomaly. A wider anomaly makes the stripes more aligned with great circles at shorter distances from the anomaly; a stronger anomaly does it as well. If the anomaly is placed closer to the source, it would be necessary to take the effects of circular wavefronts into account. This would also lead to changes of the stripe alignment.
(iv) Shift with period: In Fig. 13 (bottom plots), we use the same three anomalies as in the top plots, but we modelled the effects for the period of 50 s (rather than 100 s shown in the top plots). We see that for shorter period, the patterns changed in several ways [referring to the cases (i)-(iii) above]: (i) The arrival-angle deviations are higher. (ii) The stripes are narrower. (iii) The stripes are more aligned with the great-circle paths. Merging the width and alignment of the stripes into one observable, we may say that for longer periods, the 'U' shaped pattern gets more open (wider). This is not surprising considering that the pattern is a result of interference and that longer waves are expected to produce interference patterns with broader variations. However, the other notable feature is that lowering the period also increases the deviations, meaning, they persist to longer distances. The best example is the 800 × 1000 km anomaly in the middle of both top and bottom models in Fig. 13. While for the 100 s wave the dark red (dark blue) colour of the first positive (negative) lobe on the left (right) disappears around 7000 km, for the 50 s wave it persists to the edge of our modelling at 10 000 km distance from the anomaly.
Let us now discuss the resemblance of the models to our observation. The top middle example in Fig. 13 for a 100 s wave suggests that we can very well expect stripes with half-widths of 400-600 km at distances between 1000 and 5000 km from such an anomaly (considering the area of observation is represented by one rectangle of the meridian-parallel grid). Stronger anomaly would make the stripes even narrower. For the top right-most case in Fig. 13, the half-width of the stripes is roughly halved. Concerning the alignment of the stripes, we see that deviations of the direction of stripes of up to 20 • or 30 • from the great-circles (as observed) can easily be modelled. Compare also Fig. 13, top models for 100 s with results in Fig. 10 for 90 s wave and Fig. 13, bottom models for 50 s wave with Fig. 9 for the same 50 s wave. In other words, at a given period, the width of the stripes is controlled by the distance to the anomaly and by its width. Keeping the geometry fixed and changing the period, the modelling shows that longer periods produce broader stripes. It also means that if the width of the stripes is (almost) not changing with period for a given event, it needs to be the geometry of the anomaly which is varying with depth. The anomaly needs to be weaker or broader (or both) for shallower depths to produce the stripes of the same width as for longer waves (greater depths). The modelling shows that the same geometry produces significantly lower deviations for longer waves, which is not necessarily the case in our observation. It again means that if we see the deviations being almost the same in magnitude across the whole period range, the anomaly again needs to be weaker or broader at shallower depths to produce comparable deviations for both short and long waves.
These parameters, and especially the fact that all of them play simultaneously should be rather constraining on the type of anomaly that causes these effects. It will impose observational boundary conditions on the nature of the anomaly, for example, its location and its strength: (1) The arrival-angle deviation values give us a guess about the strength of the anomaly (stronger anomaly -higher deviations) as well as its width (wider anomaly -lower deviations) with a trade-off between the two. (2) Width of the stripes tradesoff between the strength (stronger anomaly -narrower stripes) and width (wider anomaly -broader stripes) of the anomaly as well. (3) The alignment to the great circle gives us a guess about the two parameters of strength (stronger anomaly -more aligned) and width (wider anomaly -more aligned) as well. The direction in which the stripe is tilting with respect to the great-circle determines to which of the two sides the anomaly lies with respect to the great-circle path between the earthquake and our station (array). Obviously, everything trades-off with the distance.
Our demonstration in Fig. 13 shows that already isolated anomalies can cause patterns as the ones observed. More complex heterogeneity will do this even more. Staying with the case of isolated, single and simple anomalies, the question is, how probable it is that such an anomaly is seen by waves propagating from all the events we processed. The answer is suggested by Fig. 13. Not only are the deviations pronounced out to great distances from the anomaly, but they also spread to considerable distances laterally. In terms of distances along the propagation paths, it is enough to have such an anomaly anywhere between 3000 and 10 000 km from the area of observation to get similar deviations (between ±10 • and ±20 • ) across the whole period range. The main differences in deviation magnitude depending on period occur only in the region close to the anomaly (first thousands of kilometres). If, in addition, the anomalies are weaker or broader (or both) for shallower depths, it equalizes the deviation size across the period range as discussed above. Laterally, it is enough to encounter such an anomaly anywhere within ±2000 km off the great-circle ray path to see the stripes it produces. If the anomaly is more distant, it can be place more off the path. It means that it is enough to have only a few such anomalies to produce such patterns at wavefields propagating from almost anywhere. Pedersen et al. (2015) already suggested that the term 'multipathing' may be misleading since it is mostly used for late arrivals travelling significantly different paths than the 'direct' wave. However, our observation supported by the simple model suggests that the deviations are rather a result of interference of non-plane waves scattered after passing the anomaly even though they are still propagating as a 'direct' wave. The same applies for the term 'great-circle deviation', also discussed in Pedersen et al. (2015). Again, we suggest not to assign the observed stripe-like arrival-angle deviations to 'great-circle deviations' but rather to wavefield complexity, which might be called scattering or diffraction off laterally heterogeneous structure.

Effects along propagation paths
Among our twenty events, there are examples demonstrating the difference between the great-circle deviation and diffraction of 'direct' waves. Looking at Fig. 9, right map in the middle row for the Indonesia earthquake, we see a general shift of the whole pattern by −12 • . This systematic deviation could be attributed to what is commonly considered as 'off-great-circle propagation' or 'great-circle deviation'. Pedersen et al. (2015) shows examples of the deviation from the same region as our Indonesia 2017-05-29 earthquake as well, even right in this region, a small change of the event position resulted in a deviation of opposite sign. Isse et al. (2006) discussed the ray bending in the Philippine Sea showing a systematic deviation of the true rays from the great-circles along the whole path. Cotte et al. (2000) concluded that waves from Japan propagate north of the great circles avoiding the Tibetan Plateau. Levshin et al. (1994) gives an overview of how the waves from East Asia propagate through the fast Siberian platform rather than through the slow Tibetan plateau as well. We can attribute that systematic shift to avoiding the Himalaya/Tibetan Plateau region as well (even the particular cause is not important in the moment). In addition to this systematic shift, however, we see a stripe-like pattern modulated over the −12 • deviation similarly as in other cases, namely strikingly similarly as in the case of the Tajikistan event (again Fig. 9, left of the Indonesia deviation map). These smaller-scale stripes can be attributed to interference after a diffraction caused by isolated anomaly somewhere on the way of the already greatcircle-deviated waves from Indonesia. Such a differentiation has already been proposed by Alsina et al. (1993). They observed two kinds of deviations at the shape of the wavefronts. First was the 'difference in the direction of arrival'. It was believed to be caused by the off-great-circle propagation (multipathing) exclusively due to a structural anomaly outside of the region of observation. This would correspond to the −12 • of our Indonesia example. Second, they described the 'difference in curvature' of the wavefronts as an evidence of lateral heterogeneity both below the region of observation (Iberia) and along the entire path. We admit that the discrimination between the systematic deviation and stripe-like pattern is somewhat arbitrary in our case, as, of course, even the systematic shift can be attributed as being caused by isolated anomaly of greater size far away from the array of stations so that the whole AlpArray network falls into only one stripe of negative arrival-angle deviation. Cotte et al. (2000), Pedersen et al. (2015) and Foster et al. (2014a) noted that sometimes only a small change in the location (direction) of the event (or similarly the same direction but different distance from the source) gives the opposite sign of the deviation. We can explain this by our model, referring to the two events from the South Atlantic discussed above. As we observe the interference pattern from a single anomaly somewhere on the way from both events, the pattern just shifts over the network producing an opposite sign of deviations at single subarray. The same applies also for the frequency dependence. Maupin (2011) observed again changing sign of arrival-angle deviations for different periods in southern Norway. As the interference pattern widens with increasing period (compare top and bottom models in Fig. 13), we can easily observe repeatedly changing sing of deviations as the 'U' shaped stripes move across the station (subarray). Maupin (2011) also noted that the waveforms she processed already look like being a result of interference showing the typical beating behaviour. She attributed this to multipathing discussing that several wavegroups interfere one with the other. Although for some cases she found several beams corresponding really to several distinct arrivals, for other examples she concluded that in the scale of her network, the different arrivals interfere in a way that they behave as a single wave train. Such a case corresponds to our findings. In our analysis, we have selected a single wavegroup only (of 4 periods length) and still, we observe the interference within even such a short wave train. We prefer to not call this 'multipathing' but rather interference due to scattering of waves. Foster et al. (2014a) suggested that wider stripes (bands) observed for longer periods reflected both increasing wavelengths as well as expression of smoother mantle structure. Although we agree with both parts of the statement, we need to point out that in our modelling the second is even not necessary. An isolated anomaly of the same properties over arbitrary depth (column-like shape) is enough to provide wider stripes for longer waves due to the interference. Hence wider stripes of longer waves do not necessarily mean the mantle is less heterogeneous. Wielandt (1993) has shown that it is not necessarily the structure what leads to complicated observations. The deviations we measure can be caused by the complex interference between the propagating waves themselves. We believe that understanding the wave propagation itself, namely the cause of arrival-angle deviations, is important before we make attempts to obtain the local structure.

C O N C L U S I O N S
We have introduced a technique for mapping the phase-velocity vectors by a dense broad-band seismic network such as AlpArray, using beamforming of pre-selected fundamental modes of Rayleigh waves. This study has focused on arrival-angle deviations only, to study path anomalies, which we believe will form an important ingredient for also better understanding phase velocities. In addition to benchmark the latter technique, we also investigated phase wavefronts using single-station tracking of zero-crossings of filtered quasi-harmonic signals. Twenty earthquakes were selected from the first two years of the AlpArray deployment to provide good azimuthal coverage. The (more than) 600 000 arrival-angle deviation measurements form characteristic patterns of elongated stripes that are roughly aligned with the propagation direction, for every event.
The good spatial distribution of AlpArray allows understanding their cause: using a simple isolated boxcar anomaly model we show that stripes like the ones observed can be produced by scattering, as an alternating interference patterns. That explains, beside the regular spatial stripe-like patterns also the slow and systematic variation of these patterns with period. In addition, we are able to explain abrupt changes of arrival-angle deviations observed by earlier studies from earthquakes of only slightly different location. The fact that single simple anomaly can explain the deviation patterns we observe does not exclude that other complicated geometries would produce similar deviations. However, we conclude that the complexity we see is predominantly caused by intricate wave-propagation phenomena rather than by the structure. We have shown that properties of the observed deviations (size, width, alignment, period dependence) allow guessing the properties of the anomaly. Our findings suggest that arrival-angle deviations constitute important information Downloaded from https://academic.oup.com/gji/article-abstract/218/1/115/5315763 by Biblio Planets user on 09 September 2019 on subsurface structure, beyond what is normally used in inversions. Modelling of particular anomalies in the upper mantle causing the observed arrival-angle deviations will be a topic of future study.

A C K N O W L E D G E M E N T S
We are grateful to Gidera Gröschl for downloading the data, providing the data quality control and pre-processing. Florian Fuchs contributed to the data quality report mentioned in the Appendix. We acknowledge financial support by the Austrian Science Fund (FWF) through project P 26391-AlpArray Austria and P 30707-AlpArray Austria 2. The Python Toolbox ObsPy by Beyreuther et al. (2010) was used for data pre-processing. Maps were plotted using Generic Mapping Tools by Wessel et al. (2013). We wish to acknowledge the editor Gabi Laske and three anonymous reviewers, whose comments gave important additional insight that significantly improved the manuscript.

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.
Figures S1-S19. Online-only Supporting Information consists of Figs S1-S19. They follow the same scheme as Fig. 11. Arrival-angle deviations for eight selected periods for each event are shown. Note that the number of subarrays for each period generally differs, since we only use subarrays with mean residual lower than 2.0 s for the given period. Arrows show the direction of wave propagation. White dotted lines are the great circles spaced equally with a distance corresponding approximately to a wavelength of 50 s period. The position of the great-circle lines is kept for all the periods to allow for comparison of the stripe position among the different periods. For the Greece event (close to the AlpArray region, Fig. S9), the great-circle lines are set to be 10 • apart. For the Indonesia event (Fig. S14), we added 12 • to all the arrival-angle deviations, see text and the caption to Fig. 9.
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 : DATA Q UA L I T Y
The data quality was investigated in the following five stages, from a rough inspection to very fine timing errors investigation. Our very strict quality policy led to a considerable amount of data being discarded. When any doubt about the record arose, we did not use it.
Stage 1: gap check. For each event, we cut out a record starting around the origin time and spanning from 1.5-2.5 hours depending on the distance of the event. We browsed all the three components checking for gaps, overlaps, missing traces and short traces. If there were more than 20 gaps in the record during the selected time window or if any of the gaps was longer than 3 s, we removed the station. This rule was set to allow for the maximum gap of 1/10 of the shortest period we had planned to investigate (30 s). For the event in Greece, we had planned to study also shorter periods and hence the maximum allowed gap was set to 1 s. If the gaps were shorter, the traces of the record were merged and gaps were filled with interpolated values. Overlapping traces were also merged. If a component was missing, the station was also discarded. In our study, we use only the vertical component, however, we prepared the data for future use of the horizontals as well. This stage led to the biggest loss of data, usually tens of stations were removed for each event mainly because of gaps.
Stage 2: visual inspection. If the record passed the previous stage, it was rotated for geometrical Z-R-T components, filtered between 1-200 s and resampled to 10 sps to have uniform sampling rate for all data. All records for the given event were plotted, sorted by distance to the source and visually inspected. We aimed at zero traces, electronic noise, too noisy records where the earthquake signal was not visible and records heavily shifted in time. Looking at main body wave arrivals, we were able to discover a time shift greater than 2 min (more than 100 s). During this step, just few stations were lost for each event.
Stage 3: group velocity dispersion. The group velocity dispersion itself was not needed for our study, however, we used it to taper out the fundamental mode wavegroups. The group velocity dispersion was used to compare all stations for a given event with an average group velocity dispersion curve for the region (calculated as an average of all those measured for the given event). In the range where the fundamental mode dominated the amplitudes (usually 50-150 s), we marked the stations for which the dispersion curve differed by more than 5% from the average. These stations were manually investigated to decide if a change of filtering parameters helped to improve the fundamental mode identification. If not, the station was removed. During this step, we again lost up to tens of stations for each event. The main reason why the fundamental mode could not be identified were noisy records and timing problems. For some stations, the group velocity dispersion curve was well-visible, however, significantly off in time from the other stations. By this method, we were able to remove stations with a time shift less than 100 s (passing the visual inspection). However, the group velocity dispersion is still affected by significant measurement error and so time shift less than 10 s could not be revealed.
Stage 4: wavefront maps. During the previous analysis, the zero crossings of the records (phase times) were measured. We plotted them in a map in terms of contours for every 10 s of propagation. Visually inspecting the map of the region allowed to remove stations for which the phase time did not correspond to the neighbouring stations. This way of discrimination allowed for detecting timing errors smaller than 10 s but still more than 3 s. Smaller problems were not visible in the map. Around 10-15 stations were removed during this step for each event.
Stage 5: array analysis. Records passing all the previous four steps were used for the array analysis. Each station was checked for the geometry of the surrounding stations and if it had more than five neighbours, beamforming was provided for that subarray. To evaluate the results we used the residuals mentioned in Section 2.4. First, we looked at total residuals averaged over the whole frequency range and all neighbouring stations for a given subarray. All subarrays having total residuals > 2.5 s were manually investigated. There could be two reasons why the plane-wave fit was not successful. The problem could be caused by the central station of the subarray itself or one (or more) of the neighbouring stations. So, for all these suspicious subarrays we inspected the residuals for each neighbouring station (averaged over the frequency range). If they all showed high values, the problem might be the central station. If only one (or more, but not all) neighbouring stations had high residual, we investigated what the problem was using the finest residuals depending on frequency. When these were high for only the very shortest or very longest waves, we might still decide to use the station removing the measurements for extreme periods or to narrow the period range. If it was consistently high over the whole Figure A1. The 890 stations for which the data were downloaded for our study (see also Fig. 1). Green triangles show stations used as subarray centres with colour shades marking the times each of them is used (see legend), among the (maximum) 20 events. Dark blue triangles mark the stations used at least once and always exclusively as neighbouring stations. Light blue triangles show the stations used at least once and always exclusively for the phase time contour maps but not involved in the array analysis. In total, we used 853 stations. Red triangles show the stations downloaded, but never retained for measurement. Yellow line shows the AlpArray region. range of periods, meaning the phase was simply off in time, we removed the station.
It is obvious that one station with a shifted phase destroys the results for all the subarrays it is involved in. So, by removing it, we lowered the total residuals for several subarrays around it. If it was not clear from investigating the individual station residuals, which of the neighbouring stations caused the problem, we used jackknifing. By removing the stations one by one and recalculating the plane wave fit for the given subarray, we found the problematic neighbour by watching for the significant drop in total residual.
Usually, about 10% of subarrays were listed as having high total residual, meaning roughly 50 out of 499 on average for each event. By investigating these 50 subarrays, we usually needed to remove about 1/3 of the number of stations (15)(16)(17)(18)(19)(20) to lower the total residuals for all of the 50 affected. After all these problematic records were removed, we run the array analysis again for all the remaining stations for given event. Table 1 summarizes the numbers of stations used for each event. We see that the ratio between the downloaded number of stations (749 on average) and the number of stations used for the analysis (606 on average) differs by roughly 19%. This is the price we pay for having the records with indubitably good quality. For each event, the used stations fall into three categories: (1) the stations which could be used as a centres of the subarrays ('subarr' in Table 1), (2) the stations, which contributed to the subarray calculation ('nghbrs') but they could not serve as a central stations having not enough neighbours and (3) stations which were not used in the array calculations at all ('sngls'). These were, however, still used for the phase time contours plotting. It is clear that these latter stations were not included in the last Stage 5 of the data quality assessment. The last category (4) of the stations are those which did not pass the first four steps of quality check and were not used at all for any measurement ('not' in Table 1). Fig. A1 shows the distribution of the stations by their use by colour assembled for all the 20 events. The times each of the stations was used as a centre of subarray is shown by shades of green, see the legend. Blue triangles mark the stations, which were used at least once as a neighbouring station and never as a centre of a subarray. Light blue triangles show the stations never used for array measurement, however, still passing the first four steps of the quality check and being used for the phase time contours maps. Red triangles show the stations being downloaded at least once but never used for any measurement. These did not meet our quality criteria at all. Blue (neighbours), light blue (singles) and red (not used) do not show the distribution by number of events involved. The figure shows the use of stations overlaid for all the 20 events by the hierarchy subarrays-neighbours-singles-not used. Meaning that event to event, the function of the same station could change. Some of the stations which had had the function of neighbouring stations for events with smaller number of subarrays might change into a centre of a subarray later for another event. The same applies for shifting the function from neighbours over the former singles and so on. Fig. A1 shows two aspects of the records used for subarray measurement: availability and quality. The reason, why the western part of the AlpArray region contains more dark green triangles is that these stations were deployed later than those in the eastern part and so they could not be used for all the events. The same applies also for temporary AlpArray stations, which were moved during the experiment. These are marked by different triangles considered as different stations. Obviously, the old position could not be used for the recent events and the new position could not count for the older events resulting in darker green of the triangle. However, every single station could be used for less events also because poorer data quality. The reason we show the figure using colour shades for times of use for the subarray stations is that we have presented the results for varying regions for each earthquake. For the older events, we were able to map only the lighter (green) part of the Alpine region. For more recent events, we were able to move with the array measurement to the western edge of the AlpArray project (dark green). Table 1 summarizes also the number of all above mentioned use of stations for individual events. Although we had 741 unique subarrays (stations serving as a subarray centre for at least one event), the maximum number of subarrays for single event is 588 (Komandorskiy Island). On average, the number of subarrays is 499 per event. The number of subarrays roughly follows the same distribution over time as described above for the number of downloaded stations. We also see that although we had only 37 stations which were never used for any calculation, on average, the number of stations not used per event was 143.
As a byproduct of our study, a comprehensive report of data quality has been compiled by the authors and colleagues and distributed to the AlpArray community. It is accessible at the project web page www.alparray.ethz.ch. The report includes the quality inspection described by stage 1 and stage 2 and spans over 35 earthquakes from 2016 and 2017. The 20 events presented in our study is a subselection from these 35 events.