Seasonal variations of glacier seismicity at the tongue of Rhonegletscher (Switzerland) with a focus on basal icequakes

Abstract

Fractured systems developing at the base of temperate glaciers play a significant role in the transport of subglacial water. Their influence on dynamics and hydraulics of glaciers and the seasonal alterations they experience during the hydrological year, remain however, poorly understood, as they lay beyond the reach of most conventional glaciological techniques. In this work, we investigate the seasonal evolution of glacier seismicity in the ablation zone of Rhonegletscher (Switzerland) and focus on basal fracture processes. To do this, we use seismic measurements performed between June 2012 and July 2013 and complementary measurements of surface motion and changes in melt water input inferred from degree day modeling. Results of a spectrogram analysis highlight predominant seismic emissions until end of October 2012 and from mid April 2013 in line with measured surface velocities, and indicate a rapid adaptation of the basal drainage system at initiation and end of the winter period. The analysis of icequakes emitted from two active regions shows that the basal events, due to opening and closing of cracks within fracture networks, persist after winter initiation, and points out the year-long preservation of the associated sources. In this regard, we deduce that basal crack networks aseismically evolve during the winter period.

1. INTRODUCTION

Subglacial water can have a profound impact on glacier and ice sheet dynamics. Depending on the capacity of the subglacial drainage system, incoming water masses may overwhelm or underwhelm the subglacial drainage system and thus lead to an increase or decrease in subglacial water pressure. As a result, basal sliding is either promoted or suppressed (Iken and Bindschadler, 1986; Anderson and others, 2004; Schoof, 2010).

Beneath ablation zones of temperate glaciers, subglacial drainage systems are usually characterized by a seasonal cycle. During winter an inefficient, high-pressure subglacial drainage system exists. In contrast, during summer months, surface melt water reaches the glacier bed via crevasses. At the glacier bed, this melt water flows down the hydraulic gradient dissipating part of its potential energy into heat. As a result, the ice-walled water channels quickly enlarge their diameter forming efficient drainage pipes, which operate at low water pressures. These efficient ‘R-channels’ (Röthlisberger, 1972) are only overwhelmed during sudden bursts in surface run-off produced by, e.g. warm days, rain events or the sudden drainages of ice surface lakes (Das and others, 2008). In such cases the channels inject water into the laterally linked cavity system, and thereby promote glacier sliding. The channels diminish when water input fades and creep closure overcomes melt enlargement by water flow (Iken and Truffer, 1997).

The idea of R-channel evolution explains a wide range of ice flow measurements on polar (Bartholomew and others, 2010) and non-polar glaciers (Werder and others, 2013). However, there also exists evidence that additional non-negligible water flow takes place in linked fractures near the glacier base. Such fractures can occupy up to 1% of the entire ice volume (Harper and others, 2010) and host water flow from 0.5 to 4 cm s⁻¹, which is negligible compared with turbulent R-channel flow (Fountain and others, 2005). Nevertheless, water transport through the linked-fracture drainage system may still play an important role in glacier hydraulics. In particular, it may provide an initial water flux needed for the formation of new R-channels. In this case, R-channels would not only form, where water input or basal topography (Gulley and others, 2012) favor large subglacial water fluxes. Rather, the presence of fractures near the glacier bed could be a sufficient condition for the formation of an efficient subglacial drainage system.

Due to the confining ice overburden pressure, the formation and growth of fractures near the glacier base requires pressurized water (Van der Veen, 1998). Therefore, the linked fracture drainage system most likely evolves via sudden hydro-fracture episodes (Das and others, 2008). In contrast, the slow water flows do not lead to significant melt enlargement typical for R-channels (Fountain and others, 2005). Details about the dynamics of basal fracture drainage remain elusive. For instance, it is not clear if basal fracturing is a widespread phenomenon or if there exist specific fracture birthplaces favored by local topography or strain rates. Once opened, do these fractures exist for periods exceeding a single season or do they form and diminish more frequently? If fractures remain open for longer time periods, they can be advected with glacier flow and influence glacier drainage far away from their birthplaces. Furthermore, in case of stored water, basal hydro fracturing may occur even in the absence of surface melt. This again contrasts with channelized drainage systems, whose morphology is most dynamic during times of elevated water flux. Answering these questions will provide new insights into the formation and evolution of subglacial drainage pathways. However, the answers also require a means to observe basal fracture processes, which are beyond the reach of most conventional glaciological techniques.

