Introduction
The relationship and potential feedback between changes in surface melt due to global climate change and ice velocity on the Greenland ice sheet (GrIS) may play an important role in future ice-sheet mass balance (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Rignot, Box, Burgess and HannaRignot and others, 2008). A clear link has been identified between ice velocity change and surface meltwater availability (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference DasDas and others, 2008; Reference JoughinJoughin and others, 2008; Van de Reference Van de WalWal and others, 2008; Reference Palmer, Shepherd, Nienow and JoughinPalmer and others, 2011). However, the efficiency of the subglacial drainage system determines how meltwater supply to the glacier bed enhances basal motion (Van de Reference Van de WalWal and others, 2008; Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010; Reference Sundal, Shepherd, Nienow, Hanna, Palmer and HuybrechtsSundal and others, 2011; Reference SoleSole and others, 2013). Consequently, characterizing subglacial and englacial water flow is necessary to predict the effect of increased surface melt on ice-sheet dynamics. Traditional glaciological techniques, such as tracer experiments (e.g. Reference Werder, Schuler and FunkWerder and others, 2010; Reference ChandlerChandler and others, 2013), borehole geophysics (e.g. Reference Iken, Echelmeyer, Harrison and FunkIken and others, 1993) and ground-penetrating radar (Reference Fountain, Jacobel, Schlichting and JanssonFountain and others, 2005; Reference Catania, Neumann and PriceCatania and others, 2008; Reference Catania and NeumannCatania and Neumann, 2010) can provide valuable information on englacial and subglacial water passages. However, the different techniques are mostly limited in their spatial and temporal resolution.
As an alternative, seismic monitoring allows for detection, location and characterization of seismic sources related to dynamic and hydraulic processes throughout the ice. Seismic techniques can target a glacier region, limited only by the seismic network aperture, seismic signal attenuation within the ice and seismic background noise. In previous studies, seismic measurements have already been used to monitor glacier dynamics in both non-polar (e.g. Deich-Reference Deichmann, Ansorge, Scherbaum, Aschwanden, Bernardi and Gudmundssonmann and others, 2000; Reference Walter, Deichmann and FunkWalter and others, 2008; Reference Roux, Walter, Riesen, Sugiyama and FunkRoux and others, 2010; Dalban Reference Dalban Canassy, Faillettaz, Walter and HussCanassy and others, 2012) and polar regions (e.g. Reference Winberry, Anandakrishnan and AlleyWinberry and others, 2009; Reference WalterWalter and others, 2012; Reference Zoet, Anandakrishnan, Alley, Nyblade and WiensZoet and others, 2012). Seismic signals can be interpreted in terms of source processes (Reference Walter, Clinton, Deichmann, Dreger, Minson and FunkWalter and others, 2009; Reference West, Larsen, Truffer, O’Neel and BlancWest and others, 2010) and reveal insights into subsurface glacier hydraulics. This approach is motivated by the analogy between fluid-related seismic sources in volcanic environments (magma and water) and seismic signals detected on glaciers related to water (St Reference Lawrence and QamarLawrence and Qamar, 1979; Reference West, Larsen, Truffer, O’Neel and BlancWest and others, 2010). Findings from volcano seismology and classical seismology techniques can therefore be applied for the analysis and interpretation of seismic signals on glaciers (e.g. Dalban Reference Dalban Canassy, Walter, Husen, Maurer, Faillettaz and FarinottiCanassy and others, 2013; Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others, 2013).
Using data from a temporary high-density seismic network in western Greenland’s ablation zone, we have identified seismic signals related to the presence or flow of water within the ice sheet. Besides typical surface crevasse icequakes, we were able to identify and locate icequakes at >100 m depth (deep icequakes). However, the most prominent signal is a reoccurring, sustained seismic tremor lasting several hours. Its epicentre can be reliably located to the surface entrance of an englacial channel (‘moulin’). The frequency signature of this ‘moulin tremor’ consists of discrete frequency bands, whose temporal variations correlate closely with moulin water level. Our results confirm that seismogenic processes typical for Alpine environments also occur on the GrIS. In addition, the finding of the moulin tremor signal documents a previously unreported signal. This suggests that seismic monitoring may provide new insights into the variability and evolution of the ice sheet’s subsurface drainage system in the ablation zone of the GrIS.
Study Site
Our seismic network was installed on the flowline in the Sermeq Avannarleq ablation zone (Fig. 1), ~25 km downstream of Swiss Camp (e.g. Reference Steffen and BoxSteffen and Box, 2001). This relatively small and slow-flowing region is located 30 km north of Jakobshavn Isbræ, one of Greenland’s fastest-flowing glaciers (12.6 km a–1; Reference JoughinJoughin and others, 2008) which drains ~7% of the entire GrIS (Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991). In our study area, the GrIS has a thickness of ˜620 m (determined from deep drilling as described below). Using ablation stakes, we measured a surface melt of ˜2 m w.e. for July 2011, with total ablation of 6 m w.e. a–1. The ice flows with a surface velocity of ˜100 m a–1, exhibiting meltwater-dependent diurnal and seasonal fluctuations (Reference Hoffman, Catania, Neumann, Andrews and RumrillHoffman and others, 2011; Reference McGrath, Colgan, Steffen, Lauffenburger and BalogMcGrath and others, 2011).
(a) WorldView-2 image (red band, from 20 June 2011, Polar Geospatial Center) of study area and instrumentation. Seismic stations (triangles), epicentres of the two deep icequakes discussed in the text (red stars), major moulins (red crosses), surface stream gauge (yellow cross) and location of hot-water drilling (green circles) are shown. The core seismic network (blue triangles) was arranged inside a circle of 800 m diameter and supplemented with additional stations aimed at enlarging the network’s aperture. (b) Location of study area (red circle) on the west coast of Greenland.
The seismic network is located in an area with several supraglacial streams and a few small lakes. The largest surface stream in our seismic network (average discharge ˜2.5 m3 s–1) fed into a major moulin (M1 in Fig. 1) with an extension at the surface of at least 5 m in width and 10 m in length. Another moulin (M2) is about one magnitude smaller in water discharge than moulin M1 and located ˜200 m south of M1. Surface crevassing inside and near the network was inhomogeneous, varying between hair fissures and crevasses with surface widths up to 2 m. At the study site, the ice surface is completely exposed during July and August. The ice temperature at the study site was measured with temperature sensors in boreholes. The ice column consists of ˜600 m of cold ice overlying 20 m of temperate ice at the glacier bed.
Experimental Set-Up
The seismic network contained 17 seismometers recording between 2 July and 17 August 2011 (Fig. 1). Nine 1 Hz Lennartz LE-3D surface seismometers arranged around the perimeter of an 800 m diameter circle formed the network core. Additionally, three 1 Hz Lennartz LE-3D/BH borehole seismometers were installed between 2 and 3 m depth to enlarge the network aperture. We also deployed three 8 Hz borehole seismometers (GS11-D) in boreholes at depths of 150 and 350 m. For testing purposes we installed two co-located broadband seismometers (Trillium Compact and Trillium Compact All Terrain, both 120 s corner period) inside the network. The seismometers were each equipped with a Taurus digitizer from Nanometrics, continuously recording at 500 Hz sampling frequency.
In order to reduce effects of meltout and exposure to wind and rain, the surface seismometers were installed in holes ˜0.3 m below the ice surface and covered with a fleece tarp for insulation. The installation set-up was similar to earlier deployments in the Swiss Alps (Reference Deichmann, Ansorge, Scherbaum, Aschwanden, Bernardi and GudmundssonDeichmann and others, 2000; Reference Walter, Deichmann and FunkWalter and others, 2008). Due to the high melt rate, the seismometers were re-levelled every 1–2 days (depending on weather conditions). This levelling work and the occasional exchanges of CompactFlash® cards for data backup made data gaps unavoidable. Station visits to the three LE-3D/BH borehole seismometers were less frequent (every 4–5 days) as they were installed below the ice surface. Nevertheless, high melt rates required reinstallation of these instruments on several occasions.
The seismic network was installed around one of the ROGUE (Real-time Observations of Greenland’s Under-ice Environment) drill sites. One of the major ROGUE project field activities was hot-water drilling to the glacier bed. Additional borehole geophysical and glaciological measurements of the ROGUE project were made available to the present study. These include measurements of borehole water pressures, surface velocities, surface melt and stream discharge into moulin M1 (Fig. 1). In addition, we used data from a water pressure sensor that had been lowered into the main moulin (M1). Unfortunately, the absolute height of the pressure sensor is subject to considerable uncertainty, and the moulin water level frequently fell below the sensor, resulting in data gaps. Despite these shortcomings, the sensor reliably measured moulin water-level fluctuations during high surface stream discharge.
Seismological Observations
Our seismic data document a large variety of glacier-related seismic signals, some of which had been previously observed on non-polar glaciers (e.g. Reference Deichmann, Ansorge, Scherbaum, Aschwanden, Bernardi and GudmundssonDeichmann and others, 2000; Reference Roux, Marsan, Metaxian, O’Brien and MoreauRoux and others, 2008; Reference Walter, Deichmann and FunkWalter and others, 2008; Reference West, Larsen, Truffer, O’Neel and BlancWest and others, 2010; Reference Mikesell, Van Wijk, Haney, Bradford, Marshall and HarperMikesell and others, 2012). Here we focus on those events, which can be related to (melt)water flow within the glacier. We distinguish between short-duration icequakes 0.1 s to several seconds in length (Fig. 2a and b) and longer-lasting tremor-like signals lasting minutes to several hours (Fig. 2c and d). The seismic events can also be distinguished by frequency content. Whereas icequakes show a wide frequency spectrum in the 10–200 Hz range (Fig. 2a and b), tremor-like signals are narrow-banded in the 1–10 Hz range (Fig. 2c and d).
Different types of event, with corresponding spectral content calculated with the Fourier transform in the right-hand column. Note the different timescale for the waveform (left column) and the resulting differences in total amount of energy in the power spectrum density (PSD, right column). (a) Surface crevasse icequake with P-arrival and dominant Rayleigh wave; (b1) deep icequake with dominant P-arrival, higher frequencies; (b2) deep icequake with resonance coda; (c1, c2) short-duration tremors; and (d) long-duration tremor.
Detection of short-duration signals
Short-duration icequakes occur up to ten times per minute (Fig. 3). We apply a classic short-term/long-term average (STA/LTA) trigger algorithm to detect short-duration ice-quakes in our continuous data (e.g. Reference AllenAllen, 1978; Reference WithersWithers and others, 1998). This algorithm divides a subsection of a seismic trace into a short-term time window followed by a long-term time window. Within both windows, the average of the squared amplitudes is calculated. When the ratio of the two exceeds a user-defined threshold, the trigger flag is turned on until the ratio again falls below this value. Following Reference Dalban Canassy, Faillettaz, Walter and HussDalban Canassy and others (2012), we choose 1 s and 10 s for the short-term and long-term average, respectively. We require a threshold of 3 for the STA/LTA ratio and simultaneous triggering on at least four stations for the declaration of a detection. Note that the STA/LTA trigger algorithm was applied to the vertical component seismo-grams only.
Sample waveform (vertical component), with red rectangles marking events discussed in Figure 2 (analogue a is surface crevasse icequake; b is deep icequake with coda; c is narrowband tremor).
On average, 5400 2100 icequakes per day were detected. The absolute number of triggered icequakes is highly dependent on the used threshold of the STA/LTA rate as well as the required number of triggering stations. Following Reference Walter, Deichmann and FunkWalter and others (2008), we investigate potential changes in trigger sensitivity in Figure 4. In this figure, each dot indicates one icequake detection. On the y –axis we plot the median of maximum seismogram amplitude at all recording stations (=signal strength in counts) against the trigger time on the x –axis. Different point colours mark different days of icequake detection (in Coordinated Universal Time (UTC)). Local time is UTC – 2 h, and time difference between sunrise and sunset varied between 0 hours (permanent sunshine) at the beginning of the seismic campaign (2 July 2011) and 17 hours at the end of the seismic campaign (17 August 2011).
Icequake signal strength as a function of time. Each point marks one icequake. The colours correspond to different days in UTC. Note the diurnal variability in maximum and minimum detected signal strength, also visible in zoomed window of Figure 13a.
Figure 4 shows that the maximum and minimum signal strength changes within 1 day. In particular, during the early morning hours, relatively weak events are detected, which are no longer detected during afternoon/night times. On Alpine glaciers, Reference Walter, Deichmann and FunkWalter and others (2008) and Reference Dalban Canassy, Faillettaz, Walter and HussDalban Canassy and others (2012) linked such changes in icequake detections to melt-induced seismic background noise on the glacier: as meltwater availability increases towards the afternoon, seismic background noise increases and the sensitivity of the STA/LTA trigger algorithm decreases.
With diurnal changes in trigger sensitivity, any temporal statistics of icequake detection must be interpreted with care. Figure 5a shows the daily water heights (as well as their mean) in the major surface stream feeding the main moulin, M1 (Fig. 1), which can be taken as a proxy for the availability of surface melt. For the same time window, Figure 5b shows the hourly icequake detections stacked from 4 July to 15 August 2011 visible with the blue bars in the histogram. The icequake detection (blue bars) reaches a global maximum around 7:00 UTC, when stream water height is near its daily minimum, with another local maximum near 19:00 UTC. If we consider only icequakes exceeding signal strength of 600 counts (green bars), the maximum in the early morning hours dissolves completely. Increasing the signal threshold to >800 counts (red bars) has little effect on the overall appearance of the histogram, and icequake rejections are evenly distributed throughout the day. This strongly suggests that the global detection peak in the morning hours is due to weak signals that do not trigger detection in the afternoon hours, caused by increases in background noise. By eliminating the effect of trigger sensitivity (red and blue bars, Fig. 5b), a clear diurnal variance in icequake detection correlated with stream water height is visible.
(a) Water level in stream feeding moulin M1. Red line marks mean water level of individual days (grey lines). (b) Number of triggered events (stacked) per hour of day (UTC) between 4 and 15 August 2011 (blue). A diurnal trend with two local maxima (morning and afternoon) is visible. When removing events with an amplitude lower than 600 counts (green) or 800 counts (red), a diurnal systematic with higher activity in the afternoon/night is visible.
Surface crevasse icequakes
The sample waveform in Figure 3 includes three different types of short-duration icequakes. Icequake type a (also Fig. 2a) is characterized by a dominant Rayleigh wave and few or no visible P- and S-phase arrivals (e.g. Reference Lay and WallaceLay and Wallace, 1995). The frequency spectrum of this event type peaks between 10 and 50 Hz. These characteristics are typical for near-surface events associated with surface crevasse activity as previously investigated on Alpine glaciers (Reference Deichmann, Ansorge, Scherbaum, Aschwanden, Bernardi and GudmundssonDeichmann and others, 2000; Reference Walter, Clinton, Deichmann, Dreger, Minson and FunkWalter and others, 2009). Accordingly, we refer to them as ‘surface crevasse icequakes’. Visual inspection of our continuous seismic record suggests that, similar to Alpine glaciers, these events constitute the vast majority of our short-duration seismicity. Therefore, the diurnal statistic of icequake detection reflects diurnal fluctuations in surface crevasse activity most likely depending on meltwater availability (Fig. 5).
Deep icequakes
Icequake type b in Figure 3 and shown in Figure 2b1 and b2 is characterized by impulsive first motions and dominant Sand P-phases and little or no Rayleigh wave. Its spectral content is shifted toward higher frequencies compared to surface crevasse icequakes, with peaks often beyond 100 Hz. These features resemble Alpine icequakes, which locate at depths below the crevasse zone (Reference Deichmann, Ansorge, Scherbaum, Aschwanden, Bernardi and GudmundssonDeichmann and others, 2000; Reference Walter, Clinton, Deichmann, Dreger, Minson and FunkWalter and others, 2009), and we refer to these events as ‘deep icequakes’. So far, we have identified two deep icequakes via visual inspection. Our large icequake volume demands automatic techniques for a systematic search for deep icequakes (e.g. Reference Hammer, Ohrnberger and a¨hHammer and others, 2012). This is the subject of ongoing investigation.
Our two detected deep icequakes exhibit an important difference: the waveform of the event shown in Figure 2b2 is followed by a narrow-banded coda of ˜4 Hz (Fig. 6) whereas it is missing for the second icequake shown in Figure 2b1. Impulsive high-frequency seismic events followed by a low-frequency coda often occur in volcanic environments and these are known as hybrid events (e.g. Reference ChouetChouet, 1996). The high-frequency content is likely related to brittle failure and purely elastic processes, whereas the low-frequency coda is related to resonances in liquids or gases of magmatic or geothermal origin. Reference West, Larsen, Truffer, O’Neel and BlancWest and others (2010) have proposed an analogy between glacier-related hybrid events and equivalent sources on volcanoes. They explain hybrid events as water-driven fracturing, so-called hydrofracturing, followed by resonances caused by rush of water into the newly opened space.
Spectrogram (lower) of deep event with its monochromatic coda (red rectangle, b) with the corresponding waveform (upper) filtered between 1 and 80 Hz. In addition, three surface crevasse related events are visible in the waveform and as high-energy spots (a) in the spectrogram.
Location of deep icequakes
We use the location software package NonLinLoc (Reference Lomax, Virieux, Volant and ThierryLomax and others, 2000) to locate the two deep icequakes. The impulsive onsets of P- and S-wave arrivals provide high-quality arrival-time measurements, which can be inverted for origin time and source location using a nonlinear, probabilistic approach. We assume a simple homogeneous one-dimensional velocity model with a P-wave velocity of 3.6 km s–1 and S-wave velocity of 1.8 km s–1 after Reference Walter, Deichmann and FunkWalter and others (2008) and Reference Rial, Tang and SteffenRial and others (2009). Each seismic phase arrival time estimate is associated with an individual uncertainty between 0.002 and 0.008 s accounting for different waveform quality. The arrival time of S-waves is subject to higher uncertainties than for P-waves. Furthermore, at some stations S-wave arrivals are not clearly visible.
With the arrival times and associated uncertainties, NonLinLoc computes the posterior probability density function (PDF) of the earthquake location problem, which can be irregular in shape and show multiple minima (Reference Lomax, Virieux, Volant and ThierryLomax and others, 2000). Reference Dalban Canassy, Walter, Husen, Maurer, Faillettaz and FarinottiDalban Canassy and others (2013) recently applied this technique successfully to icequakes recorded on Triftgletscher, Swiss Alps.
Location results of NonLinLoc are shown by the maximum likelihood hypocentre location, the scatter density cloud and a traditional 68% confidence ellipsoid (Fig. 7). The latter two describe location uncertainties as given by the network geometry and uncertainties in arrival time (Reference Lomax, Virieux, Volant and ThierryLomax and others, 2000; Reference Husen, Kissling, Deichmann, Wiemer, Giardini and BaerHusen and others, 2003). Both deep icequakes locate to the southwest of our core network at depths of 100 and 160 m below the ice surface (Fig. 7). They both have a well-defined hypocentre location, as indicated by confidence regions that are ellipsoidal and compact. Moreover, the confidence regions agree well with the 68% confidence ellipsoids showing maximum semi-axes of 20 and 13 m, respectively. The location uncertainty of the western icequake is slightly larger as it locates further away from the network centre, in particular in the radial direction. As the icequakes locate outside the network, azimuthal coverage is poor, resulting in limited resolution in depth. However, a hypocentre location clearly beneath the ice surface is in good agreement with the observed impulsive P- and S-wave arrivals (Reference Walter, Deichmann and FunkWalter and others, 2008; Dalban Reference Dalban Canassy, Walter, Husen, Maurer, Faillettaz and FarinottiCanassy and others, 2013). The icequakes are both located outside the network and clearly separate in their locations but are only 30 min apart in time. Clustering of deep icequakes has been observed on Triftgletscher and Gornergletscher, Swiss Alps (Reference Walter, Deichmann and FunkWalter and others, 2008; Dalban Reference Dalban Canassy, Walter, Husen, Maurer, Faillettaz and FarinottiCanassy and others, 2013). It remains unclear at this stage whether deep icequakes beneath the GrIS will cluster in a similar way to that observed on Alpine glaciers. Detection of deep icequakes requires sophisticated algorithms beyond simple STA/LTA trigger algorithms since surface crevasse icequakes constitute the vast majority of the detected events.
Probabilistic location results of two deep icequakes shown in horizontal plane view (east–north) and two vertical cross sections (east–depth, north–depth). Results are shown as scatter density cloud (green points) and projection of the 68% confidence ellipsoid (blue). The maximum likelihood hypocentre locations (black circles) of the icequakes locate at 100 and 160 m below the ice surface, respectively, and are located outside the network.
Short-duration narrowband tremor
We observe tremor-like signals that last tens of seconds (Figs 2c1 and c2 and 3c), significantly longer than the above-described surface crevasse and deep icequakes. Their energy concentrates in narrow frequency ranges below 10 Hz, and the waveforms lack impulsive onsets and distinct body or surface wave arrivals.
Preliminary analysis of selected events indicates that they originate from outside the seismometer network. Similar events have been observed in volcano seismology (e.g. Reference ChouetChouet, 1996) and various glaciers (e.g. Reference Métaxian, Araujo, Mora and LesageMétaxian and others, 2003; Reference O’Neel and PfefferO’Neel and Pfeffer, 2007) where they were explained as resonances in fluid-filled cracks. This is also a reasonable explanation for our study site in the ablation zone of the GrIS. It indicates that fluids play an important role on glaciers not only for short-duration narrowband tremors but also for the tremor-like coda of deep icequakes (Fig. 2b2).
Long-duration tremor
Visual inspection of spectrograms of the entire continuous seismic dataset revealed hour-long tremor signals (Fig. 2d) typically starting in the afternoon hours. These events occur on 29 days out of the 45-day-long monitoring period. Signal energy is concentrated in the 2–11 Hz range within distinct frequency bands (Fig. 8c). Duration and signal intensity vary between days of observation. The earliest tremor begins at 11:00 UTC. The latest ends at 5:30 UTC. Tremor duration is 4–15 hours, with an average duration of 6 hours. Figure 8 reveals a clear correlation between water level in M1 and start and end times of long-duration tremor: tremor starts while water level in the moulin rises, and stops during falling water level. Therefore, we refer to this type of seismic signal as ‘moulin tremor’.
(a) Water level in M1. Dashed line marks time when no data are available. (b) Seismic waveform filtered between 2 and 5 Hz. (c) Spectrogram of 24 hours of observation showing the frequency content of the tremor. Black rectangle marks frequency and time range used for tremor location. 
marks the band of dominating energy at around 4 Hz. 
marks the frequency range where the energy anticorrelates to water level in the moulin (a).
Seismic records of moulin tremors, bandpass-filtered in the 2–5 Hz range, show a cigar-shaped envelope with no clear onsets or phase arrivals (Fig. 8b) interrupted by frequent icequakes. The concentration of energy in low frequencies (2–6 Hz), smaller tremor amplitudes at the 350 m deep borehole sensor compared to the sensor at 150 m depth (FX13/FX14) and horizontal particle motion concentrated in the radial component suggest that Rayleigh waves dominate the tremor signal. In the spectrogram, energy concentrates in discrete frequency bands. Below ˜6 Hz these frequency bands do not change over time ( in Fig. 8c). Above 6 Hz, however, the frequency content changes over time, with increased energy at the beginning and end of the tremor period and weak energy in between ( in Fig. 8c). As water level in M1 rises (Fig. 8a), these high frequencies quickly decay until the tremor reaches a stable mode, where the frequency band below 6 Hz stays mostly stable. Thus, the tremor spectrogram appears symmetric with respect to its midpoint in time. This multilevel U-shaped temporal evolution of energy during a single tremor episode points to a complex source mechanism including amplification, resonance and absorption of the seismic signal. The symmetry and the change in frequencies is apparent to a varying degree in all moulin tremors.
Before we can discuss possible mechanisms for the observed moulin tremor and their correlation with the water level measured in M1 we need to confirm that the observed tremor originates close to M1. This is done in the next section by employing a grid search location algorithm with an amplitude decay model.
Moulin Activity as Source of Long-Duration Tremor
As the moulin tremor lacks clear onsets or phase arrivals, we cannot employ traditional seismic source location methods based on arrival time measurements. Consequently, we determine tremor epicentre using the decay of signal amplitude attenuation with distance. This technique was developed for the location of rockfalls, long-period events and eruption tremor sources at Piton de la Fournaise volcano, La Réunion, Indian Ocean (Reference Battaglia and AkiBattaglia and Aki, 2003), and recently adapted for application on glaciers by Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others (2013) (Russell Glacier, Greenland). Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others (2013) estimated location uncertainties with the help of seismic reflection shots. As we did not conduct calibration shots in our experiment, first a synthetic data approach is needed for testing the appropriateness of the procedure and for error estimation.
Tremor location by amplitude decay model
The amplitude decay of a seismic signal measured at different stations is primarily a consequence of geometrical spreading and of the anelasticity of the medium neglecting any local variations. It can be described by Eqn (1) approximating the study volume by a homogeneous half-space model (Reference Battaglia and AkiBattaglia and Aki, 2003; Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others, 2013):
Equation (1) expresses geometrical spreading and anelasti-city of the medium as a power-law and exponential dependence on distance ri
, respectively. Here, A(ri) is the seismic amplitude measured at station i at distance ri
to the source. Equation (1) relates A (ri) to the amplitude A
0 at the source and to the anelastic attenuation coefficient , depending on the seismic quality factor Q and wave velocity β valid for a certain frequency of interest f. Seismic quality factor Q accounts for the fractional loss depending on properties of material between source and receivers (e.g. Reference Lay and WallaceLay and Wallace, 1995). Parameter n refers to geometrical spreading for body waves (n = 1) or surface waves 
Local tremor amplitude derivation
The amplitude of the tremor signal cannot simply be characterized by the maximum amplitude at a specific time as different phase arrivals cannot be picked and correlated between the stations (Fig. 2d). Rather, we calculate the representative tremor amplitude by averaging over a 30 min time window during a strong and constant tremor signal at as many seismometers as possible. In order to smooth the waveform, the root-mean-squared (RMS) envelope is calculated for non-overlapping windows of 60 s. This 1 min window is chosen as a compromise between reduction of noise and smearing of information. The envelope function is calculated for the 2–5 Hz bandpass-filtered signal (two-pole Butterworth filter) and its Hilbert transform (Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others, 2013). This frequency band exhibits only marginal variation during tremor episodes (black rectangle in spectrogram, Fig. 8c). The observation A ðri Þ is determined by calculating the mean value of the RMS envelopes for 30 min of consecutive 1 min windows. In this study, only measurements of the surface seismometer LE3D of the core network (FX01–FX09) are used, because they have the same response function and were installed in a similar way, i.e. similar coupling to the ice. A longer time window than 30 min is not feasible as each station underwent regular maintenance work introducing brief interruptions in recording.
Figure 9 shows an example of the calculated RMS waveform envelope during tremor activity recorded at all stations. It suggests that the epicentre lies closest to station FX01, and, hence, near moulin M1, because FX01 shows the highest RMS envelope. The time-dependent amplitudes show a very high correlation between different stations and thus reveal that they are dominated by energy emitted from the same source. Scattering and three-dimensional path effects are comparatively minor as the RMS envelopes change uniformly for all stations. Hence, for source location procedure we may safely approximate the study volume by a homogeneous half-space model as implied by Eqn (1). In order to compare different signal strengths for further computations, RMS envelopes are normalized by those of the strongest station. We calculate the ratio between the station with highest amplitude (closest station) and every other station for each 1 min window. Ideally, if we could separate the observed signal into tremor signal and ambient noise (including, e.g., icequakes and site effects), we could use the former for epicentral distance calculation and the latter to estimate the observation uncertainty. If the amplitude ratio for a specific station relative to the signal at the epicentre-nearest station remains constant, while the absolute amplitude varies slightly over time, we interpret these relative amplitude variations as noise and use this noise as observation uncertainties. With the estimated uncertainties for all our epicentre locations, we want to see if we can assign the tremor signal unambiguously to the location of one particular moulin. Therefore, we use the maximum observed variation of amplitude ratio of 9% as a conservative estimate of the average observation error in our synthetic tests.
RMS envelopes (1 min windows) calculated during high tremor activity for all stations. RMS envelopes show a strong correlation between different stations, indicating that the signals come continuously from the same source with slight changes in emitted energy.
Grid search location algorithm
We derive the epicentre location of the moulin tremor using in Eqn (1), as our tremor is dominated by surface waves. Since wave velocity β, quality factor Q and source emitted amplitude A
0 are unknowns, when locating the epicentre of the tremor source we effectively solve a coupled inversion problem. We optimize for the model parameters (α (β, Q, f), A
0) and the two epicentre parameters (lat./long.), that denote r towards zero. As documented in several studies, Q and wave velocity can change significantly depending on type of glacier, temperature, water content, frequency and water saturation (e.g. Reference KohnenKohnen, 1974; Reference Gusmeroli, Clark, Murray, Booth, Kulessa and BarrettGusmeroli and others, 2010; Reference Peters, Anandakrishnan, Alley and VoigtPeters and others, 2012). All further estimations of Q are made with wave velocity for surface waves of 1650 m s−1 taken from experiments on Alpine glaciers (Reference Roux, Walter, Riesen, Sugiyama and FunkRoux and others, 2010). A frequency f of 3.5 Hz is chosen because it is the centre frequency of the analysed frequency band of 2–5 Hz. Using our RMS envelope A ðri
Þ, the epicentre is located by a double-nested grid search procedure following Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others (2013). We compare sets of observed and calculated amplitudes at each station i by means of a least-square misfit function. The two-step grid search procedure is more adept at finding the global minimum. The solution for epicentre is found by curve fitting (see Eqn (1)) of the observed amplitudes with the Levenberg–Marquardt algorithm (Reference MoréMoré, 1978) where the final location of the epicentre is chosen with the minimal residuals of misfit calculated as L-2 norm.
Assessing location uncertainties, we performed a series of tests with synthetic data mimicking the real situation at hand but with a priori known source location and seismic model parameters as routinely done in earthquake tomography (e.g. Reference KisslingKissling, 1988). The synthetic source location has been chosen inside the network in the vicinity of moulin M1 in order to simulate the real dataset as closely as possible. Only the most reliable eight stations (FX01–FX08) are used for processing. Station FX09 ran for only 50% of the deployment period and hence was not included in the synthetic tests. Theoretical amplitudes A (ri ) for each station are calculated with Eqn (1) using specific choices of model parameters A 0, Q, f and β.
Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others (2013) estimated Q = 35 for surface waves for a frequency of 20 Hz on the Greenland ice sheet, and Reference Métaxian, Araujo, Mora and LesageMétaxian and others (2003) estimated Q = 3.4 for water for low-frequency events (<5 Hz) on the ice cap of Cotopaxi volcano, Ecuador. Since for our study site the exact Q value is not known, first we tested the location procedure and its sensitivity to Q values with synthetic tests for variable Q values in the range 2–35. Figure 10 shows the result of the inversion with different Q values in the top row without error and in the bottom row with a large error added to station FX06. We note that while in all cases the region of minimal normalized residuals encompasses the true location, the size of this region systematically increases with increasing Q value. With increasing Q value, the exponential part of Eqn (1) representing the inelasticity of the medium rises in weight of total attenuation of the amplitude. Due to station geometry, the shape of this region of equal likelihood for source location is non-elliptic. We further observe an accuracy estimate of 20–40 m for single large observation error (of 9%) added to FX06, which is a crucial station in the location procedure.
(a–c) Synthetic location tests with different Q values as indicated in the forward model for data without observation error. (d–f) Synthetic location tests with different Q values and an error of 9% added to amplitude for station FX06 to study effects of large error for a critical station. White cross marks the synthetic source location (forward model), and the star marks the derived solution with offset dR. The normalized L-2 norm of residual (colour-coded) marks area of similar fit quality. Changes in Q value do not affect the capability of finding the absolute minima within the study area for epicentre location.
To obtain a quantitative uncertainty estimate of location procedure we performed a Monte Carlo simulation. For this test, synthetic amplitudes have been calculated with Q = 4 as this value denotes the best average fitting of all real moulin tremor data and our test documented no significant dependence of epicentre location on Q. Finally, Gaussian distributed observation errors in the order of the above estimated maximum real observation error (9%) are added to obtain synthetic datasets ready for inversion. One hundred such datasets with variable noise added were calculated and inverted for source locations. The results are visible in Figure 11 as a point cloud of blue dots. The circumcircle encompassing all resulting locations has a radius extension of 63.9 m. The slight asymmetric clustering of resulting locations is explained by the geometrical distribution of the stations around the synthetic source. However, we calculate the circumcircle as we neglect the influence of the station geometries in order to assure having an uncertainty including worst cases. From these results we derive an uncertainty estimate of 65 m in latitude and longitude for any single location.
Results of Monte Carlo simulation to estimate location uncertainties. One hundred simulations (blue dots) were computed using f = 3.5 Hz, a surface wave velocity β = 1.65 km s−1 and Q = 4. True simulated source location is marked with a black cross. The circumcircle including all epicentres has a diameter of 127.8 m.
Tremor source location
With the described method it was possible to locate 23 out of 29 detected tremors (Fig. 12). Six tremors could not be located due to low signal-to-noise ratios or due to a low number (<5) of recording stations. Grey dashed circles around each location represent the maximum 65 m of uncertainty. The epicentres of all tremors cluster around the major moulin (M1) inside the station network. The precise extent of this moulin is unknown, though at least 5 m in width and 10 m in length. The moulin positions (black crosses) were measured with a handheld GPS at some distance from the moulin as it is dangerous to go too close to the edge of a moulin. The horizontal location of the moulin entrance is estimated with a precision of about ±20 m (visualized with black circle in Fig. 12). Little is known about the moulin subsurface geometry; it is most likely a chain of waterfalls with intermediate horizontal conduits or chambers.
Epicentre location results (blue dots) of all locatable moulin tremors using the amplitude decay location method. Grey circles show expected uncertainty of ±65 m as derived from synthetic tests. Black crosses mark the positions of M1 and M2 with their expected uncertainty (black circle). Within their uncertainties all moulin tremors locate close to M1 inside the seismic core network.
All inverted locations except one are located within the accuracy of the moulin location (Fig. 12). The one outlier distinctively southeast of the other locations belongs to a day when station FX01 had higher noise than usual and FX05 was missing due to outage. However, we note that this tremor signal is clearly not connected to any other moulin in the region and can still be assigned to M1 with an estimated uncertainty of ±65 m. Moulin M1 inside the network was equipped with the pressure sensor, and the correlation between water level and tremor visible in Figure 8 is therefore no coincidence.
The approach used in Eqn (1) revealed a quality factor Q between 2.2 and 7.7. This is rather low compared to other experiments (e.g. Reference Peters, Anandakrishnan, Alley and VoigtPeters and others, 2012; Reference Jones, Kulessa, Doyle, Dow and HubbardJones and others, 2013). However, we are analysing low frequencies (f = 3.5 Hz), and our calculated Q values also directly depend on the (unknown) phase velocity (α (β, Q, f) is estimated). The maximal change between calculations of Q with phase velocity β = 1500 m s−1 and 1800 m s−1, respectively, reveals a difference of 1.4 for the smallest α derived by the fitting procedure. Therefore, we conclude, that the sensitivity to phase velocity β is small. Our revealed Q values (with β = 1650 m s−1) are similar to results observed on Cotopaxi volcano (Reference Métaxian, Araujo, Mora and LesageMétaxian and others, 2003) for water-filled cracks. Thus, we assume that our bulk Q value is influenced by site effects of crevasses, possible water or air cavities and inhomogeneity in the ice properties lowering the Q value. Nevertheless, the absolute Q value is insignificant for our epicentre location at M1 due to its low influence on the absolute location (Fig. 10).
Discussion
We have presented a variety of seismic signals from the GrIS (Fig. 2), whose sources are likely directly or indirectly related to the presence or flow of water. The moulin tremor is the most prominent signal, and its epicentre can be reliably located within a few tens of meters of M1. This moulin tremor can therefore be regarded as a seismic signature of large-volume englacial water flow whose identification was the main goal of the present study. To our knowledge, this type of seismic signal has not yet been reported in the literature. A possible reason may be that advances in instrumentation have only recently permitted recording of continuous seismic data in difficult terrain such as glacial ablation zones.
At this point we can only speculate about the source mechanisms of the moulin tremors. A potentially important fact is that we did not observe a tremor signal that can be associated with the smaller moulin M2 at the edge of our seismic network. This suggests a critical inflow, moulin conduit size and/or geometry is necessary to generate a sustained tremor signal above the background noise. In view of the observed relationship between frequency content of moulin tremors and moulin water level (Fig. 8), the analogy to volcanic tremor should be noted. Interaction between solid walls of magma chambers and cracks, fluids and gases has been proposed as the source mechanism of volcanic tremor (e.g. Reference JulianJulian, 1994; Reference ChouetChouet, 1996; Reference ChouetChouet and others, 1997; Reference Lesage, Glangeaud and MarsLesage and others, 2002). Equivalently, we suggest that in our case, the water-filled moulin acts as a resonating body for acoustic and seismic waves possibly triggered by a waterfall inside the moulin. Changes in water level within resonating chambers as well as changes in air–water mixture may explain the temporal changes in tremor spectra.
Before applying numerical models of water resonances, hypocentres of moulin tremors as well as changes thereof should be constrained. Our location procedure via amplitude attenuation is not suited for this task, as it is based on surface waves, which cannot resolve source depth. Alternatively, one could employ matched filter techniques using body waves (e.g. Reference Corciulo, Roux, Campillo, Dubucq and KupermanCorciulo and others, 2012), which we plan to apply in the near future.
Tremor spectra changed not only during a single tremor episode but also over the course of our study period. Interpreting these temporal changes in terms of moulin geometry may provide important insights into the morphology of englacial water channels. Geometries of englacial channels are rarely at static equilibrium, but permanently evolve via the two competing mechanisms of melt enlargement and creep closure (Reference NyeNye, 1953). Temporal changes in tremor spectra would therefore highlight the relative efficiency of these processes on both small (hourly) and seasonal timescales. Consequently, future investigation should focus on moulin tremor sources.
Like moulin tremor, certain icequake sources can also be related to meltwater. So far, we have located only two icequakes with maximal likelihood depths 100–160 m below the ice surface. We suggest hydrofracturing as a source mechanism for these events for several reasons: First, reducing the effective pressure, water can induce tensile fracturing even at depths where the ice overburden usually inhibits crevasse formation (e.g. Van der Reference Van der VeenVeen, 1998). Second, one of the events is followed by a low-frequency coda indicative for resonances of a water-filled cavity (Reference West, Larsen, Truffer, O’Neel and BlancWest and others, 2010). And third, the waveform pattern corresponds to icequakes at intermediate depth detected by Reference Walter, Clinton, Deichmann, Dreger, Minson and FunkWalter and others (2009) on Gornergletscher, which are characterized by moment tensor inversion as tensile fractures.
On the other hand, surface crevasse icequakes are only indirectly related to meltwater. As shown by the blue and red histograms in Figure 5, their activity tends to be highest when surface melt reaches a maximum. Furthermore, flow velocity of the glacier likely increases with melt, which promotes fracturing (e.g. Reference Walter, Deichmann and FunkWalter and others, 2008). We observe strong similarity in the structure and timing of moulin water level, RMS envelope and icequake signal strength (Fig. 13; 26–31 July 2011). The signal strength in Figure 13a is a zoom of Figure 4 illustrating the daily fluctuations in minimum signal strength. Absence of weak icequakes is clearly visible in the afternoon/evening hours. The RMS envelope (Fig. 13b) for the continuous seismic record is calculated in the dominant tremor frequency range and identifies four distinct tremor episodes. The comparison between icequake signal strengths and RMS envelope shows that non-detection of weak icequakes coincides with tremor episodes. Consequently, when no tremor activity is observed, the threshold for non-detection is significantly lower throughout the day. Therefore, we conclude that surface melt in the afternoon hours increases the seismic noise level; however, the moulin tremors influence the detection capability of our network critically. As diurnal fluctuations in icequake detection capability have also been reported on Gornergletscher (Reference Walter, Deichmann and FunkWalter and others, 2008), our observations may indicate that Alpine glaciers also produce tremor signals influencing the seismic background noise.
(a) A zoom into 5 days of observation shown in Figure 4, with each point representing a triggered icequake with corresponding signal strength. (b) RMS envelope of the waveform filtered between 2 and 5 Hz. Each rise in envelope corresponds to a moulin tremor with different intensity and duration. (c) Moulin water level and water level measured in the stream feeding into moulin M1. Higher water level in the moulin coincides with the occurance of the moulin tremor, but the correlation with the water level inside the stream has lower influence on signal strength than high moulin water level.
Figure 13c is a comparison of the water level in the moulin and in the surface stream, with the latter representing meltwater availability at the ice surface. The stream water level shows a continuous diurnal pattern, whereas the moulin water level rises only occasionally above the pressure sensor. Possible explanations of this are temporal changes in the complex drainage system consisting of crevasses and subglacial cavities as well as subsurface streams feeding into the moulin (e.g. Reference McGrath, Colgan, Steffen, Lauffenburger and BalogMcGrath and others, 2011; Reference Gulley, Grabiec, Martin, Jania, Catania and GlowackiGulley and others, 2012). In any case, the smaller influence of diurnal stream water level and RMS envelope suggests that the surface water availability has a relatively minor influence on seismic background noise. Instead, higher RMS envelopes correspond directly to tremor activity and show daily changes in duration and intensity. They are directly correlated to high moulin water level, and their spectral characteristic is most likely also related to inflowing water and water level inside the moulin.
Conclusion
We successfully operated a dense surface and borehole seismic network on the GrIS in boreal summer 2011. Regular maintenance and re-levelling of the stations ensured high-quality data with few data gaps. Our continuous data contain a large variety of seismic signals, including shallow icequakes related to surface crevassing, intermediate-depth icequakes (100–160 m beneath the ice surface) related to hydrofracturing, and tremor lasting from several minutes to several hours. Qualitatively, our data confirm similarities to the seismic signal variety measured on volcanoes (Reference West, Larsen, Truffer, O’Neel and BlancWest and others, 2010). Our icequake catalogue includes seismic signals whose sources are well understood from previous studies (e.g. Reference Deichmann, Ansorge, Scherbaum, Aschwanden, Bernardi and GudmundssonDeichmann and others, 2000; Reference Walter, Clinton, Deichmann, Dreger, Minson and FunkWalter and others, 2009). The moulin tremor, on the other hand, is a new observation relating seismic signal to englacial water flow. The tremor-generating moulin constitutes a main hydraulic connection between the glacier surface and its bed. With a total of 150 hours of tremor signal, we plan to study the conditions necessary to generate the observed characteristics of the tremor spectrum and the triggering mechanisms of such long-duration signal. From these studies we expect to gain a more in-depth understanding of tremors on glaciers and, hence, of englacial and subglacial water flow. Measuring englacial and subglacial water flow remains a challenge. However, seismic monitoring in combination with other glacial observations (such as done in the ROGUE project) broadens the prospects for studying hydrologic processes and provides valuable insights into the inaccessible regions of the ice sheet.
Acknowledgements
The project was funded by Swiss Federal Institute of Technology Zürich (ETH) under grant ETH-27 10–3. The salary of F.W. was partially funded by the European Union Seventh Framework Programme (FP7-PEOPLE-2011-IEF) under grant agreement No. 29919. This work was partially supported by grants from the Swiss National Science Foundation (200021_127197 SNE-ETH (ROGUE)) and the Arctic Natural Sciences Program of the US National Science Foundation (OPP-0909454 and OPP-0908156). The satellite image captured by the WorldView-2 satellite in Figure 1 was provided by Polar Geospatial Center. Seismogram picking was done with the software package SNAP written by M. Baer. Waveform processing and plots were done with MathWorks MATLAB, Obspy (A Python Toolbox for seismology/seismological observatories) and ESRI ArcMap. We are grateful to many members of the Institute of Geophysics, the Swiss Seismological Service and the Laboratory of Hydraulics, Hydrology and Glaciology of ETH. We also thank Matthew J. Hoffman for providing an additional dataset within the ROGUE project for further analysis. The reviews of the two anonymous reviewers were extremely helpful for improving the quality of the paper. We also thank M. Funk, C. Ryser, T. Wyder, C. Senn, K. Plenkers, S. Hiemer, M. Meier and A. Bauder for their help during the field campaign in Greenland and making it all possible.
