Localization of seismic events produced by avalanches using multiple signal classiﬁcation

S U M M A R Y Snow avalanches generate seismic signals as many other mass movements. Therefore, detection of avalanches by seismic monitoring systems is highly relevant to assess the avalanche danger level. Apart from estimating the overall activity, information on the location of individual events is also important. We deployed a seismic monitoring system consisting of seven vertical geophones at a site above Davos, Switzerland. During a major snowfall on 2017 March 9 and 10, we visually identified 19 events in the seismic data that we considered as avalanches. For all these events, we calculated the backazimuth angles as well as the apparent velocity of the incoming wavefields using two different array processing techniques: a commonly used beamforming approach and multiple signal classification (MUSIC). One of the avalanches released directly above the array and the path was well documented in a field survey. We analysed this event in more detail and compared both array processing techniques. Using the MUSIC method we were able to reproduce the avalanche path with reasonable accuracy, whereas the beam-forming method only provided a rough direction. The results obtained using MUSIC suggest that with our seismic array mainly infrasonic waves were recorded during the event. From the remaining 18 events, 15 had a clear path and were classified as avalanches. The other three events were dismissed as falsely classified avalanches. By comparing the results for these 15 avalanches events with our field survey, we were able to map 11 avalanches within a distance of 3 km of the array.


