A benchmark case study for seismic event relative location

, ray from the source region. We demonstrate that, had such corrections been available, a signiﬁcant improvement in the relative location estimates would have resulted. In practice we would of course need to solve for event hypocentres and slowness corrections simultaneously, and signiﬁcant work will be needed to upgrade relative location algorithms to accommodate uncertainty in the form of the outgoing waveﬁeld. We present this data set, together with GT coordinates, raw waveforms for all events on six regional stations, and tables of time-delay measurements, as a reference benchmark by which relative location algorithms and software can be evaluated.

greater improvement is possible using accurately measured relative times (e.g. Lin et al. 2007). The time-evolution and operational significance of induced seismicity can also be revealed using precision seismology (e.g. Goertz-Allmann et al. 2017;Yu et al. 2019). In forensic seismology, accurate relative location estimates of underground nuclear tests (e.g. Fisk 2002; Waldhauser et al. 2004) can place constraints on emplacement and depth with possible consequences for yield estimation (e.g. Pasyanos & Myers 2018;Voytan et al. 2019).
There are sources of error in relative event location which are often overlooked by practitioners. First is the assumption that the time of maximum correlation between two waveforms is an accurate indicator of the time separating two equivalent segments of the wave trains from different events. Simply displaying the signals from the 3 September 2007 North Korean underground nuclear test together with the signals from the significantly smaller previous explosions at the same site made it clear that a correlation-based alignment of waveform features following the signal onset could result in qualitatively misleading estimates of the relative onset-time (Gibbons et al. 2018b). The limitations of such estimates were discussed in greater depth by Bachura & Fischer (2019) for earthquakes in Czechia. Secondly is the assumed mapping between a relative time-delay and an interevent distance. Although a bias in the traveltime model may be cancelled out by the simultaneous location of two events, there is an underlying assumption of how far a P wave, S wave, or surface wave has travelled in a given time. Selby (2010) and Wen & Long (2010) had both calculated the location of the 25 May 2009 North Korean test relative to the 9 October 2006 test and found rather different locations from different sets of measurements: one teleseismic and one regional. Gibbons et al. (2017b) argued that the outgoing seismic wavefield from the test-site must be rather more complicated than could be explained using a 1-D velocity model and estimated an azimuthally dependent set of slowness corrections which made the regional and teleseismic estimates more internally consistent. The slowness corrections obtained appeared to vary approximately sinusoidally with azimuth, as would be observed given a dipping layer beneath the test-site. However, we cannot confirm a dipping layer hypothesis on the basis of these observations since the formulation of Gibbons et al. (2017b) only allowed the slowness, and not the backazimuth, to vary.
There are few situations in which the results of double-difference location estimates can be confirmed absolutely. Cross-validation studies can be performed by using different subsets of measurements (e.g. Murphy et al. 2013;Gibbons et al. 2017b) or, for example, body wave relative location estimates can be compared with surface wave relative location estimates (e.g. Ichinose et al. 2017). However, ultimately, assumptions about the propagation of seismic waves will have been made. In this paper we present a set of 55 surface explosions carried out at Hukkakero in Finland in August and September 2007. The explosions all take place within 300 m of one another and all were well recorded by many seismic stations within a distance of 2 • . The exact coordinates of the explosions were provided by the Finnish military, allowing for an absolute check on quality of double-difference location estimates.
We select a set of six seismic stations: three permanent and three temporary. The temporary stations were part of the LAP-NET/POLENET deployment as part of the International Polar Year in 2007 and 2008 (e.g. Kozlovskaya et al. 2016;Usoltseva & Kozlovskaya 2016;Kozlovskaya 2007). Raw seismic waveforms from these stations are provided in miniSEED format for all events in the Supporting Information. In addition, cross-correlation calculations have been performed for both P and S arrivals on all six stations for every pair of events. The results of these calculations, with delaytime estimates and correlation coefficient values, are provided in a single ASCII file in the Supporting Information. Details of the Supporting Information are provided in the Appendix. With the raw waveforms, readers are encouraged to test relative event location algorithms and software through measurement of the delay times and subsequent location estimates. We describe the explosions and the corresponding set of seismic signals and discuss some of the issues involved in our attempts to locate the events relative to each other.

