A Two-Station Seismic Method to Localize Glacier Calving

A method of determining glacier calving location using seismic wave arrival times from paired local seismic stations is presented. The difference in surface wave arrival times for each pair is used to define a locus (hyperbola) of possible origin. With multiple pairs, this can be used to triangulate for the origin of the seismic wave, which is interpreted as the calving location. This method is motivated by difficulties with traditional seismic location methods that fail due to the emergent nature of calving, which obscures the P and S-wave onsets, and the proximity of the seismometers, which combines body and surface 5 waves into one arrival. Human observed calving events are used to calibrate the seismic velocity for the method, which is then applied to other calving events from August 2014 to August 2015. From this, a catalogue of calving locations is generated, which shows that calving preferentially happens at the northern end of Helheim Glacier.


Introduction
The calving of marine-terminating grounded glaciers is a significant contributor to rising sea levels worldwide due to the massive volumes of ice involved that can suddenly be discharged into the sea.Depending on the exact location, the contribution from calving to sea level rise can be equal to the contribution from melt processes (Rignot et al., 2013;Depoorter et al., 2013).
However, the lack of understanding of the physical and mathematical principles that cause these events means that it is difficult to precisely forecast their contribution to expected sea level rise in the near future (e.g.Pfeffer et al. (2008), Meier et al. (2007)).
Calving glaciers can retreat in response to climate signals, and this can rapidly change the sea level (Meier and Post, 1987;Nick et al., 2013).A better understanding of calving processes is vital to developing accurate predictions of sea level rise.
The lack of understanding of why and how calving events happen makes it hard to create a general 'calving law' (Amundson and Truffer, 2010;Bassis, 2011).Recently, seismic arrays have been deployed to monitor glaciers and to detect calving (e.g.Walter et al. (2013); Amundson et al. (2012)).A common automated calving detection method is to take ratios of short-timeaverage and long-time-average seismic signals (STA/LTA).Large calving events can also generate glacial earthquakes, with surface waves detectable at a teleseismic level (Nettles et al., 2008;Nettles and Ekström, 2010;Tsai and Ekström, 2007).
Currently, most localization methods require visual confirmation of the calving location, unless they are sufficiently large to be seen by satellite imagery.There have not been enough direct observations of these smaller calving events (e.g.Qamar 1 The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -85, 2016 Manuscript under review for journal The Cryosphere Published: 20 April 2016 c Author(s) 2016.CC-BY 3.0 License.(1988), Amundson et al. (2008)) to form a general calving law.Calving events are somewhat intermittent, even if they also exhibit seasonality due to the seasonality of the mélange ice (Foga et al., 2014;Joughin et al., 2008), so monitoring equipment has to be deployed long-term in order to capture these events.Automatic methods like STA/LTA can help narrow down the manual search in satellite and camera imagery for calving, but ultimately, locating a calving event requires clear weather and well-lit conditions (O'Neel et al., 2007).An exception to this is radar, but radar cannot be deployed year-round without constant refueling and swapping out of data drives, and also has problems seeing through atmospheric precipitation.Recently, high frequency pressure meters, such as Sea-Bird Electronics tsunameters, have been deployed to monitor calving at Helheim (Vankova and Holland, in review).
In particular, land-based seismometers (providing seismic data) are much more useful than simple camera or satellite imagery because seismic arrays are not limited to daylight hours, are not affected by snow, can be deployed year-round without maintenance and also provide quantitative data to help estimate the magnitude of calving events.Seismic studies of calving have been done at the regional (<200 km) as well as the teleseismic level.Generally, teleseismic detections of calving are done via low frequency surface waves (e.g.Walter et al. (2012); O'Neel and Pfeffer (2007); Chen et al. ( 2011)), while local detections are done at some range of frequencies within 1-10 Hz (e.g.(Bartholomaus et al., 2012;Amundson et al., 2008Amundson et al., , 2012;;O'Neel et al., 2007)).
Until recently, seismic signals generated by glacial calving were believed to be caused either by capsizing icebergs striking the fjord bottom (Amundson et al., 2012;Tsai et al., 2008), or by sliding glaciers that speed up after calving (Tsai et al., 2008).However, Murray et al. (2015) found that glacial earthquakes at Helheim Glacier are caused by glaciers temporarily moving backward and downward during a large calving event.It is not yet known how to fully categorize and characterize different calving events.For example, Nettles and Ekström (2010) found that only capsizing icebergs generate observable low-frequency surface wave energy, with calving events creating tabular, wide icebergs not generating glacial earthquakes.
Calving seismic signals typically have emergent onsets with dominating frequencies around the order of 1-10 Hz (Richardson et al., 2010;O'Neel et al., 2007;Amundson et al., 2012).The emergent signals make it hard to accurately identify a P-wave onset time, let alone a S-wave onset time, which prevents the traditional seismic triangulation method that takes the time difference between the P-and S-waves to generate a distance to the epicenter (Spence, 1980).The other main method involves calculating backazimuths from a ratio of easting and northing amplitudes of P-waves from a broadband seismic station (e.g.Jurkevics (1988)); this also fails due to the proximity of the station and the high speed of the sound waves (4 km/s through ice) which makes the waves arrive near-simultaneously.Although an emergent signal does not prevent the detection of calving as the STA/LTA method mentioned above shows, it does thwart the localization of the calving as the different seismic waves cannot be separated.
A more recent method for localizing calving events is the use of frequency dispersion, which uses a regional array (100-200 km away) of hydroacoustic stations to estimate a range between event and detector and combines this with an azimuth to create a unique intersection (Li and Gavrilov, 2008), as the stations are sufficiently far to distinguish P-and S-waves.This method has similar precision to using intersecting azimuths from two remote stations, which is enough to identify which glacier the calving occurred at, but not enough to localize the event within the glacier.This motivates the creation of a more precise localization method that does not require the separation of different seismic waves.We denote this method the "Hyperbolic Method", as it uses the difference between the signal onset times at two nearby stations to generate a unique hyperbola, on which the event is localized.Calving events were observed with the naked eye on August 12 and August 13 2014 at Helheim Glacier (in East Greenland) and on June 10 at Jakobshavn Glacier (in West Greenland).These events are used to calibrate the local surface wave velocity on the glacier surface.The Hyperbolic Method is then applied to calving events at Helheim to localize the epicenters of the seismic signals generated during calving events.The signals for our observed calving events match those of Amundson et al. (2012) very well in both frequency distribution and shape, with an emergent onset and a power spectrum dominating in the 1-10 Hz range (see Figure 2).In contrast, teleseismic events from regional earthquakes have much lower frequency signals, below 0.1 Hz.In Figure 3, we show a M = 7.0 earthquake that happened in Atka, Alaska, USA on August 30 2013 1 where the dominant frequencies received at the HEL seismometers are all well below 1 Hz.This means we can easily separate calving events from regional seismic activity by using a bandpass filter (Butterworth, two-pole and zero-phased).We bandpass filter between 2-18 Hz based off these spectrograms in order to maximize the signal-to-noise ratio for our hyperbolic location method.From this we are able to create a catalogue of calving events to run our Hyperbolic Method algorithm.Calving events, with the exception of events in January/February 2015 for 5 which imagery is too snow-covered to use, are confirmed with local camera imagery and satellite imagery from the Rapid Ice Sheet Change Observatory (RISCO) website 2 .

Localization Methods and Results
After isolating the calving events, we now present the Hyperbolic Method and apply it to generate a catalogue of calving locations.

Hyperbolic Method
This method relies on the fact that a hyperbola can be geometrically defined as the locus (set of points) with a constant path difference relative to two foci, as seen in Figure 4.In our case, each pair of seismometers act as foci.
Assuming that the speed of seismic waves across Helheim does not vary horizontally, a calving event that happened exactly at the midpoint of the two seismometers (or indeed, any point along the perpendicular bisector of the two seismometers) would register simultaneously at the two seismometers.Similarly, if the event happened closer to HEL1, the seismic waves would arrive slightly earlier to HEL1, and the locus of possible calving locations would instead be the set of all points whose distance from HEL1 is shorter than HEL2 by a fixed length, 2a as seen in Fig 4, which would be the product of the speed of the waves through the glacier (v seismic ) and the time lag in signal arrival (∆t).This difference of 2a is defined for a hyperbola with equation x 2 /a 2 − y 2 /b 2 = 1: see Figure 4. We may use the time lag of the signal arrivals at the two seismometers (which become the foci of the hyperbola) to determine the path difference of the signals, and from this deduce the locus of possible signal sources.
One of the arms (in Figure 4, either the left curve or the right curve) of the hyperbola may always be eliminated, as we can always see which seismometer the event was closer to by seeing which station has the first signal onset.This remaining arm will intersect uniquely with the calving front, which will give the location of the calving.If the calving front is not known, using multiple pairings of the four stations, the calving location can be estimated using triangulation.
One key unknown in this method is the speed of the seismic waves through the glacier.We do not attempt to identify whether the seismic waves are P, S or surface waves, except to note that it appears that there is a mixture of all of these.We do assume that the seismic wave travels at the same speed from the epicenter (the vertical projection of the true seismic source to the surface) to each station.The dependence of wave speed on glacier depth is not important for this method as long as the effective lateral speed to each seismometer is the same in each direction via symmetry from the source.We also assume that the glacier surface, calving epicenter and seismometers are all coplanar, so that the hyperbolas can be kept two-dimensional for simplicity.In reality, there is some elevation between the seismometers and the glacier surface, though this distance (<300 m) is so much shorter than the seismometer separation (>6500 m) so that refraction at the ice/rock boundary is likely negligible for characterizing the hyperbola.However, this method would become more powerful with three-dimensional hyperboloids instead of two-dimensional hyperbolas.
This method therefore requires evaluating the time lag between the signal arrival times at each seismometer (see Figure 6), and obtaining the speed of the seismic waves through the glacier.As we are not distinguishing whether the seismic waves are body or surface ones, we rename the variable as v eff , which is the effective speed of the seismic packet through the glacier surface using the above assumptions.
To identify the time lag, we first try using cross-correlation of the signals.However, as seen in Figure 6, the signals are sometimes too different to generate an accurate lag time.For the top two subplots, the cross-correlation gives 0.47 s, which is a plausible value by eye, but for the panels HEL2 and HEL3, cross-correlation gives a value of 2.4 s, which is not plausible by eye.Instead of using cross-correlation, we use an automated script that searches through the signal for the first instance of a slope exceeding 44% of the standard deviation of all slopes at for a time step of 0.025 s (corresponding to the sampling frequency of 40 Hz).This value of 44% was empirically determined as this produced the closest match to cross-correlation for signals that were qualitatively similar.
To calibrate the seismic velocity v eff , we use two events that were observed locally at Helheim in August 2014 (with only two seismometers) as well as one event at Jakobshavn Glacier in West Greenland (where there is a similar setup with three seismomters) in June 2015.In both cases, the approximate seismic velocity corresponding to the observed location of calving was 1.17 ± 0.1 km/s (see Figure 5).For all further plots, we therefore use v eff = 1.17 km/s.
Once we generate four hyperbolas, we may take their intersection to be an estimate of the calving area, as seen in Figure 7.
Applying this to our entire catalogue of calving events (excluding four out of fifteen due to unclear and noisy signal onsets), we have Figure 8.

Discussion
The Hyperbolic Method described in this paper offers a powerful alternative to traditional seismic location techniques, which are more suited for regional seismic arrays that can distinguish between the different seismic wave types (e.g.O'Neel et al. (2007)).Moreover, these distant arrays do not give the kind of precision that local arrays would have, as small errors on a regional azimuth translate to a large area of uncertainty on the local glacier surface.The Hyperbolic Method takes advantage of the stations' proximity to calving events and does not require separating out the different wave phases, thus solving the P-wave identification problem that hampered localization techniques from Amundson et al. (2008) and Richardson et al. (2010).The method also offers advantages over traditional calving detection methods, which require the use of a local camera and/or satellite data to visually confirm that calving took place.As seen in (Amundson et al., 2012), calving generates a characteristic seismic signal, which we also see in Figures 2 and 3, where we directly compare teleseism with calving and note that teleseism does not have any energy above 5 Hz or so and have most of their energy <1 Hz.This is likely because higher frequency signals are severely attenuated by the time they reach the seismometers; in any case, this means that seismometers can be used to monitor glaciers and quickly identify calving when power in the 2-18 Hz range exceeds some ratio above the ambient noise.Importantly, this monitoring could take place year-round during the night and also cloudy days, thus replacing satellite imagery, camera imagery and radar monitoring as the primary method for locating calving.
The seismic wave detected here does not appear to be a P-wave.Its speed of 1.17 km/s is much lower than would be expected for a P-wave in ice, which is typically larger than 3 km/s (e.g.3.25 km/s in O' Neel et al. (2007) or 3.8 km/s in Robinson (1968)).
Moreover, the particle plots in Figure 9 clearly show the characteristic elliptical shape of a Rayleigh wave (a surface wave).
The Rayleigh waves, which are in theory parallel to the vertical axis, appear slanted in Figure 9.It is possible that the mix of different wave phases (e.g.Love waves, also a surface wave) has sheared the Rayleigh wave such that it is no longer parallel to the vertical axis.In any case, the lack of a linear polarization (as would be expected for a P-wave) is clear.This means that the wave packet is likely dominated by surface waves.Given our characteristic velocity on the order of 1 km/s with frequencies The Cryosphere Discuss., doi:10.5194/tc-2016-85,2016 Manuscript under review for journal The Cryosphere Published: 20 April 2016 c Author(s) 2016.CC-BY 3.0 License. Figure 6.Seismic signals for a calving event at Helheim Glacier on September 2, 2014.The signal onset times are determined using an automated script that searches for the first instance of a gradient exceeding a particular threshold as defined in the text.The differences in the wave onset times is then used to generate a characteristic path difference for each hyperbola.
of the order of 10 Hz (see Figure 2), this corresponds to a surface wavelength of order 100 m, which is small enough to be affected by crevasses along the surface of the glacier which are of depths of order 50 m.This means that we can reasonably expect these air-filled crevasses to affect the surface wave velocity, making our reported speed of 1.17 km/s a plausible value.
Because we are only working with surface waves, this limits our localization technique to just the epicenter of a calving event, with no suggestion of a focal depth.Moreover, we have assumed a planar ice front for simplicity.We cannot, at this point, use this method to locate any events happening below the glacier surface.It is possible that this method could be extended to determine the depth at which calving occurs by using a 3-D hyperboloid instead of 2-D hyperbolas.However, the main limitation of this is that the density of glaciers changes with depth and so the seismic wave speed itself should be a function of depth, and so the locus of possible origins would not form a hyperbola.
The calculation method we have used ignores the presence of the rock, as the proximity of the seismometers to the glacier 10 makes the time the wave spends travelling across rock negligible compared to the time spent on the glacier.The localization of edge events, for which the rock surface contributes a greater proportion of the path of the surface wave, are not significantly affected by including the rock medium into the surface velocity calculation.However, our method does not take into account the refraction at the ice-rock interface.Again, due to the ice dominating the wave path from the source to the seismometers,  we may assume that the refraction has a negligible affect on the azimuthal measurement.The rock also gives elevation to the seismometers, so that a linearly polarized P-wave would also have some vertical component.However, we do not see linear polarization in our particle plots in Figure 9.
Calving appears to preferentially happen at the northern end of the calving front at Helheim.This is consistent with the topography of the bedrock at Helheim as seen in Figure 10, where the northern half is on the order of ∼200 m deeper than the southern half (Leuschen and Allen, 2013).It is possible that the deeper the ice, the higher the freeboard of the ice front and the greater the stresses that affect the calving front.Camera imagery further suggests a structural difference between the northern and southern halves of Helheim, as seen in Figure 11, where there are wider ridges in the north (left side of image) as compared to the south.This may also mean that the surface velocities are different in each half, which would affect the localization results.The topographic differences of both the ice surface and ice bottom may contribute to why we see calving primarily in the northern half of Helheim.
It is possible to constrain the fault size of the rupture caused by calving.Using a shear model from Brune (1970), the radius r 0 of a circular fault is inversely proportional to the corner frequency f c of a S-wave is given by where β is the shear velocity, K c is a constant, equal to 2.34 for Brune's source model (Gibowicz and Kijko, 2013).From Figure 12, the corner frequency is approximately bounded between 2.5 and 10 Hz.Using a shear velocity of 1.3 km/s (given that the Rayleigh wave velocity is typically 92% of the S-wave velocity as per Sheriff and Geldart (1995)), this bounds the     magnitude.100 m is more on the order of a crevasse, which also occur during/before events, so it is possible that crevassing events continue to happen during the calving event and obscure the power spectrum seen in Figure 12.It is also possible that calving events are a sequence of small cracks on the order of 100 m that quickly cascade en masse into a larger event.A fracture size of order 1 km would require a corner frequency of order 0.1 -1 Hz, which we do not observe.

5
Our results show that there are ways to get around the emergent P-wave problem for glacial calving, which prevents traditional seismic location methods, via the development of a hyperbolic method that simply measures the time delay in the signal arrival times at two seismometers.This model can be made more complex by using hyperboloids of revolution in place of the hyperbolas, which would also give insight into the depth of the calving location.With three seismometers, triangulation becomes possible, and calving events can be detected and located even without satellite or camera imagery.This method can 10 be automated for spectra with good signal-to-noise spectra, but signal onset detection of noisy spectra still requires manual  inspection at this point.Our catalogue of calving events at Helheim suggests that calving typically initiates at the northern half of the calving front, which will help to constrain model simulations of glacier dynamics at Helheim.
seismometers (HEL1: Nanometrics Trillium 120, HEL2-HEL4: Nanometrics Trillium 240) with dual sampling rates of 40 Hz and 200 Hz were deployed around the mouth of Helheim Glacier, as seen in Figure 1.HEL1 and HEL2 were deployed in August 2013, while HEL3 and HEL4 were deployed in August 2014.They were synchronized with Coordinated Universal Time.These stations picked up seismic activity from both calving as well as distant earthquakes, so we first inspect the frequency distributions of the signals to isolate calving events.

Figure 2 .
Figure 2. Spectrograms for CE-I, the calving event of August 12 (left), and of CE-II on August 13 (right) in 2014.The y-axis of the top panels shows counts, a dimensionless quantity showing relative ground motion of the instrument measured within the digitizing system.The easting amplitude of the seismometers is depicted here.The seismogram (top) and spectrogram (bottom) of each event share the same time axis for direct comparison.

Figure 3 .
Figure 3. Spectrograms for the teleseismic event from Atka, AK, United States discussed in from HEL1 (L) and HEL2 (R), showing the different frequency distribution when compared to the calving event in the previous figure.
Figure 5.A calving event (yellow diamond) at Jakobshavn Glacier in West Greenland as shown by Google Earth.The black circle indicates the uncertainly of the observed calving location.The same method as at Helheim is used to generate a time lag of signal onset, which is combined with a range of seismic velocities to determine which hyperbola pair intersects closest to the observed location.Note that the satellite imagery here does not show an up-to-date depiction of the calving front, and so the calving event falsely appears to originate in the ocean.

Figure 7 .
Figure7.Hyperbolas corresponding to four pairings of seismometers with signal lags as seen in Figure6.Of the six possible pairings of four seismometers, we use HEL1-HEL2, HEL1-HEL3, HEL2-HEL4 and HEL3-HEL4 as these have the greatest distance between the stations to improve precision and also avoid traversing large areas of rock as in the HEL2-HEL3 and HEL1-HEl4 pairings.The overlapping area, shown here in red, is assumed to be the location of the calving event.The centroid of the overlapping area is indicated with a red 'x'.

Figure 8 .
Figure 8. Catalogue of all calving events with clear signal onsets at Helheim Glacier from August 2014-August 2015 overlaid on Google Earth imagery of Helheim Glacier.Each color corresponds to a calving event, with only the area of overlap as shown in Figure 7 depicted.Note that the calving front is shown only for July 2014.

Figure 9 .
Figure 9. Particle plots of surface wave arrivals for the calving event of July 7 2015, split into radial and transverse components.The characteristic elliptic shape of the Rayleigh wave is clearly visible in the radial component of the particle plot.

Figure 10 .
Figure10.The calving events from Figure8shown with the bedrock topography from NSIDC IceBridge MCoRDS(Leuschen and Allen, 2013) overlain.From depth soundings made in the mélange we know that the depth is approximately 600 m, and from the MCoRDS overflight of the ice front we know that some of the ice upstream of the ice front is deeper than 600 m, and so we only show depth values that are near or below 600 m.

Figure 11 .
Figure 11.Helheim Glacier viewed from the air; the left side of the image is north and the right side of the image is south.Image provided by the Geological Survey of Denmark and Greenland.

Figure 12 .
Figure12.A typical power spectrum for the two observed calving event in August 2014, for the time period indicated between the two red lines.The shaded inset in each subplot shows the zoomed-in view of the section of the seismic trace that was used to make the power spectrum.

Acknowledgements.
The fieldwork necessary to collect this seismic data was made possible by the Center for Global Sea Level Change, grant G1204 of the NYU Abu Dhabi Research Institute and the Undergraduate Research Fund at NYU Abu Dhabi.Denise Holland provided the field logistics support in Greenland.5 13 The Cryosphere Discuss., doi:10.5194/tc-2016-85,2016 Manuscript under review for journal The Cryosphere Published: 20 April 2016 c Author(s) 2016.CC-BY 3.0 License.