In this concern passive seismic investigations represent an innovative technique for continuous monitoring of near-surface, intermediate and subglacial environment. This has been done on various glaciers in the Alps (Deichmann and others, 2000; Roux and others, 2008; Walter and others, 2008; Helmstetter and others, 2015), in Alaska (O'Neel and others, 2007; Bartholomaus and others, 2008; West and others, 2010; Mikesell and others, 2012), in Greenland (Rial and others, 2009; Jones and others, 2013; Roeoesli and others, 2014) and in Antarctica (Wiens and others, 2008; Winberry and others, 2009; Zoet and others, 2012). Among these studies, seismometers installations in direct contact with glacier ice have provided new insights into glacier dynamic and fracturing at all depths. However, restricted investigation periods have typically affected the conclusiveness of these studies. More particularly, harsh climatic conditions and abundant summer surface melt on glacier ablation zones inhibit extended monitoring periods. Similarly, during the winter season the deployed seismic equipment remains under a snow cover for weeks to months and cannot be visited frequently. The limited availability of long-term seismic records from on-ice installations thus demands further investigation of seasonal variations in glacier seismicity, in particular related to changing melt conditions.

Here we present a multi-seasonal seismic analysis performed between May 2012 and July 2013 at the tongue of the temperate Rhonegletscher (Switzerland). Our work addresses two main points: first we give a year-long overview of the continuous glacier seismicity, i.e of all kinds of seismic signals recorded in the study area, and discuss the observed seasonal patterns; second we constrain the role of basal fracture processes in the dynamics and hydraulics of the glacier by focusing on near-bedrock icequakes indicative of opening and closing of basal fractures. To our knowledge we produce the first year-long glacier seismicity record on the ablation zone of an alpine glacier.

2. STUDY SITE

Rhonegletscher is located at the eastern end of Valais in the Swiss Alps (Fig. 1a). It flows from 3556 m a.s.l. to 2250 m a.s.l. and covers a surface of ~16 km² (Farinotti and others, 2009a). Over the past 10 a, Rhonegletscher's tongue retreated substantially and a proglacial lake formed at 2208 m a.s.l. (Tsutaki and others, 2013), acting as a spring of the Rhone river. The area of investigation chosen in this study ranges from 2450 to 2250 m a.s.l. (Fig. 1a), with moraines on both east and west sides and the proglacial lake ~250 m downstream.

Fig. 1. (a) Aerial view of the study area at the tongue of Rhonegletscher (coordinates system CH1903/LV03) and situation map in the small frame. The dashed red marks delimit the surface considered for the icequakes location, sized to 600m × 750 m. Ice thickness was obtained from helicopter-based radar survey. Automatic GPS measurements were performed at two stakes s23 and s33 at daily or sub-daily time scale during the investigated period. (b) Zoom of the study area with identical layout. Colored arrows depict horizontal surface motion measured between various time intervals indicated in the legend. For more clarity, displacement values are enhanced by a factor of 20. In both frames red and green contours depict the two basal active regions R1 and R2.

3. FIELD MEASUREMENTS

A number of seismometers ranging between 3 and 13 as well as 2 GPS stations were continuously operated in the study area from beginning of May 2012 to end of July 2013. Sampling characteristics and associated monitoring periods are summarized in Figure 2. Indicated monitoring periods correspond to time when sensors were installed on the glacier and do not take into account data gaps. Note that seismic data recorded before 10 June 2012 were not considered in this study.

Fig. 2. Sampling characteristics for seismometers and GPS measurements according to monitoring periods. Active sensors are referred with names (R1 and R2) used in Figure 1a. Note that data recorded before 10 June 2012 were not considered in this study.

3.1. Seismic monitoring and blasts

Between 10 June 2012 and 25 July 2013, three seismic sensors were installed into 2 m deep boreholes at altitudes between 2450 and 2250 m a.s.l. (RA11, RA12, RA13; Fig. 1a), for maximal north and east extents of 548 and 356 m, respectively. The instrumentation consisted of three 3-component Lennartz (LE3D) seismometers with 1 Hz eigen frequency and flat response on the 1–80 Hz interval linked to Nanometrics Taurus data loggers. Power supply was provided by batteries coupled with solar panels, and 3 components ground velocity (vertical and two horizontal components) was continuously recorded on memory cards with 500 or 250 Hz sampling rate (Fig. 2). In principle, the system had an autonomy of ~3 months assuming solar panels remained free of snow. Maintenance operations were performed at least once every 3 months. An additional array composed of ten LE3D with 1 Hz eigen frequency, a flat response on the 1–80 Hz interval and positioned on tripods at the ice surface, was set up in the study area from 28 August to 2 October 2012 (RA01 to RA10; Fig. 1a). Its altitude ranges between 2352 and 2294 m a.s.l., with maximal north and east apertures of 578 and 567 m, respectively. We used the same data recording systems as mentioned above, with a 500 Hz sampling rate. The array benefited from daily or sub-daily maintenance visits.

In order to constrain englacial and subglacial seismic velocities, nine explosive charges (Riodin; referred to bl1 to bl9) from 0.5 to 8 kg were detonated within boreholes at depths of 2 m (bl2, bl4–bl9), 100 m (bl1) and 130 m (bl2) (Fig. 3) when the 13 sensors array was operational.

Fig. 3. Aerial view of the study area (coordinates system CH1903/LV03) with surface and borehole sensors denoted by green triangles and blue sensors, respectively. Sensor positions correspond to coordinates measures in September 2012. Red and cyan points depict horizontal projection of the probability density functions (pdf) associated with surface and near bedrock blast relocations (bl1 to bl9; Section 3.1), respectively. (b) Vertical transversal cross sections of pdf relocations at blasts epicenters. Blue curves are glacier surface and bedrock topography. Horizontal and vertical units are indicated in the top frame. Blue and yellow stars identify maximum likelihood point and true blast location, respectively. Epicentral and vertical discrepancies are indicated in (a) and (b) together with the blast labels.

3.2. Melt and accumulation modeling

Snow accumulation and ice/snow melting were locally computed at point s23 (Fig. 1a) over the entire study period using an enhanced temperature-index melt model (Pellicciotti and others, 2005). This model includes the shortwave radiation balance and distinguishes between energy supplied by solar radiation and sensible heat flux. Model input data were taken from the MeteoSwiss weather station, Grimsel Hospiz (1980 m a.s.l.) located ~4 km west of the study area. Collected parameters were air temperature, precipitation, global radiation, humidity, wind speed and wind direction in hourly resolution. Differences in elevation were compensated applying precipitation and temperature lapse rates (Gabbi and others, 2014). Model calibration was achieved using in situ mass balance measurements, and temporal extent of modeled snow cover was also verified with the help of daily photographs of the study area from an automatic camera.

3.3. Glacier and bedrock geometry in the study area

In order to approximate the glacier's surface geometry we use a DEM derived from aerial photogrammetry (Swisstopo pictures) (Fig. 1a) taken in September 2012. At a grid spacing of 5 m the estimated uncertainty is ±0.2 m in height. Bedrock topography was determined from helicopter-based radar surveys carried out in spring 2008 and winter 2012. The data were acquired from the University of Münster Airborne Ice Radar (Blindow, 2009). We subsequently applied the Ice Thickness Estimation Method in order to interpolate bedrock topography between radar measurements. This technique makes use of surface data, which are inverted by means of an ice flow model. The associated uncertainty is ±15% of the ice thickness (Farinotti and others, 2009b).

The results (Fig. 1a) reveal a bedrock overdeepening in the upstream part of the investigated area, with a progressive ice thinning along the flow direction, relatively high slopes on the glacier banks and a maximal ice thickness up to 120 m in the central part of our study area (Fig. 1a).

3.4. Surface motion

Glacier surface motion in the study area was derived from automatic differential GPS measurements (Leica System 500) at two stakes (s23 and s33) set up at respective altitudes of ~2330 and 2270 m a.s.l. (Fig. 1a). Stake s23 was operational from 14 June 2012 to 4 July 2013 and s23 from 20 June to 2 October 2012. Both resulting datasets exhibit time gaps, due to power supply issues or too large uncertainties. We measured positions at two sampling rates: twice a day at 0700 h and 1900 h (referred as static measurements) from October 2012 to July 2013, and every 2 min (referred as kinematic measurements) during the rest of the investigated period (Fig. 2).

The Static dataset was processed using the Leica's software and resulting positions show estimated accuracy of 3 × 10⁻³ and 5 × 10⁻³ m in the horizontal and vertical direction, respectively. Kinematic data were retrieved and processed using the Track software (Chen, 1998; King, 2004) according to the procedure followed in Riesen and others (2011). Results present the same uncertainty ranges as mentioned above. Stake position information is available from 19 June 2012 to 4 July 2013 for s23 and from 19 June to 2 October for s33. Additional measurements of positions of the 3 borehole sensors were also manually conducted in May, September and October 2012 using a Leica Viva system offering comparable accuracy.

4. SEISMIC DATA PROCESSING

The seismic data processing has three steps: (1) computation of spectrograms using the time series obtained by the reduced array composed of RA11-RA13 (Figs 1a and 2). These spectrograms encompass the whole glacier seismicity (i.e. the spectral energy of all the seismic signals recorded in the study area including non-glacial events); (2) the location of some icequakes detected from 28 August to 2 October 2012, when the dense array composed of RA01–RA13 was maintained and selected as candidate template location; and (3) the investigation of the year-long continuous dataset obtained by the reduced array using cross correlation techniques.

From 10 June 2012 to 25 July 2013, rates of data completeness for RA11, RA12 and RA13 were 59, 70 and 28%, respectively (Figs 7a and 6). Reasons for such data gaps were mostly power supply deficiency and data corruption. However, we note that there was permanently at least one operational borehole sensor from 10 June to 28 December 2012 and 12 January to 25 July 2013.

4.1. Daily spectrograms computation

Daily spectrograms were computed for each of the 3 borehole sensors using windows of 2¹² samples and an 85% overlap. Frequency ranges between 1 and 125 Hz, equal to the Nyquist frequency. Note here that the 136 dB dynamic range Lennartz 1 Hz used in this study has a corner frequency of 100 Hz, implying a loss of sensitivity for values >115 Hz when a 250 Hz sampling rate was chosen.

4.2. Icequakes location

4.2.1. Probabilistic nonlinear procedure

A nonlinear probabilistic procedure approach as implemented in the software package NonLinLoc (NLLoc) (Lomax and others, 2000) was applied to locate icequakes in our study. NLLoc computes the posterior probability density function (pdf) of seismic source location as defined by Tarantola and Valette (1982) and Moser and others (1992). The solution accounts for location uncertainties due to the geometry of the network, measurement errors of the observed arrival times and errors in the calculation of the theoretical travel times in ice (Lomax and others, 2000; Husen and others, 2003; Husen and Hardebeck, 2010). The posterior pdf is computed using the Oct-Tree importance sampling algorithm (Curtis and Lomax, 2001), which provides an efficient and reliable sampling of the solution space by looking for the maximum probability of location (Husen and others, 2003). Dalban Canassy and others (2013) provide further details on probabilistic icequake location. Location results are represented by means of a maximum likelihood hypocenter location and associated pdf describing the location uncertainty. We used manually picked p- and s-wave arrival times of investigated events as input data in the inversion procedure, with measurement error individually estimated for each picked value in the range ±1 to ±4 ms.

Linear regression between differences in arrival times of surface blasts measured at each of the sensors and epicentral distances determined a p-wave velocity of 3.33 km⁻¹ ± 0.08 km⁻¹, in agreement with values found in previous studies (Gischig, 2007; Dalban Canassy and others, 2013). For s-wave velocity we used 1.8 km s⁻¹ according to Deichmann and others (2000).

4.2.2. Velocity model validation

Theoretical ice travel times in the study area were computed with a homogeneous velocity model. The relevancy of such a simple model for investigating sources of icequakes detected in our study area was assessed by relocating the above mentioned near bedrock and surface blasts. This was performed using arrival times of direct p-waves picked at the 13 sensors. Results are shown in Figure 3. Blast relocation accuracy is expressed by means of epicentral (Fig. 3a) and vertical discrepancies (Fig. 3b) between maximum likelihood points of relocated sources and true position. The associated precision is depicted by the three-dimensional (3-D) volume described by the resulting pdf, which is bounded by the 100% confidence level.

Blast relocations exhibit epicentral and vertical accuracy ranging from 1.6 to 6.5 m and 2.5 to 5.1 m, respectively. Such high values can be explained by the dense array coverage, by the large number of sensors used for the relocation and by the high signal to noise ratio characterizing blast waveforms. The best epicentral accuracy is found in the upstream part of the study area, whereas relocation in the vicinity of RA10 and RA13 point out the biggest errors (6.5 m). Observed vertical differences remain below 10% of the ice thickness measured at the source (Fig. 1a), for both the near surface and near bedrock blasts.

Best epicentral precisions are found in the middle of the array, where pdf shapes closely resemble traditional ellipsoidal uncertainties. On the contrary, pdfs computed for relocation of bl5, bl8 and bl9 present more non-ellipsoidal shapes resulting from azimuthal station coverage. Vertical pdf extent present a close constraint of hypocenter locations, which allows us to unambiguously locate basal events in the study area.

The accuracy and precision of blast relocations justify the use of a homogeneous velocity model in our study area to locate detected icequakes. More particularly, such a model provides vertical constraints on the epicenters within <10% of the ice thickness, which lies well within the bedrock depth uncertainty of 15%. In this context, we conclude that our homogeneous velocity model offers a reliable means of identifying basal icequakes.

4.3. Definition of deep icequake dataset

As a first step, icequakes were detected in the continuous seismic data recorded with the 13 sensors array from 28 August to 3 October 2012. This was achieved using a standard short-term window over long-term window average trigger algorithm (STA/LTA trigger algorithm; e.g. Allen, 1978) with window lengths of 80 and 800 ms, respectively. Trigger threshold was set to 3, and an event was declared if the trigger threshold exceeded on the vertical component of at least five stations within a 1 s time period. As a second step, the resulting icequake waveform time series was visually scrutinized, and deep icequakes were tracked considering characteristics mentioned by Walter and others (2008) and Dalban Canassy and others (2012). This analysis identified two deep event templates referred to as Ev1 and Ev2, with impulsive p- and s-wave onsets, reduced Rayleigh waves, duration below 0.1 s (Fig. 4a) and focal depths near the bedrock. When the 13 sensors array was not operational, Ev1 was visible on RA11 and RA13, whereas Ev2 only appeared in RA11 data. Figures 4a and b respectively, show the unfiltered vertical component of Ev1 recorded at RA13 and the associated power spectral density (PSD). PSD is found in the 8–50 Hz interval, while a clear corner frequency not visible in the spectrum of noise portions (Figs 5c and d) can be observed within the 50–70 Hz frequency band.

Fig. 4. (a) Near bedrock icequake clusters resulting from the cross correlation analysis using Ev1 (cl1–cl6) and Ev2 (cl7–cl20). For each cluster individual waveforms are plotted in gray and thick red (R1) or green (R2) depicting a stacks of all waveforms. Note that deep events in cl3 were detected using templates recorded on RA13 before 11 September 2012 and on RA11 for the following period. cl7–cl20 were defined using RA11 only. Events of cl17–cl20, as well as icequakes of cl3 detected after 3 October 2012 were sampled at 250 Hz. (b) Aerial and cross sectional views of hypocenter locations of cl1, cl3 (red pdfs) and cl8–cl14 (green pdfs). Layouts are the same as in in Figure 3. Cross sections positions are depicted by dashed black line in (b). (c) Normalized vertical components of basal events of cl1 and cl10 recorded on multiple seismometers with layout similar to (a).

Fig. 5. (a, c) Unfiltered vertical component of Ev1 and of 1 s noise portions recorded at RA13. (b, d) Associated spectrums with logarithmic scales. A corner frequency can be discerned in the spectrum of the basal events between 50 and 70 Hz.

Ev1 and Ev2 were subsequently used as templates to search for similar events in the continuous data recorded over the entire study period. This was achieved by means of a two steps cross correlation procedure based on the method developed by Helmstetter and others (2015): (1) each component of the high pass filtered (20 Hz) seismic signal of deep events templates and continuous seismic data was cross correlated with each other. For Ev1, data for RA13 and RA11 from 10 June to 11 September 2012 and 12 September 2012 to 24 July 2013 were, respectively, used. For Ev2, data for RA11 from 10 June 2012 to 24 July 2013 were used. Moreover, a decrease of sampling rate from 500 to 250 Hz between 2 October 2012 and 10 July 2013 combined with a sensor re-installation performed on 2 October made necessary to find analogous 250 Hz templates for both Ev1 and Ev2 for pursuing the analysis during this period. We retained this approach because seismogram downsampling did not provide satisfactory results. Finally, 3 s window of seismic signals were retained when associated cross correlation coefficient exceeded 0.5; (2) the vertical component of events in the resulting time series was thereafter cross correlated, and occurrences sharing cross correlation coefficient equal to 0.8 assumed to originate from the same source. Finally, clusters including icequakes recorded between 28 August and 2 October 2012 were tentatively located by picking p- and s-wave arrival times on stacked vertical components of the aligned waveform.

Results of the cross correlation procedure are presented in Figures 6 and 4a. Events similar to Ev1 and Ev2 were, respectively, found to be distributed into 6 (cl1–cl6) and 14 (cl7–cl20) clusters, including between 9 and 193 events for a total of 292 and 682 basal icequakes. For each template, the different clusters exhibit highly similar patterns with very close p- to s-wave arrival time intervals, indicating their proximity in space. Following the change of template required by the sampling rate decrease from 500 to 250 Hz at the end of September 2012, new clusters were defined (cl17–cl20). Although the hypothesis of new sources initiated exactly at this time cannot be completely ruled out, we reasonably assume here that activity periods of clusters found before and after sampling rate alteration were overlapping. In the case of Ev1, cl3 is found to unambiguously include events recorded with 500 Hz on RA13 and RA11 before 2 October 2012 and events recorded on RA11 with 250 Hz thereafter (respectively upper, middle and bottom frames labeled cl3 in Fig. 4a). In the case of Ev2, the low amount of samples comprising the 250 Hz signal in cl17–cl19 makes the comparison with cl7–cl16 uncertain. Despite the close similarity shared by waveforms of these clusters, we prefer in this context to separately analyze clusters defined after decrease in sampling rate and we are aware that our determined cluster activity represents a lower limit. Our results show a wide range of emissions duration characterizing the clusters (Fig. 6), from sub-hourly time spans (cl4, cl5) to several days/weeks (cl1–cl2, cl6–cl16) or inter seasonal (cl3, cl17–cl19) activity period up to 300 d.

Fig. 6. Clusters characteristics for both region R1 (solid red) and R2 (solid green). Up or down black arrows in the column entitled ‘polarity’ depict single first motion polarity of associated clusters. Dash indicate cluster, where no reliable determination was possible. Vertical black lines stand for icequakes occurrences, while technical data gaps are indicated by colored patches. Dashed red lines point out the times when sampling rate was altered.

Reliable locations could be derived for cl1, cl3 and cl7–cl13 (Fig. 4b). Hypocenters associated with cl1 and cl3 show overlapping pdfs covering, at the surface, an area of ~1500 m² eastward from RA10 with a maximum likelihood depth of ~85 m. In a similar way, clusters cl8–cl14 exhibit grouped hypocenters with overlapping pdfs covering a surface of ~4500 m² upstream from RA11 and focal depth ~80 m. In both cases, reliable constraints of focal depth, of ~10–20% of the ice thickness, allow us to unambiguously associate detected icequakes with subglacial icequakes released near the glacier bed. However, overlapping epicenter pdfs do not allow us to distinguish all the different clusters by means of location. In this regard we chose to define two regions, R1 and R2, respectively, delimited by epicentral pdfs, cl1 and cl3, and cl7–cl13 (Fig. 4b).

Reliable locations for other clusters were prevented due to an insufficient number of recording seismometers or poor quality in time arrival pickings. However, we argue that they most probably locate R1 and R2, first because of the high coherence threshold (0.8) used in the cross correlation procedure, and second due to the remarkable similarity of the time interval measured between p- and s-wave arrivals observed for located and unlocated clusters in both region. Considering this, we tentatively conclude that cl2 and cl4–cl6 locate within R1 and cl14–cl20 locates within R2. Respective icequake time series are consequently included in the two events lists illustrated in Figure 7b.

Fig. 7. Time series of the seismic, geodetic and climatic parameters studied in the investigated area from 10 June 2012 to 25 July 2013. (a) End-to-end overview of all computed daily spectrogramms for RA11, RA12 and RA13. Horizontal black lines indicate data gaps. RA11 and RA13 were respectively used for detection of events released from active regions R1/R2 and R1. Red circle denotes some of the temporary releases discussed in Section 6.1. (b) Inverse waiting time (iwt) for deep events time series released from R1 (red crosses) and R2 (green crosses) with semi log scale. Colored patches refer to periods with technical data gaps preventing the detection of icequakes from R1 (yellow), R2 (blue) and both R1 and R2 (magenta). Dashed red vertical lines point out the times when sampling rate was altered. (c) Horizontal surface velocity at stakes s23 (green) and s33 (red); (d) Hourly evolution of ice and snow ablation (red) and snow accumulation (blue) computed at s23. Yellow curve depicts ablation time series with 2 d smoothed values. Climatic data for the investigated area were obtained from Meteo Swiss measurements performed at Grimselpass.

R1 locates near the central flowline of the glacier where underlying bedrock slopes increase in the north-south direction from 7° to 13°; on the other hand R2 locates ~100 m from the lateral moraine where bedrock slope ranges from 8° to 20° in the west-east direction (Fig. 1a). We investigated the first motion polarities of our clusters by scrutinizing the stacked waveforms derived from each recording sensors (Figs 6 and 3a for stations RA11 and RA13 and Fig. 4c for example of clusters detected on multiple stations). The analysis indicates an unambiguous dilatational first motion, i.e. down, observed at all recording up- and down-glacier stations for cl1–cl3 and cl6 in R1. A single polarity is also confirmed on all up- and down-glacier recording sensors for some of the clusters located in R2, with a compressionnal up motion for cl7–cl16, and a down polarity for cl20 (Fig. 6). In this case the determination of the first motion polarity was sometimes made uneasy due to poor signal to noise ratio, as shown for the traces of events released from cl10 and recorded on RA04 (Fig. 4c).

5. DYNAMIC AND SEISMIC SEASONAL VARIATIONS

5.1. Changes in glacier dynamics

Horizontal surface displacements are shown in Figure 1b. Between 19 June and 2 October 2012 (black arrows in Fig. 1b), total horizontal displacements measured at stakes s33 and s23 were 4.53 and 9.20 m, respectively, with daily mean values of 0.068 and 0.082 m d⁻¹. Between 19 June 2012 and 4 July 2013 (purple arrow), displacement at s23 amounts up to 25 m with a mean velocity value, 0.074 m d⁻¹. Surface motions derived from boreholes sensor positions from 6 May to 2 October 2012 (red arrows) shows values equal to 4.50, 13.95 and 7.83 m at RA11, RA12 and RA13, respectively. A comparable range of velocity values was found at both stakes, from a few cm up to 0.35 m d⁻¹ (Fig. 7c). Associated flowlines were found oriented north-west/south-east at s23 and RA12 and north-east/south-west at s33 and RA13, while RA11 followed an intermediate direction, north/south. Velocity differences observed between upstream (RA12 and s23) and downstream (RA13 and s23) measurements reveal a significant longitudinal gradient in the study area, combined moreover with a strong change in surface motion orientations most probably between RA05 and RA10.

We identify two main regimes of surface velocity, largely aligned with melting and snow accumulation periods (Fig. 7d). Highest velocities are found before 29 November 2012 and after mid April 2013, together with enhanced melting conditions. During this period, referred to hereafter as the summer regime, position measurements at s23 show a mean daily velocity equal to 0.081 m d⁻¹ and a horizontal displacement of 15.75 m, which accounts for 64% of the motion recorded over the entire study period. Moreover, diurnal patterns in surface motion aligned with diurnal melting cycles (insets Figs 7c and d) are generally observed.

In contrast, low horizontal displacements with a mean daily velocity value (s23 position) of 0.060 m d⁻¹ prevailed for the rest of the investigated period. This regime, referred to hereafter as the winter velocity regime accounts for 36% of the total displacement. It is characterized by no or very low melting combined most of the time with a substantial snow cover. It begins on 29 October 2012 and ends, most probably, in the second half of April, when unfortunately, data are missing. Such an assumption is supported by two main arguments: first considering the recovery of snow melting on 15 April 2013 (Fig. 7d), and second given the increased velocity measured shortly thereafter on 28 April 2015. Moreover, such surface motion increase in the presence of snow cover (Spring event) has already been documented (e.g. Kavanaugh and Clarke, 2001; Anderson and others, 2004). The authors pointed out a 1 month time delay between first rise of borehole pressure and rise of velocity, whereas this delay is of 10 d maximum in our study if we take the data gaps into account. Winter velocity regime does not exhibit any particular variations at diurnal or longer time scales. The increase in mean daily velocities from winter to summer regimes is ~25%, comparable with what Harbor and others (1997) measured on Haut Glacier d'Arolla (Switzerland), but remains smaller than observations of Macgregor and others (2005) on Bench Glacier (Alaska).

The evolution towards a winter or summer regime (29 October 2012, mid April 2013) shows very rapid transition phases; ~1 d at the start of winter and 10 d at the start of the summer regime. These two moments closely coincide with changes in melting rate (Fig. 7d). Indeed, winter regime initiation happens 2 d after a drop of melting on 27 November 2012 and ends no more than 10 d after it suddenly recovers on 13 April 2013. Such seasonal correlation, together with the high correlation between diurnal variations of surface motion and melting rate, confirms the control of liquid water supply on surface motion in our study area.

5.2. Evolution of seismicity

5.2.1. Changes in glacier seismicity

An end to end overview of computed spectrograms is shown in Figure 7a. Within the entire 1–115 Hz frequency band, there exists a general agreement among the temporal distribution and intensity of spectral density of the 3 seismometers. This consistency is for instance, illustrated by the strong increase of spectral energy noticed on the 3 sensors on 13 June 2012, as well as by the coherency observed between RA11 and RA12 from October 2012 to June 2013. It therefore shows that data recorded by a single sensor did not consist of site specific noise but accounts for the large scale seismicity (Roeoesli and others, 2014).

Most of the seismic energy recorded in the study area was released before 27 October 2012 and after 15 April 2013, in agreement with periods of summer velocity regime. Highest spectral density values up to 40 dB and broadest emission spectra were however found before 27 October and after 11 June. Two predominant frequency intervals are roughly evidenced during this period: a quasi continuous and largely predominant seismic release ~25–30 dB recorded on the three sensors between 1–15 Hz and 1–10 Hz for RA11, RA13 and RA12, and less energetic emissions characterized by density values ~15 dB in a range up to 80 Hz.

On the contrary, few and distinct seismic energy releases were recorded in the time interval between 27 October 2012 and 15 April 2013, coinciding with the winter velocity regime. The consistent 1–15 Hz emissions now constitutes the main signal, however with reduced density values ~10 dB. Other temporary releases lasting no more than a couple of days are also reported on very narrow frequency bands between 40 and 60 Hz and ~100 Hz with density values ~5 dB, mostly from 27 November to 31 March (red circle in Fig. 7a).

5.2.2. Changes in basal seismicity

Besides data gaps, changes in ambient background noise level are likely to affect deep icequake detection sensitivity and consequently the completeness of the R1 and R2 time series. We propose two arguments indicating that this may be neglected: First we underline that our high performance cross correlation method allows the detection of events with amplitudes below the noise level (Schaff, 2008) and second, our results point out a high activity from R2 throughout October together with changing melting rates (0–4 mm w.e h⁻¹; Fig. 7e) and therefore significant background noise variations (Dalban Canassy and others, 2011). In this context, we say that both deep events times series can reasonably be assumed not significantly affected by noise changes within intervals with continuous data availability (Fig. 7b).

Basal icequakes within R1 were detected from 19 July 2012 to 4 July 2013, mostly as sporadic bursts with duration ranging from a few hours to 1–2 weeks (Figs 6 and 7b). Taking into account the data gap from 23 December 2012 to 23 March 2013, this region is found to be inactive during at least 6.5 months between 18 October 2012 and 30 June 2013, when associated sources are suddenly reactivated.

In contrast, R2 is associated with a more continuous activity. The icequakes dataset ranges from 10 September 2012 to the end of the investigated period. It is characterized by high activity between 22 September 2012 and the beginning of the winter velocity regime, with activity bursts from 22 to 26 September 2012. A rapid decrease in activity then follows between 1 and 12 November 2012 and a stabilization until the 21 December 2012 is however interspersed by an intense 2 d peak of seismicity on 9–10 December 2012. Seismicity restarts on 10 April 2013, after an interruption of at least 18 and 107 d or at more considering data gap. The activity progressively increases until the end of the study period with a peak on 3 July 2013.

The data gaps makes an analysis of seasonal variations of basal seismicity difficult. However, we note that activity in R2 tends to follow the seasonal trend of glacier seismicity, mimicing transition periods between summer and winter velocity regimes. Variation in R1 activity does not exhibit a particular trend, except that associated icequakes were exclusively detected during the summer regime. We also notice that periods of enhanced velocity are not necessarily accompanied by enhanced basal seismicity.

6. DISCUSSION

Over the 410 d monitoring period on the tongue of Rhonegletscher, our study highlights the following results:

• • the surface motion in the study area exhibits two distinct velocity regimes closely correlated with changes in liquid water supply. The sudden variations at the beginning and end of the winter velocity regime suggests a rapid adaptation of the basal water drainage system;

• • 20 clusters of basal icequakes were found within two active regions. They are characterized by activity duration ranging from a few hours to more than 300 d, and by waveforms with single polarity compressive or dilatational first motions;

• • the fluctuations in surface velocity correlate with the distribution of glacier seismicity over the study period as well as with variations in seismic activity at R2;

• • near bedrock seismicity from R2 persists during a period of at least 2 months after the beginning of the winter velocity regime;

• • within the 1–15 Hz frequency band a quasi continuous seismic energy release is recorded during the whole hydrological year.

In the following, we interpret the variations in both the spectral energy and the activity of basal icequakes in the light of changes in water input and surface velocity in the study area.

6.1. Glacier seismicity

The seasonal correlation between changes in teglacier velocity regime and spectral power of our continuous seismic recordings suggests a predominant control of ice dynamics on the glacier seismicity recorded in the study area (Figs 7a and d). This relationship had also been observed in previous studies over shorter periods (Neave and Savage, 1970; Roux and others, 2008; Mikesell and others, 2012; Roeoesli and others, 2014). We adopt the view of these previous studies that melt-driven ice flow acceleration promotes tensile fracturing and thus englacial seismicity, and we emphasize that this relation is also reflected at a seasonal time scale. At the same time, the seismic signature of englacial tensile fracture events covers frequencies between 10 and 100 Hz (Faillettaz and others, 2008; Roux and others, 2008). Consequently, fracture events alone cannot explain the persistently recorded seismic energy between 0 and 15 Hz (Fig. 7a).

In order to investigate the possible sources, we focus on temporal changes in spectral intensity within the 1–15 Hz band at the 3 borehole sensors together with smoothed modeled surface melting rate (Fig. 8). During the summer velocity regime, periods of enhanced surface melt coincide with elevated seismic energy. This can be explained by water-flow induced seismic tremor (Lawrence and Quamar, 1979; Métaxian and others, 2003; Roeoesli and others, 2014) in englacial conduits and fractures. On the contrary, spectral peaks are suppressed during the winter velocity regime. At this time, cyclic noise bursts below 5 Hz between 0600 and 1800 h during weekdays dominate the continuous signal. Considering that three hydro-powerplants are located with a radius of 10 km, we suggest these <5 Hz vibrations are induced by turbines. Monochromatic spectral peaks, such at ~12 Hz are likely related to electronic noise within the seismic station.

Fig. 8. Yearly evolution of the 1–15 frequency band for RA11, RA12 and RA13 and variations in melt water supply in the study area smoothed with a 24 h sliding window (solid red). Red rectangles indicate the monochromatic releases ~12 Hz as discussed in Section 6.1.

Another observation deserving comment is the very temporary and narrow band signal ~40–60 and 100 Hz during the winter velocity regime (red circles in Fig. 7a). As they were mainly recorded concurrently with enhanced solid precipitation on a growing or already thick snow cover (Fig. 7d), we propose fresh snow transport at the snow surface due to wind effect as the most likely source.

6.2. Interpretation of near bedrock seismicity variations

Two main source mechanisms are currently considered to explain basal seismicity. First the opening and closing of cracks in the bedrock vicinity (Walter and others, 2009; Dalban Canassy and others, 2013) and second the release of seismic energy due to shear faulting during stick-slip motions, first observed in Antarctica (Anandakrishnan and Alley, 1997; Wiens and others, 2008; Winberry and others, 2009; Zoet and others, 2012) and more recently in the Alps (Helmstetter and others, 2015).

We interpret the single first motion polarity of 15 of the 20 clusters located in R1 and R2 as evidence for strong isotropic components of the underlying moment tensors (Walter and others, 2013). This suggests that tensile fracturing is the most likely source model, with fracture opening or closing corresponding respectively to up and down polarity (Walter and others, 2008, 2009; Dalban Canassy and others, 2013). We do acknowledge that due to limited coverage, the possibility for shear motion cannot be fully excluded. However, we remark that the surface velocities (Fig. 1b), indicate potential stick-slip in the study area would preferentially occur in the north-south to north east-south west direction. If we assume our basal events had double-couple mechanisms, the following first motion polarities can be expected: up at RA13 and R10 and down at R04 and R05 for clusters included in R1, up on RA11 and RA10 and down at RA06 for clusters included in R2. Looking at the traces of events belonging to cl1 (R1) and cl10 (R2) and detected on multiple stations located up- and down-glacier (Fig. 4c), we observe (1) a clear negative onset at RA13 and RA10 for cl1 and therefore an indisputable down motion single polarity for R1, and (2) a positive onset at RA06, RA01 and RA11 and for cl10, which similarly points out a single up polarity for R2. In this context, we emphasize the lack of any mixed polarity observations and the single polarity observed at up- and down-glacier seismometers for both active regions, concluding against significant shear faulting. Additional information on the sources dimensions could possibly be derived from the rate of decrease of the PSD following the corner frequency highlighted for Ev1 (Fig. 5; Section 4.3). However, we remark that the standard approach assuming a shear fracture occurring within a simple fault geometry (Stein and Wysession, 2009) needs to be adapted to tensile crack icequakes before application to our case. This is planned in further investigations.

Several previous studies have focused on activity variations in basal fracturing during high melt periods. Walter and others (2008), Walter and others (2009) and Dalban Canassy and others (2013) found that tensile fracturing occurs primarily during low or decreasing basal water pressures. This led Walter and others (2013) to suggest that such basal fracturing is a result of growth and shrinkage of subglacial cavities during large pressure fluctuations in adjacent R-channels. However, in our present study we observe basal fracturing during the winter season when melt-induced fluctuations in basal water pressure are absent. Moreover, sources located in cavity walls should be affected by collapse and shrinkage during winter, which contradicts the >280 d activity periods spanning over winter 2012, observed for some clusters. This suggests that our events are not related to the hydraulics of subglacial cavities.

Here we propose that our basal events are independent of subglacial cavity evolution and are a manifestation of a dynamic fracture network near the glacier base (Van der Veen, 1998). The role of such a network has been acknowledged for temperate glaciers by Fountain and others (2005), Hock and others (1999) and Stuart and others (2005) and for the Greenland ice sheet by Das and others (2008) and Andrews and others (2014). We interpret the coexistence of different clusters active in both regions R1 and R2 as neighbouring cracks presenting a very close configuration, which explains the high similarity shared by the associated waveforms. Furthermore, our results indicate that fracture networks evolve and consequently exist over the entire year, which suggests non-negligible water movement even for no or very reduced surface melting. This is in line with previous findings by Fountain and others (2005). Consequently, fracture networks developing during the winter relaxation phase of basal drainage system may serve as a starting point for R-channel formation as previously suggested (Walder, 1982; Fountain and others, 2005). While the origin of liquid water during the winter period remains unclear, we surmise that such a supply depends on the glacier characteristics and on local conditions, such as moraines or glacier lakes. In our case, we propose water storage and strain heating as potential water sources during the winter season. Additional inputs may also be released from snow melting on the lateral moraines and infiltrating the basal ice layers.

We remark that our basal fractures were found near moderately steep bedrock (10^o–20^o), in places characterized by pronounced velocity gradients favoring large strain rates (Figs 1a and b). These two particularities may thus constitute some of the characteristics of birthplaces for near bedrock fractures. We investigated the possibility that fracture networks are advected under the influence of glacier flow by analyzing potential changes in distances between icequake sources located in R2 and station RA11. To do so we measured p- and s-wave arrival times at RA11 belonging to the events of cl17 and cl18. We then assumed that a change in p–s delay times from the first to the last event within these clusters is due to a source displacement relative to station RA11. Over our measurement period, the arrival time lags of p- and s-wave changed by no >4 ms, corresponding to only two sampling intervals. Whereas such a change corresponds to a source displacement of ~15 m, it also lies within the uncertainties of our arrival time picks. With our available data we can therefore make no assertions about temporal changes in cluster locations.

7. CONCLUSION

Our year-long monitoring on Rhonegletscher's tongue revealed seasonal changes in the glacier and near bedrock seismicity with elevated activity during the summer velocity regime. While correlation between fracturing activity and surface motion has verified a seasonal time scale, our results also point out the non-negligible role played by sources related to flowing water and non-glacial mechanisms in the content of the yearly spectrogram.

Perhaps most importantly, our data reveals icequake sources, which were active for an entire year between summer 2012 and summer 2013. These icequakes are most likely a manifestation of hydrofracturing and their occurrence during winter months suggests non-negligible water movement even in the absence of surface melt. Based on this, we also suggest that such crack-networks may serve as a starting point for R-channel formation. To our knowledge, this is the first documentation of glacier and basal seismicity fluctuations in the ablation zone of an Alpine glacier during a full hydrological year.

Our data was insufficient to study advection of fracture sources in glacier flow. This matter could be investigated in the future with seismic monitoring programs lasting several years. Moreover, spectral analysis of icequake signals, including methods based on corner frequency, and sustained tremors could provide further insights into the volume and extent of subglacial fracture networks and their capability to transport water and to interact with channels of more efficient subglacial drainage systems. These investigations might be combined with geophysical exploration of the subglacial water drainage system by means of ground-penetrating radar measurements, as shown for instance for a low thickness glacier in the work of Gabbud and others (2015).