T H E G RO U N D -T RU T H DATA S E T
For many decades, the Finnish military have disposed of expired ammunition at a site in Central Lapland (northern Finland) in a sequence of controlled surface explosions between August and September every year. These events, referred to as the Hukkakero explosions, have received a lot of attention due to the fact that they generate strong infrasound signals recorded at distances up to several thousand kilometres (e.g. Gibbons et al. 2018a). It was demonstrated by Gibbons et al. (2007) that the seismic signals generated were essentially identical from explosion to explosion, meaning that the significant differences in the infrasound signals from event to event were due entirely to changes in the atmospheric conditions. Accurate measurements of backazimuth and phase velocity could be made of the infrasonic signals at distances of up to several hundred kilometres and Blixt et al. (2019) demonstrated that the azimuthal deviation in these measurements could be used to validate cross-wind predictions from atmospheric models.
Whereas the events have been of greatest interest for the infrasound signals, the repeatability of the seismic signals has been exploited for studies of the stability of direction estimates using f-k analysis on seismic arrays. Gibbons et al. (2010) measured the direction of arrival (DOA) of the initial P-wave arrival at the ARCES seismic array for many Hukkakero events in different frequency bands. It could be demonstrated that the frequency band in which the most stable backazimuth estimates were made was the 2-4 Hz band. In bands of increasing dominant frequency, the signal-to-noise ratio (SNR) improved while the stability of backazimuth estimates became worse. This was the exact opposite to what was observed for several other more distant sites, for which both the stability and the SNR were optimal at a higher frequency. It should be noted that at ARCES, the Pg, Pb and Pn phases (see Storchak et al. 2003) from Hukkakero events arrive essentially simultaneously. These phases are associated with slightly different slowness vectors and dominant frequencies and this may also have consequences for the relative event location procedure.
Most Hukkakero events have an explosive yield of between 15 and 30 tonnes TNT equivalent, although a number of far lower yield explosions have also been associated with the same site (Gibbons et al. 2018a). The explosives are raised above ground level and are not well-coupled to the bedrock. This limits damage to ground, such that the same sites can be reused repeatedly. This is the reason that so much of the energy in the explosions is converted to infrasound, with the seismic source resulting from the downward shock wave. Table 1 gives the coordinates of 55 explosions from 2007 as provided by the Finnish military. We do not have exact values for the elevation at each location but differences in elevation from site to site are small compared with the lateral distances. Entering the coordinates into the satellite map system https://www.bing.com/maps (last accessed 16 June 2020) allows the reader to identify visible ground Downloaded from https://academic.oup.com/gji/article/223/2/1313/5879484 by guest on 21 January 2021  Fig. 1 shows the location of the explosion site in relation to the six selected seismic stations. Although other stations were and are operating in this region, we deemed it desirable to present a realistically sparse regional network. The azimuthal coverage of the source region is good, although the distances range from about 0.5 • (with the crustal Pg-type phases anticipated) to about 2.0 • (with the Moho and uppermost mantle Pn-type phases anticipated). We saw little value in including stations at greater distances as these are often unlikely to record events of magnitudes of interest when monitoring induced seismicity or earthquake swarms. We also avoided adding more stations at short distances as this may give the impression of applicability only in the case of far denser networks than are likely to be found in regions of interest. LP34, LP53, and LP61 were among the first stations to be deployed in the temporary LAPNET deployment. All three components of each station are displayed in Fig. 2 for event H19, which took place at 07:30 UT on 24 August 2007. The stations are ordered by epicentral distance and it is possible to see the waveform morphology at each of the stations. At the shortest distance (LP53 at 59 km), only 7 s of P arrival and coda are available before the S-wave arrival. At ARCES and KEV, over 20 s of P-wave energy is observed before the S arrival. The length of window to be used for time-delay calculations in doubledifference locations is a theme which needs careful consideration. The desire for a long time-bandwidth product for an accurate delaytime measurement must be balanced against the understanding that the phase velocity of the coda waves decreases with time after the signal onset. Ideally, only a few cycles at the start of the signal would be used for the correlation calculation. In practice, a somewhat longer window is likely to make a clearer and less ambiguous peak in the cross-correlation trace. All calculations performed here used a template with length 6.0 s and a frequency band 2-5 Hz. Fig. 3(a) displays the locations of the individual explosions. They take place in three or four clusters in a near ring with a diameter of between 250 and 300 m. Since the exact origin times of the explosions were not available in the ground truth (GT) information provided by the Finnish military, a correlation detector (using the formulation of Gibbons & Ringdal 2006) was run. A full-waveform signal for event H04 (2007-08-16T12.00.00) on the ARCES array was selected as a template and correlated against continuous ARCES array data for all of August and September 2007. This was done to confirm the presence of explosions at all the times indicated, at only the times indicated, and also to refine the origin time estimates. Note that the origin time estimates obtained from this detector assume that all events are exactly colocated and so will be inaccurate by the short differences in traveltime over the extent of the explosion site. In addition to the time of maximum correlation for each detection, the detector writes out the maximum value of the correlation coefficient at this time. The location symbols in Fig. 3 are coloured according to these values, providing a visual validation of the coordinates provided. The highest values of the correlation coefficients obtained are for the explosions set off closest to the Reference Event. The full waveform correlation values diminish with distance from this location. Fig. 3(b) displays waveforms at ARCES from four events as indicated from different clusters. It is clear how the waveform similarity is higher for events in the same cluster than for events in different clusters.