I N T RO D U C T I O N
Seismic monitoring methods are well suited for the detection of snow avalanches. Since the 1970s, various studies have shown that snow avalanches generate seismic events strong enough to be recorded using seismometers and simple geophones (Harrison 1976;St. Lawrence & Williams 1976). The focus has been mostly on understanding signal characteristics and estimating avalanche characteristics from the seismic signal (Sabot et al. 1998;Suriñach et al. 2000;Suriñach et al. 2001;Vilajosana et al. 2007), or on monitoring avalanche activity (Leprettre et al. 1996(Leprettre et al. , 1998Bessason et al. 2007;Rubin et al. 2012;Heck et al. 2018a). While it is clear that seismic methods can be used to monitor avalanche activity, information about the type, size and location of the avalanche would be an added value.
Snow avalanches also generate infrasound signals (Bedard 1989) and those signals can be used to detect avalanches (e.g. Marchetti et al. 2015;Thüring et al. 2015). During windy conditions, however, avalanche detection based on one single infrasound sensor is prone to errors (Scott et al. 2007). Hence, for the infrasonic monitoring of avalanches arrays of several sensors are used. Apart from the increasing robustness of detection, array analysis methods can be used to derive additional information about the avalanche event, such as the velocity of the avalanche front (Havens et al. 2014) as well as the direction of the source (Marchetti et al. 2015). Kogelnig et al. (2011) combined seismic and infrasound data to identify different avalanche regimes and by combining both data they identified additional characteristics.
Only a very few studies have focused on estimating the size and type of snow avalanches from seismic signals. To do this, the focus was mostly on signal length and frequency content (van Herwijnen & Schweizer 2011a;van Herwijnen et al. 2013;Pérez-Guillén et al. 2016). While these approaches can give indirect information on the size and type of the avalanches, the location of the path still remains unknown. Furthermore, these studies mainly used data from single stations or one single sensor of a seismic array. As for infrasound monitoring, seismic array analysis can also be used to obtain more robust results of avalanches.
Source localization through array processing is commonly used in seismology, for instance, to determine the epicentre of earthquakes or the localization of volcanic tremors (e.g. Almendros et al. 2014). One of the most commonly used techniques is beam-forming, which identifies coherent phases in the seismic signal by applying crosscorrelation (Rost & Thomas 2002). The beam-forming method is perfectly suited for sharp phase arrivals, e.g. the P-wave onset of an earthquake or to identify the arrival of surface waves and the apparent velocity as well as the backazimuth of the incident seismic wavefield can be determined.
An alternative approach is the multiple signal classification (MU-SIC) (Schmidt 1986). In contrast to beam-forming, the MUSIC method does not stack the traces but separates the incident wavefield in a signal and a noise subspace. The eigenvectors of the signal subspace, which are normal to the noise subspace, contain information on the apparent velocity and backazimuth of the incident wavefield. Each eigenvector describes the features of an incident wavefield generated by a single source, so the number of eigenvectors corresponds to the number of recorded sources. The MU-SIC algorithm was first introduced in seismology by Goldstein & Archuleta (1987). Goldstein & Archuleta (1991) used the MUSIC method to determine the parameters of the rupture propagation during an earthquake. Bokelmann & Baisch (1999) located a narrow band signal recorded at the German Experimental Seismic System array. By advancing the MUSIC algorithm to three component data, Hobiger et al. (2016) estimated the backazimuth and velocity of the incident wavefields, and also estimated polarization parameters of Rayleigh waves. Lacroix et al. (2012) monitored avalanches in the French Alps using a seismic array and located the events using beam-forming. They visually identified 80 events, and were able to assign the avalanche events to known avalanche paths. By calculating a correlation coefficient, they filtered the event directions obtained by beam-forming. Based on these results they determined the direction of the source, i.e. the avalanche front, and using additional knowledge on past avalanche events, they were able to determine the path of the event. Furthermore, they derived the speed of the avalanche front and therefore obtained indirect information on avalanche type, as dry-snow avalanches are generally faster than wet-snow avalanches (Nishimura & Izumi 1997;Vilajosana et al. 2007).
Array analysis methods provide additional information on seismic sources and hence may provide further insight into avalanche characteristics. To perform a reliable array processing the exact positions of the sensors have to be known, which can be easily achieved using readily available GPS devices. However, the resolution of the array is limited by the array geometry (Wathelet et al. 2008;Maranò et al. 2014). When planning array geometry it is crucial to roughly know the parameters of the monitored signals in advance.
In terms of localization of mass movements, especially snow avalanches, the MUSIC algorithm has not been used yet. In this work, we therefore compare results from beam-forming and MUSIC applied to seismic signals generated by snow avalanches recorded at a small aperture geophone array (maximum distance between geophones ≤ 70 m).

F I E L D S I T E A N D I N S T RU M E N TAT I O N
We used a seismic monitoring system consisting of seven vertical geophones ( Fig. 1b) with an eigenfrequency of 4.5 Hz and sampled at a rate of 500 Hz (van Herwijnen & Schweizer 2011a). Data were recorded using a 24-bit acquisition system and a GPS device provided accurate timing. To increase the signal-to-noise ratio sensors were buried ∼ 50 cm deep as suggested by Heck et al. (2018a). The average interstation distance was 30 m and the maximum distance between two sensors in the array was 64 m. The resolution of the array was thus rather limited (Supporting Information Fig. S2). For a frequency band of 4 Hz, the array was only capable of resolving wavefields with velocities between 120 and 450 m s −1 . For increasing frequencies the range of resolvable wave-velocities increased, e.g. for 10 Hz to velocities between 300 and 1200 m s −1 .
The instrumentation was deployed at the Dürrboden field site at the end of the Dischma valley near Davos (Eastern Swiss Alps;46.27 • N,9.92 • E, see Fig. 1a) from 2016 November to 2017 May. The field site is a flat meadow at an elevation of 2000 m a.s.l. surrounded by mountain peaks that rise up to 3000 m. The field site is also equipped with several automatic weather stations measuring the snow height using a ultrasonic distance sensor and one station measured the precipitation using a heated precipitation sensor. The precipitation is measured in mm and by assuming a density of 100 kg m −3 for new snow, 1 mm corresponds to 1 cm of new snow. Furthermore, eight automatic cameras were installed at two different locations for visual monitoring of the adjacent slopes. Images were recorded every 10 min throughout the winter.

M E T H O D S
3.1 Seismic signal processing and identification of avalanche signals van Herwijnen & Schweizer (2011a) showed that continuous seismic data contain various types of events, such as airplanes, helicopters and avalanches. To limit the amount of data to process, only events in the seismic data with amplitudes higher than a threshold value of five times the daily mean amplitude were used (Heck et al. 2018a). To identify seismic signals generated by avalanches within these events, we then performed a visual inspection of the waveform and the spectrogram of the seismic events as proposed by van Herwijnen & Schweizer (2011b), Lacroix et al. (2012), van Herwijnen et al. (2016, Heck et al. (2018a). Although we installed automatic cameras to monitor a wide area around the field site (red polygons in Fig. 1), only two of the visually identified events in the seismic data for the whole season 2016-2017 could be confirmed as avalanches on the images of the cameras. Both events released in May. The other events could not be confirmed either because they occurred at night or during snowfall when visibility was poor.
Once we identified events in the seismic data, we visually picked the start and end times of the seismic signals. As seismic signals of avalanches do not have a sharp onset, it is difficult to determine the exact onset. However, the exact onset of the avalanche is not critical as it is more important to analyse the entire signal. We therefore allowed some noise before the onset and after the end around the obvious signals.

Beam-forming
Beam-forming is a commonly used technique in array processing where all the traces are stacked. The idea is that coherent phases in the signal are amplified whereas incoherent noise is reduced (Rost & Thomas 2002). Since the same signal arrives at different sensors of the array at different times, the seismic traces are shifted in time before stacking. By combining the time-shift resulting in the maximum stacked amplitude of the signals and the positions of the sensors, it is possible to calculate the apparent velocity of the seismic wave as well as the backazimuth of the signal. A seismic signal recorded at sensor i of the array is represented by a time-series x i (t): with f(t) the seismic time-series shifted for r i · s hor , where r i is the distance to an arbitrary reference point and s hor the horizontal slowness of the wavefield with and w i incoherent noise. The absolute value |s hor | = u is the slowness of the wavefield, which is the inverse of the velocity v hor = 1/u. The angle φ is the direction of propagation of the wavefield. The angle θ = 180 • + φ points towards the source of the signal, i.e. the backazimuth.
The time-series are shifted in a manner that coherent phases arrive at the same time at a reference station r = r 0 : The stacked trace for N sensors for an event j is To find the correct slowness vector s hor , the total energy for the stacked trace b j (t) is maximized. In our case, the slowness vector s hor is unknown. A 2-D grid search is therefore performed by varying the direction φ and the slowness, and thus the velocity of the incident wavefield. The values maximizing the energy of the resulting stacked trace define the incident wavefield. In our analysis we use the beam-forming method provided by the ObsPy toolbox using a window length of 1 s and no overlap (Beyreuther et al. 2010). By expecting at least three wave cycles within the window this restricts the analysis to frequency bands above 3 Hz. We analysed the beam-forming results for two different filtered data sets. First we analysed the data using a fourth-order Butterworth low-pass filter with a corner frequency of 50 Hz. In a second analysis we filtered the data using a fourth-order Butterworth bandpass filter with corner frequencies of 1.5 and 14 Hz.

Multiple signal classification (MUSIC)
The MUSIC algorithm, as introduced by Schmidt (1986), is a method to separate the signal and noise subspaces of incident wavefields. The signal X recorded at an arbitrary array consisting of n vertical geophones is a linear combination of noise (W) and the incident wavefields (AF) The elements of A are known functions of the signal arrival angle and the array location and F represent the amplitude and phase of the signal. The covariance matrix S of size n × n of the signal is computed by with * denoting the complex conjugate of a matrix. The eigenvectors and eigenvalues of the matrix S are calculated. The eigenvectors with eigenvalues exceeding a predefined dynamic threshold value span the signal subspace, whereas the remaining eigenvectors define the noise subspace G. Depending on the threshold value, the MUSIC algorithm can identify several signals with different backazimuth θ and incident velocity v, that is, the algorithm can distinguish multiple sources. For an incident wave, the MUSIC functional is with a(k) the steering vector, which represents the set of phasedelays at each sensor of the array for an incoming wave with wavenumber vector k. The wavenumber vector k is given by where s(f) describes the slowness of the propagating wave with frequency f and θ the backazimuth of the wave. The values s(f) and θ maximizing P MU are determined using a 2-D grid search.
As suggested by Hobiger et al. (2016), the MUSIC algorithm was applied to several frequency bands of the signal. We defined 38 central frequencies ranging from 3 to 21.5 Hz with an equidistant spacing of 0.5 Hz. The data were divided into 38 frequency bands with these central frequencies and a bandwidth of 0.3 Hz using a first-order Chebyshev filter. The MUSIC method was then applied to non-overlapping data windows for each frequency band. The length of a window was chosen such that five cycles of the lower corner of the analysed frequency band fitted. This means that the data were divided into more windows for higher frequencies and hence resulting in more backazimuth and slowness results. Doing so, we were able to investigate the frequency dependence of the backazimuth and apparent velocity. Due to the small array, we were limited in the resolvable velocities by array limits between k min = 0.111 and k max = 0.208. We therefore adjusted the minimal and maximal velocities for the grid search depending on the analysed frequency band.

P O S T -P RO C E S S I N G
We implemented a moving median filter to smoothen the backazimuth values obtained with MUSIC and beam-forming. To calculate the median of several angles, we assumed the backazimuth values as polar-coordinates on the unit circle and transformed them to the Cartesian system. For all backazimuth values within a sliding window of 8 s we performed this transformation and determined the median valuesx andỹ for the x and y component for this time window. By using the atan2 function the θ coordinate of φ is calculated, which is the median backazimuth for a time window. This step was performed for all time windows, resulting in more robust backazimuth values.

Winter season 2016-2017 and avalanche period in 2017 March
The winter season of 2016-2017 was relatively short and characterized by a below-average snow depth (Fig. 2a). The first snowfall in December was followed by a four week period without substantial precipitation. During this period, a strong temperature gradient prevailed within the snow cover. Due to the ongoing snow metamorphism, by the end of December 2016 the entire snow cover consisted of poorly bonded depth hoar crystals. While several small precipitation events in January and early February caused avalanching in the region of Davos (red bars in Fig. 2a), the total snow depth remained below 100 cm and the avalanche activity at our field site remained very modest. Several snow storms between March 6 and 10 deposited around 80 cm of new snow and resulted in a widespread avalanche cycle with avalanches stepping down into the weak depth hoar layer at the base of the snow cover. The visibility during this early March storm was poor, and it was not possible to obtain accurate avalanche release times from the images of the automatic cameras. Once the weather cleared the intensity of the avalanche cycle at our site became evident as the number of avalanches that had released was very high. During a field survey on 2017 March 15 we identified and mapped 24 avalanches in a 4 km radius of the seismic array (Fig. 2b). We are confident that all these avalanches released during the intense snowfall period, as older avalanches would have been covered by the large snowfall.
As we did not know the accurate release time of the avalanches, we visually inspected the seismic data to identify avalanches (see Section 3.1). We thus identified 19 seismic events within a 24 hr period starting on 2017 March 9. In Table 1 we classified these seismic events according to the event duration in the seismic data, which provides a rough estimate of the avalanche size (van Herwijnen et al. 2013;Pérez-Guillén et al. 2016). Of particular interest is an event on 2017 March 10 at 03:47, corresponding to an avalanche that released directly above the seismic array, approaching to a distance of less than 100 m. Based on the characteristics of the seismic event we were confident of linking it with the observed avalanche.

Avalanche on 2017 March 10 at 03:47
The avalanche that released on 2017 March 10 at 03:47 generated a seismic signal with a duration of approximately 100 s and can be considered a medium to large avalanche (Fig. 3). It released at an elevation of 2600 m a.s.l near the top of Chlein Sattelhorn, travelled a distance of approximately 1.2 km and came to a halt approximately 80 m from the seismic array (Fig. 3d). As the avalanche moved towards the seismic array, the seismic signals became more energetic and the high frequency content (> 30 Hz) considerably increased, in particular after 40 s.

Path reconstruction applying MUSIC
We first applied the MUSIC algorithm and analysed several frequency bands. In Fig. 4 the percentage of backazimuth angles for a maximum of two sources resolved by the MUSIC algorithm are shown for each analysed frequency band. For the frequency bands between 3 and 12 Hz two backazimuth angles dominated at 210 • and at 60 • . Between 12 and 17 Hz these two dominant directions shifted to 50 • and 230 • . These two directions correspond to the directions of the sides of the valley and were found in the analysis of several avalanche events. For frequency bands above 17 Hz, there were no clear maxima as scattering dominated wave propagation resulting in more evenly distributed backazimuth angles.
Since the path of the avalanche was known (Fig. 3), backazimuth values lower than 150 • and higher than 210 • were not associated with the event. Frequency bands between 3 and 12 Hz therefore contained information about the event source, although there also was a second source.
For a central frequency of 3.5 Hz, the source for the strongest signal dominated (red dots in Figs 5a and b), a secondary source generated by weaker sources was almost not resolved (low number of grey asterisks in Figs 5a and b). For the 6.5 Hz frequency band, a secondary source became more evident (grey asterisks in Figs 5c and d).
The backazimuths for the 3.5 and 6.5 Hz frequency bands started at 210 • and decreased with time towards 160 • (Figs 5a and c). In the 3.5 Hz frequency band the decrease was consistent and showed little variation. However, at higher frequencies strong variations in backazimuth occurred between 30 s and 55 s.
To determine the most probable path of the avalanche, we combined the backazimuth angles for the dominant source of all analysed bands with a central frequency ≤ 12 Hz and used a polar plot  to better visualize changes in the backazimuth (Fig. 6). In a polar plot, the radius represents time with 0 s in the middle; the angle represents the backazimuth and therefore points towards the source. By plotting the results in such a manner, changes of the backazimuth become more clearly visible (compare Fig. 5 with Fig. 6). We also calculated the median filtered backazimuth of the avalanche (blue line in Fig. 6). Initially, the backazimuth angle was around 210 • and slightly decreased with time. After about 50 s, the angle sharply changed towards 150 • and around 80 s the backazimuth became more randomly distributed suggesting that only noise was recorded. Additionally, a second direction given by the second most dominant source emerged with backazimuth angles between 45 • and 60 • between 0 and 70 s. Changes of these backazimuth angles correspond with changes of the median backazimuth path.
The apparent velocity of the incoming wavefield for the most dominant source did not vary much for the different frequency bands and was generally between 300 and 400 m s −1 (Figs 5b and  d). In the beginning of the signal, the velocities were somewhat higher and decreased with time. For the second source, however, velocities were more variable.

Path reconstruction applying beam-forming
We also reconstructed the path of the avalanche using the beamforming approach (Section 3.2). Before we analysed the data using the beam-forming algorithm, we filtered the data using a fourthorder Butterworth low-pass filter with a corner frequency of 50 Hz. However, as shown above, the most promising frequencies were between 3 and 12 Hz. We therefore also applied a fourth-order Butterworth bandpass filter with corner frequencies of 1.5 and 14 Hz.
For the low-pass filtered data, a constant angle at 210 • was visible for the first 20 s. Afterwards variations in backazimuth became larger as time increased (Fig. 7a). For t > 80 s variations in backazimuth became even stronger suggesting that noise was recorded. The same behaviour of the backazimuth was also observed for the     bandpass filtered data. However, the variations were not as strong. For both filter methods the velocities were around 330 m s −1 for the first 20 s and between 60 and 80 s and more erratic for the rest of the signal.
Since the variations of the backazimuth were smaller for the bandpass filtered data, these data were used for further analysis. The trend in the median filtered backazimuth was similar to the MUSIC results (Fig. 8). However, due to the strong variations, the reconstructed path of the avalanche was less accurate and only the initial direction was well resolved. Changes in the backazimuth corresponding to a source approaching the array were not well observed.

Apparent velocity
In both the MUSIC and beam-forming approaches the apparent velocity of the incident wavefield was estimated.
The apparent velocity corresponds to the velocity of incoming wavefields in the sensor plane for an incident angle of 0 • , but increases as soon as the waves have an incident angle >0 • . As it is unclear if seismic or acoustic waves were measured by the array, we therefore analysed the velocities in more detail. Since we were only interested in the signals generated by the avalanche, we only kept signals with a backazimuth angle within 10 per cent of the calculated median path of the event (blue lines in Figs 6 and 8).
While there was a large spread in the velocity values obtained with MUSIC, the majority of the values were around 330 m s −1 , which is the speed of sound in air at 0 • C (red line in Fig. 9a). Two clusters were clearly present, approximately the first 30 s and from 55 to 80 s. For both these clusters, the velocities were around the speed of sound, with slightly higher values for the first cluster. For the period in between the clusters, from 30 to 55 s, the velocities were between 330 and 1500 m s −1 , the maximum value for the grid search.
The velocity values obtained using the beam-forming method were similar and both clusters were also present in the results, although less clear (Fig. 9b). For the first cluster the values were higher than the speed of sound and for the second cluster around 330 m s −1 . Apart from these two clusters, the velocities were higher, between 1000 and 2000 m s −1 .

Reconstructing the median filtered path
The rough direction of both array processing techniques was the same. However, due to the strong variations of the median backazimuth obtained with beam-forming (Fig. 8), the source was not as well defined as with the MUSIC method. We therefore analysed the resulting median event path using the MUSIC method and projected it onto a map.
We divided the time-series into four sections as indicated in Fig. 10(a), depending on the energy received at the sensors. The first section corresponds to the beginning of the signal with a slight increase in energy from 10 to 35 s. Then there was a short section of less than 5 s when almost no energy was transmitted to the array. For the third section, from 42 to 57 s, the energy increased again. The fourth section contains the highest amplitude signals and from 57 s to the end of the signal at 80 s most energy arrived at the array.
The signal in the first section had an almost constant backazimuth angle of about 210 • and an apparent velocity for the incoming wavefields around 450 m s −1 (Figs 10b and c). A projection of this almost constant direction onto the map corresponds to the slopes north of the Chlein Sattelhorn (Fig. 10d). As the avalanche moved towards the array, the backazimuth remained about the same for the first 30 s.
The approximate direction of the second section corresponds to the first terrain break the avalanche was flowing over (2 in Fig. 10d). Shortly after this second section, the energy increased again and in the steeper terrain of Section 3, the backazimuth as well as the apparent velocity could be determined (3 in Figs 10b and c).
The fourth section was characterized by strong signals providing a clear backazimuth and apparent velocities (4 in Figs 10b and c). The median backazimuth showed a clear change from 190 • to 160 • indicating that the front was passing close to the array. The apparent velocities were around 330 m s −1 suggesting that the source had reached the same elevation as the array.
Although we observed strong variations for the backazimuths, the actual path of the avalanche was an almost straight line. The strong changes in the backazimuth were a consequence of the mutual position of the moving front of the avalanche and the array.
The visually observed avalanche path (Fig. 3d) corresponded well with the reconstructed event path obtained with MUSIC. From the point of release of the avalanche around 10 s to the arrival at the valley bottom at around 70 s, the front of the avalanche travelled about 1100 m in 60 s. Although the avalanche front came to a halt, the tail of the avalanche continued to generate signals. A duration of 60 s corresponds to an average speed of the avalanche front of approximately 70 km hr −1 .

F U RT H E R E V E N T S O N 2 0 1 7 M A RC H 9 A N D 1 0
We analysed the remaining 18 avalanches, which we had visually identified in the seismic data from March 9 to 10 and investigated if an obvious path was visible in the polar plots. For 4 of the 18 events, there was no obvious path and the events were probably falsely identified as an avalanche. It is likely that those events were generated by environmental noise such as airplanes. For instance, the median filtered path for the event on 2017 March 9 at 18:48 is undefined (Fig. 11a). However, visually several backazimuth directions were identified, especially between 60 • and 90 • . Based on the visual analysis, we assumed that we observed an avalanche coming from this direction.
For the remaining 14 events, a median path was clearly visible, using the same procedure as described above (Fig. 12).
For all of the analysed events, except the avalanche that released on 2017 March 10 at 00:28 (Fig. 12b), directions pointing towards the side of the valley, between 45 • and 60 • and between 220 • and 240 • , were also visible. In addition for the events shown in Figs 11(b) and 12(a), an additional artefact was observed.
To relate the seismic events to avalanches observed in the field, we determined a median path for all mapped avalanches from our field survey. We then searched for all intersections of these median paths with the line starting at the middle point of the array with the direction of the backazimuth for each event. This allowed us to automatically match seven events with avalanches observed in our field survey.
By visually comparing the remaining events with the mapped avalanches, we were able to relate an additional four events with avalanches. Thus, in total we matched 11 of the 14 events with visually observed avalanches (Fig. 13a). The remaining three events did  not match any observed avalanche in the vicinity of the array. Nevertheless, we assumed that these events corresponded to avalanches not mapped during our field survey as the median paths were clear and did not point to the sides of the valley.
It is well known that most natural dry-snow avalanches release during or directly after a snow storm Perla & Martinelli (1978, p. 35). However, little is known about the exact timing. We therefore compared the amount of fallen snow with the release times of the   identified avalanches (Fig. 13b). The first avalanches released after only 20 cm of snow had fallen and avalanche activity at our site continued to the end of the snowfall. Once the snowfall stopped, no further avalanches were identified at our field site.

W E T -S N O W AVA L A N C H E S I N M AY
The avalanches during the period in March were mostly dry-snow avalanches. The apparent velocities of the incoming wavefields for these event were all around 330 m s −1 , corresponding to the speed of sound. To investigate the localization of wet-snow avalanches, we searched for further events in April and May. However the avalanche activity during these months was very low (Fig. 2a) and we only identified two wet-snow avalanches on the images from the cameras on 2017 May 6 at 14:13 and on 2016 May 14 at 19:35 (red areas in Fig. 2b). Furthermore, we identified a second potential avalanche in the seismic data on 2017 April 5 at 5:41. We analysed these two avalanches using the MUSIC method.
According to the images of the automatic cameras (larger red area in Fig. 2b), the backazimuth of the event on May 6 should start at around 270 • and end around 320 • . However, the median path of the event as shown in Fig. 14(a) had a direction around 320 • with only slight changes in backazimuth. This direction corresponds to the opening of the valley towards the city of Davos. Furthermore, there were strong variations in the apparent velocity and no clear trend was visible. Although we identified a clear signal in the seismic data, it was not possible to clearly locate the event.
This location of the avalanche was known in advance and it was unambiguously identified as a wet-snow avalanche on the pictures of an automatic camera since it released to the ground. Such events are expected to generate strong seismic signals, and indeed clear waveforms were recorded by the array (Fig. 14c). Nevertheless, it was not possible to determine a backazimuth corresponding to the visually observed avalanche that was between 270 • and 320 • . Even more detailed analyses using beam-forming and MUSIC with a wider frequency band did not reveal any significant backazimuth corresponding to the event (Supporting Information Fig. S2). In addition, we performed an fk analysis for this event to identify wavefields (Supporting Information Fig. S4). Apart from a fundamental mode, as also visible for the avalanche event of 2017 March 10 on 03:47 shown in Supporting Information Fig. S3, no clear wavefields were visible. We assume that this event mainly generated seismic signals that could not be resolved by our array due to the limited interstation distance. Hence, the amplitude of the infrasound signals was very low and signals were likely to be superposed by environmental noise coming from Davos making it impossible for the MUSIC method to distinguish between two different types of source.
A second wet-snow avalanche was monitored by the automatic cameras on 2017 May 14 at 19:35 (smaller red area in Fig. 2 and Fig. 15). Similar to the first observed wet-snow avalanche in May, it was not possible to resolve a median backazimuth path, which corresponded to the visually observed avalanche path, although a decrease in the apparent velocity between 40 and 60 s was observed.

D I S C U S S I O N
In this study we used two seismic array processing methods, namely MUSIC and beam-forming, to calculate the backazimuth angles as well as the apparent velocity of incident wavefields generated by snow avalanches. An analysis of a representative avalanche event that released close to our seismic array on 2017 March 10 at 03:47 revealed that the most robust results for the localization were obtained by bandpass filtering the data between 3.5 and 12 Hz, regardless of the method (Figs 5 and 7b). However, results obtained with the beam-forming method showed strong variations for the backazimuth angles preventing an unambiguous reconstruction of the avalanche path, while the results obtained with the MUSIC algorithm were more robust. In addition to the backazimuth angles, we also obtained values for the apparent velocity of the incoming wavefields, which were mostly close to the speed of sound in air (Fig. 9).
Overall, our results suggest that with our seismic array, we mostly resolved infrasound signals generated by the avalanche. While these results are somewhat surprising, this is mainly due to limitations related to our deployment geometry (Fig. 1b). Indeed, the distance between adjacent sensors was 30 m, resulting in a rather limited maximum distance between sensors of 64 m. For a frequency band of 4 Hz, the array could only resolve wavefields travelling at a speed between 120 and 450 ms −1 . For a frequency band of 10 Hz, on the other hand, waves travelling between 300 and 1200 ms −1 could be resolved. Our array was therefore very well suited to resolve infrasonic wavefields, while seismic waves, which generally propagate at velocities above 1200 ms −1 can only poorly be resolved at low frequencies.
Although the installed array showed limitations, the MUSIC algorithm was capable of resolving the infrasonic wavefield and estimate the backazimuth angle of the source whereas the beamforming method worked less well. Our results are in contradiction with those from Lacroix et al. (2012), who successfully used the beam-forming method to locate avalanches with a seismic array in the French Alps. However, their array was substantially larger as the minimum distance between the sensors was around 45 m, resulting in an array with a better resolution. Based on the backazimuth results, Lacroix et al. (2012) were able to assign signals to known avalanche paths and by combining the information of the path, backazimuth and time, they were able to estimate the velocity of the avalanche front. Based on visual observations obtained during a field survey, we were able to link certain seismic events with specific avalanche paths, which also allowed us to estimate the velocity of two avalanche events.
The results obtained with the MUSIC method showed that for many avalanches there was also a strong secondary source, with backazimuth angles either between 45 • and 60 • or between 220 • and 240 • (Fig. 4). An fk spectrum analysis of an avalanche event (Supporting Information Fig. S1) showed that these directions were likely associated with a fundamental Rayleigh wave mode stimulated by the seismic events. Indeed, for frequencies between 12 and 17 Hz incoming wavefields travelled at a speed between 400 and 1200 ms −1 , whereas the speed decreased with increasing frequency. These signals were not observed for noise. The backazimuth angles associated with these waves were either pointing towards the northeast between 45 • and 60 • or were shifted by almost 180 • pointing to the south-west between 220 • and 235 • . These directions correspond to the sides of the valley. The analysis of the avalanche events was therefore performed within a range of 3.5 to 12 Hz to reduce the influence of this mode.
The avalanche period in March, when most of the avalanches were unambiguously identified, mostly dry-snow avalanches were monitored and it was possible to determine the source direction using the MUSIC method. For several wet-snow avalanche in May, which were visually identified in the pictures, it was not possible to estimate the source direction using the MUSIC method.
Based on images from the automatic cameras, the release time of these avalanches was well constrained. Furthermore, there were clear signals in the seismic data. We therefore assume that these wet-snow avalanches mostly generated seismic signals and did not generate energetic infrasound signals. Given the limited resolution of our seismic array, we could therefore not identify a clear backazimuth direction. These results confirm that the geometry of our seismic array was indeed not well suited to resolve seismic waves. Furthermore, due to the rapid ablation of the snow cover in May (Fig. 2a), the overall signal-to-noise ratio decreased, further complicating the identification of a clear backazimuth direction, in particular at higher frequencies where our array could potentially resolve seismic wavefields.
The main goal of this work was to detect the seismic signals of avalanches using a geophone array and applying array processing techniques to locate the avalanche path. In order to reduce environmental noise as much as possible, the sensors were buried 50 cm deep in the ground, as suggested by Heck et al. (2018a). We first manually identified 19 probable avalanche events in the continuous seismic data during a high avalanche period from 2017 March 9 and 10, and analysed these events using the MUSIC method. Using the backazimuth results, we were able to reject four of these events as these were not associated with a clear direction (e.g. Fig. 11a) ). Furthermore, 12 of the events were linked to avalanches observed during a field survey (Fig. 13a), suggesting that results from array-based source localization could be a valuable for the automatic detection of avalanches, as is commonly done for infrasound avalanche detection (Marchetti et al. 2015). However, as shown by Heck et al. (2018b), the MUSIC method is not appropriate as a standalone automatic classification technique as the computational costs are too high. Nevertheless, this method is suited as an additional tool to confirm or deny already detected events with different automatic classification methods.
With the advent of automatic classification method for continuous seismic data, obtaining accurate information on avalanches is becoming more tangible (Hammer et al. 2017;Heck et al. 2018a). Such data are very valuable to improve our understanding of avalanche formation processes (e.g. van Herwijnen et al. 2016). As mentioned by Crouzy et al. (2014), the synchronization between snowfall and avalanches still remains a fundamental problem that is difficult to solve with a stochastic model. They suggested a model to predict the frequency of small avalanches depending on the slope angle and precipitation at a well-monitored avalanche area in the Vallée de la Sionne in Western Switzerland. To apply this model at different locations, avalanche activity data containing information of the exact release time and location of the avalanches are required. In our study we compared the release times of avalanches with the observed precipitation measured at a nearby automatic weather station to investigate the timing of avalanche release during the snow storm. Results obtained during such an analysis are valuable data for stochastic models. Based on this analysis we observed that the first avalanches released after approximately 20 cm of new snow had fallen. Thereafter, the activity was relatively equally distributed throughout the precipitation event, that is, there was no clear peak in the avalanche activity. Finally, as soon as the precipitation ceased, no more avalanches were observed. During a major snowfall, Lacroix et al. (2012) also reported similar results. While these results provide new insight into the timing of avalanches during a snow storm, more data are required to confirm whether such avalanche activity patterns are generally observed.

S U M M A RY A N D O U T L O O K
We used two different array processing techniques to reconstruct the path of 19 avalanche events recorded by a seismic array during a period of high avalanche activity on 2017 March 9 and 10. Both methods, the MUSIC as well as the beam-forming approach, allow calculating the backazimuth of the source and the apparent velocity of the incident wavefields.
The MUSIC algorithm provided more reliable results with less variations for the backazimuth angles than the beam-forming method. A more detailed analyses of an unambiguous avalanche recorded on 2017 March 10 at 03:47 showed, that based on the MUSIC results a reconstruction of the avalanche path was possible. Using the beam-forming results, however, only a rough direction could be estimated. By combining results for the 19 events using the MUSIC method with observations from a field survey performed several days after the avalanche period, we were able to match 12 events with observed avalanches. These results show that the detection range of our seismic monitoring system was approximately 2-3 km. The obtained velocities of the incident wavefields revealed that with our array we mainly resolved infrasonic waves rather than seismic waves. This was attributed to the specific array geometry, as the maximum distance between sensors was 64 m. Finally, we investigated two wet-snow avalanches that released in May, showing no clear backazimuth direction. These results suggest that such avalanches mostly generated seismic signals, which could not be resolved with our array.
To improve the localization of snow avalanches through seismic monitoring, we intend to modify our instrumentation to increase the overall aperture of the seismic array, as used by Lacroix et al. (2012). Moreover, using additional arrays could improve the information on the distance of the events as triangulation methods could be applied (e.g. Métaxian et al. 2002). Based on the localization using the MUSIC method, it is possible to determine whether a seismic event has randomly distributed backazimuth angles or a median path is visible. Therefore the localization of avalanches seems a helpful tool for the classification of seismic events and we plan to implement it into our previously developed automatic classification approach (Heck et al. 2018b).

A C K N O W L E D G E M E N T S
M.H. was supported by a grant of the Swiss National Science Foundation (200021 149329). We thank numerous colleagues from SLF for help with field work and maintaining the instrumentation.

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. Figure S1. Results of the fk analysis. The spectral density is indicated by different colours. Signals within the solid and dashed line can be resolved. Signals below the dashed line and above the straight line are due to aliasing. Figure S2. Resolution of the array used by Lacroix et al. (2012) (red) compared to the resolution of our array. Figure S3. Dominant backazimuth values for the avalanche released on 2017 May 6 at 14:13 for different frequency bands obtained with the MUSIC method. Figure S4. Results of the fk analysis for the avalanche released on 2017 May 6 at 14:13. The spectral density is indicated by different colours. Signals within the solid and dashed line can be resolved. Signals below the dashed line and above the straight line are due to aliasing. 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.