S E I S M I C R E L AT I V E L O C AT I O N E S T I M AT E S
As a first step to attempting a relative location calculation for these events, manual arrival-time estimates were made for both first P and first S arrivals for a single event (H04) on each of the available stations: 12 onset-time estimates in total. For each station, the P-arrival time was estimated on the vertical component trace filtered into many different frequency bands. The corresponding S-wave arrival time was made on the transverse rotation of the horizontal component traces (i.e. the ground motion perpendicular to the great-circle path connecting source and station), also observed in many different bands simultaneously. For the ARCES array, a beam was formed for both P and S arrivals using steering-parameters chosen to optimize the SNR for each arrival. These manually picked times for event H04 were used as the starting times for provisional waveform templates, and a set of reference times was obtained for the same phases for all other events using a simple time-domain correlation detector on the appropriate waveform segments. All waveforms were upsampled to 1000 Hz in order to perform the time-domain crosscorrelation calculation; this is not necessary for signal detection but mitigates a source of uncertainty in the time-delay estimates. The upsampling was performed using a time-domain sinc-interpolation (Schanze 1995). With provisional reference P and S times for all events on all stations, an exhaustive all-versus-all correlation calculation was performed in order to measure all necessary relative time delays. All of these correlation results are provided in the file absolute CC times.txt in the Supporting Information.
We follow the procedure described by Gibbons et al. (2017b) for calculating relative location estimates. This algorithm considers only the differences between pairs of arrival-time differences, eliminating the need to solve for the origin time. In addition to the time-delay measurement for each phase, we need a model prediction for the traveltime for that phase from any given source location to the station. For the initial calculation, we use two different 1-D models: the AK135 global model (Kennett et al. 1995) and the Fescan, or Fennoscandian, model (Mykkeltveit & Ringdal 1981).
The depth is held fixed in all calculations performed here. While there is no difficulty in scanning hypocentres at different depths, the traveltimes to stations at regional distances provide very little vertical resolution (see also Pasyanos & Myers 2018). Moving an event hypocentre laterally in one direction decreases the travel time to stations in that direction and increases the travel times to stations in the opposite direction; moving an event hypocentre vertically increases or decreases the travel times to all stations approximately equally. In this study, where all explosions are carried out on the surface, the fixed depth assumption is approximately correct. This may not be the case for natural or induced seismicity, although a fixed depth assumption is still recommended unless data is available which actually provides resolution in depth differences. Fig. 4 displays contour plots of the double-difference misfit residual norms for six different calculations with different velocity models and sets of measurements used. In each panel, the fixed location of the reference event (H03) is displayed and the (assumed unknown) location of event H02 is assumed to be at the location for which the double-difference misfit residual norm is at a minimum. In each of the calculations, the residual minimum (indicated by the red square) coincides quite well with the GT location (indicated by the white circle). The results for the two different models, AK135 in panels (a), (c) and (e) and fescan in panels (b), (d) and (f), are very similar. While the crustal and Moho velocities in fescan are rather faster than in AK135, the favourable azimuthal coverage means that the effect of a shorter or longer distance indicated in one direction is approximately balanced by a similar effect in the opposite direction. One of the most important visual lessons to be learned from Fig. 4 is the spatial resolution for different combinations of phases.
The residual contour lines for the (slower) S-phases (panels e and f) are denser than for the (faster) P-phases (panels c and d). The calculations with both P-and S-phases (panels a and b) show a density of the contour lines somewhere between those for the P-only and S-only calculations. The temptation is to claim that it would be best to use S-phases in isolation for improved spatial resolution. In practice, the S-phase may be poorly observed or competing with significant P-wave coda energy. While the location of the residual minimum is approximately in the correct location, the uncertainty is significant relative to the interevent distance.
We can perform the same calculation displayed in Fig. 4 with all of the events in the GT dataset. That is to say that we hold the location of event H03 fixed and attempt to locate the other 54 events relative using the misfit norm minimization technique. Fig. 5 displays the relative location estimates together with the events' GT locations and the corresponding mislocation vectors. For all phase combinations and both 1-D velocity models, the length of the mislocation vectors varies from close to zero to over 50 m (over 20 per cent of the true interevent distances). The mislocation vectors are encouragingly consistent for neighbouring events, providing confidence that the mislocations are due to fundamental model assumptions about the outgoing wavefield rather than, for example errors in measurement of the time-delays. The mislocation vectors for the northernmost cluster (in which the reference event H03 is located) are essentially zero, give or take noise within the region of uncertainty. The spread here is likely due to measurement errors, possibly from low SNR signals. For the other clusters of explosions, the mislocation vectors generally point towards the reference event. This means that the interevent distances are underestimated, and that the velocity with which (at least some) phases have left the explosion site have been underestimated.
In this section we have demonstrated that, using the available waveforms and simple time-domain cross-correlation time-delay estimates, we have largely managed to reconstruct the large-scale structure of the event cluster using a simple residual minimization algorithm. While the spatial resolution appears coarse from the contour plots, the locations of the minimum misfit are quite stable. The relationship between the absolute levels of the misfit residual norms and size of the regions enclosed by each contour is a clue as to how stable the calculation would be to a given time-measurement error, for example due to a low SNR signal or a significant waveform dissimilarity introduced by a difference in the source. The consistency within clusters of the mislocation vectors is an indicator of robustness of the procedure itself and identifies assumptions about the velocity of the phases leaving the explosion site as the most probable source of error.

Gibbons et al. (2017b) sought the relative locations of the declared underground nuclear tests carried out under Mount Mantap in North
Korea. They recognized that errors in the assumed velocities of the outgoing wave fronts were the likely source of discrepancy between earlier relative location estimates. An empirical algorithm was devised for estimating correction terms for each of the predicted traveltime differences. The correction terms obtained were subject to considerable uncertainty; the time-differences themselves are tiny fractions of a second and the correction terms estimate marginal adjustments to these delays. In addition, a dipping layer below the test-site would result in a deviation in both the apparent velocity and the direction of the outgoing waves; any azimuthal deviation in an outgoing wavefront was neglected. In spite of these difficulties, a set of corrections was obtained which greatly improved the consistency between estimates made using different sets of seismic phases. The obtained corrections showed a near sinusoidal azimuthal variability resembling that which would be observed given a dipping layer beneath the source.      Table 1 relative to the fixed event H03 (red star) using velocity models and phase subsets as indicated. For each event, the black circle indicates the Ground Truth location the white square at the end of the red line indicates the location for which the double-difference time residual norm was a minimum. Note that the scale of the maps is larger than in Fig. 4 in order to better display the mislocation vectors.
With this benchmark data set, unlike for the North Korean events, we have the exact coordinates of the sources. This may allow us to make direct measurements of the slowness vectors of the outgoing wavefield using so-called source-array analysis. Source-array analysis can be viewed as treating each event as a sensor observation at the coordinate of the event, where the data first are time-shifted to a common origin time. In this way, standard array processing methods can be applied. However, we have the disadvantage of not having exact origin times in the GT data and hence will have to solve for them. The initial origin time estimates we have from the correlation detectors are quite accurate, especially compared with say those obtained from classical event hypocentre estimation. However, they are based only on measurements on the ARCES array and will need to be adjusted for the short differences in the traveltime to ARCES based upon the locations within the explosion site.
The Bayesloc program (Myers et al. 2007) estimates joint probability distributions of hypocentre locations, origin times, arrival time uncertainty and corrections to modelled traveltimes for multiple events. If prior information is available, this can be applied to constrain the solution. In particular, we know the exact geographical hypocentres of these events and these can therefore be fixed in the input to Bayesloc. We have very accurate relative arrival time estimates from the cross-correlation calculations and the set of input files to run Bayesloc to solve for the event origin times is included in the Supporting Information. The Bayesloc solution uses traveltime tables as input and it could be argued that this introduces circularity, since the slowness vectors for the outgoing wavefield are implicit in the gradients of the traveltime tables. However, the ability of Bayesloc to compensate for deficiencies in traveltime predictions is exploited (e.g. Gibbons et al. 2017a) and we anticipate that this will mitigate any complications of circularity. To address the sensitivity of the origin time estimates to the traveltimes used, we ran the Bayesloc program using both the global AK135 velocity model and the regionally calibrated Fennoscandian model. The origin times solved for using the Fennoscandian model are given in Table 2. The AK135 model, with slightly slower crustal velocities, returned origin time estimates that were earlier than the estimates in Table 2 by an average of 0.0941 s (94.1 ms). However, the values of (t fescan orig − t AK135 orig ) were highly consistent from event to event with a standard deviation of 0.00145 s (1.45 ms). Since this variability is of the same order as the precision at which the traveltime differences are recorded, we are confident that the relative origin time estimates obtained are as accurate as we can demand given the accuracy of the differential traveltime data. The user is encouraged to investigate the stability of these origin time estimates by experimenting with parameters for the Bayesloc program and repeating the calculations as required.
Any error in the origin times given in Table 2 which is constant for each event is of no consequence to the source-array analysis. This is because the origin times are used only as a reference point for cutting data segments for each of the stations. Fig. 6 displays the waveforms for all 55 events on station SGF aligned according to these origin time estimates (panel a); five signals are singled out for closer inspection (panel b). The waveforms were upsampled to a sampling interval of 0.001 s in order to mitigate any inaccuracy due to the temporal discretization when the traces are cut according to the origin time estimates. The direction of the outgoing P-and S-wave fronts at SGF are estimated using broadband frequencywavenumber (f-k) analysis in the windows as indicated. Panels (c) and (d) show the results of f-k analysis in the two different timewindows together with the AK135-predicted slowness vectors. For both P-and S-wave arrivals, the f-k maximum falls slightly closer to Table 2. Origin time estimates for the Hukkakero explosion listed in Table 1 using Bayesloc using the Fennoscandian velocity model (Mykkeltveit & Ringdal 1981 the centre of the grid than the model-predicted slowness. This is to say that the outgoing wavefield has a slightly faster apparent velocity than is predicted by the velocity model-and it is the apparent velocity (the speed at which the wavefront appears to cross the array) which is critical for distance determination. This procedure was repeated for each of the six stations, for both P and S arrivals. Panel (a) of Fig. 7 displays the AK135-predicted slowness vectors for the outgoing P and S waves to each of the stations, while panel (b) shows the corresponding slowness vectors estimated empirically using the array analysis displayed in Fig. 6. The differences, analogous to the slowness corrections used to calibrate array measurements of incoming wavefronts (e.g. Schweitzer 2001), are displayed in panels (c) and (d) for P-and S-wave fronts, respectively. In panel (a), we see a clear difference between the slowness vectors for the the more distant stations (ARCES and KEV) and those for the closer stations (LP53, LP61, SGF and LP34). For ARCES and KEV, the predicted apparent velocities of P and S arrivals are about 8.0 km s -1 and 4.5 km s -1 , respectively. For the other stations, the predicted apparent velocities for P and S arrivals are about 5.9 and 3.5 km s -1 , respectively.  Table 2. The P-wave and S-wave time intervals are shaded blue and red, respectively. (b) Zoom of five signals from panel (a). f-k analysis performed on virtual source array, with each trace ascribed the Ground Truth coordinates of the corresponding explosion, for P window (c) and for S window (d). The small grey circles in panels (c) and (d) indicate the slowness vector with the optimal relative power and the blue and red circles indicate the theoretical P-and S-slowness vectors for the wavefield leaving the source under the AK135 velocity model and great-circle path assumption. Both measured P and S arrivals appear closer to the centre of the plots than the AK135 predictions indicating faster apparent velocity. models predict a Pg-type phase to this closest station, the apparent velocities are higher. The apparent velocities for station LP61, to the West, are significantly higher than predicted. However, we urge some caution regarding the accuracy of the empirical slowness vectors; some uncertainty will result from the relatively low resolution of the calculations. (The peaks in the relative power functions in Fig. 6 are relatively broad since our virtual source array has aperture 300 m, considerably smaller than most operational seismic arrays.) We recall also the results of Gibbons et al. (2010) in which significantly different slowness vectors for an incoming wavefield were obtained when using different frequency bands in the f-k analysis.
In several cases, a small azimuthal deviation is observed in the measured slowness vector. This was ignored in the subsequent relative event location calculation since the procedure described by Gibbons et al. (2017b) assumes great-circle propagation; a factor α multiplies the traveltime gradient. The α were calculated using the apparent velocities displayed in Fig. 7; α > 1 implies a slower-thanpredicted wave front while α < 1 indicates a faster-than-predicted wave front. . Slowness maps for the outgoing wavefield at Hukkakero using (a) the AK135 velocity model and great-circle path assumption and (b) the source-array analysis as displayed in Fig. 6. Panels (c) and (d) display the corrections for P and S arrivals, respectively, needed to move from the AK135-predicted slowness vector to the empirically determined slowness vector.
The relative event location results, with and without the empirical slowness corrections, are displayed in Fig. 8. The relative event locations using the corrections (displayed in panel b) have much smaller mislocation vectors than those without (displayed in panel a). The empirical slowness corrections have increased the apparent velocities of the outgoing wavefield to South, East, and West. The resulting location estimates are far closer to the GT locations than the uncalibrated estimates. Significantly, whereas the mislocation vectors displayed in panel (a) of Fig. 8 are relatively uniform within the clusters of explosions, the mislocation vectors in panel (b) are almost random with no clear pattern. We conclude that the event mislocations after the application of slowness corrections are primarily due to accuracy of time-delay measurements rather than the fundamental assumptions regarding the propagation of the outgoing seismic wavefield.
While we have obtained provisionally satisfactory relative location estimates for the 2007 Hukkakero explosions, we sit with the fundamental caveat that we used the GT coordinates to calculate the necessary slowness corrections. In other words, there is a circularity inherent in these calculations and we would not be able to follow an identical procedure for a set of earthquakes or explosions for which we had no a priori knowledge. What we have demonstrated is that systematic mislocation is most likely dominated by incorrect assumptions about the speed of the outgoing wavefield to each of the stations used. All approaches to accurate relative event location need to take this into account.

C O N C L U S I O N S
We present a set of GT explosions, well-recorded seismically up to distances of over 200 km, for which the source coordinates are known exactly. We provide a set of waveforms for all events and accurate time-delay measurements for all event and phase pairs. The distances between events are relatively short and are of relevance to interevent distances of small earthquakes and induced events. Assessing both the accuracy and precision of an algorithm's ability to locate the events presented here using the available data will be of great benefit to assessing the likely fidelity of maps of relocated seismicity in case studies of natural and induced seismicity. We demonstrate that, under the assumptions of two 1-D velocity models, the event-to-event distances are underestimated. From a source-array analysis, we obtain empirical slowness corrections which result in a significant improvement of the relative location estimates. In subsequent practical situations, we will not have enough information to perform the source-array analysis described here. Event mislocations are likely to be dominated by assumptions about the outgoing wavefield and we advocate theoretical advancements and implementations which address this issue. The GT data set includes 55 events, providing the possibility of working on many different subsets of events. We recommend that future studies address how much information about the relative location and the outgoing wavefield can be obtained in situations in which far fewer events and far less data is available. The GT will allow assessment and validation of any subsequent algorithms and software implementations. We hope that this study will motivate the publication of other GT data sets, with potentially very different network configurations and types of seismic source.

A C K N O W L E D G E M E N T S
We are grateful to Dr Stephen C. Myers and Dr Nima Nooshiri for thoughtful and thorough reviews of this manuscript which we believe have made significant improvements.
The ARCES seismic array is operated by NORSAR.
The station KEV (Kevo) is a collaboration between the Global Seismograph Nework (GSN)/ Incorporated Research Institutes for Seismology (IRIS), the United States Geological Survey (USGS), GEOFON, GEUS and the University of Helsinki. The waveforms were obtained via IRIS.
The station SGF is operated by Sodankylä Geophysical Observatory of the University of Oulu and the waveforms are available via GEOFON.
The POLENET/LAPNET project is a part of the International Polar Year 2007-2009 and a part of the POLENET consortium. The study was financed by the Academy of Finland (grant no. 122762) and the University of Oulu (Finland), the BEGDY programme of the Agence Nationale de la Recherche, Institute Paul Emil Victor (France) and the ILP (International Lithosphere Program) task force VIII, grant no. IAA300120709 of the Grant Agency of the Czech Academy of Sciences, , and by the Russian Academy of Sciences (programmes nos. 5 and 9).
We are grateful to Finnish Defence Forces for providing the coordinates of these events and for granting permission for them to be used in scientific research.
Most maps and figures in this paper are created using GMT software (Wessel & Smith 1995).