3.1 Discharge

Over the course of the 2016 melt season, Lac des Faverges steadily filled (orange dashed line in Fig. 2a) and reached a maximum volume of approximately 2 million cubic meters of water at the end of August. In the evening of 27 August (day 240), a lake drainage through the moulin marked in Fig. 1a was initiated and emptied the lake basin in approximately 6 d. The moulin routed the water to the subglacial environment, and it escaped the glacier beneath the Rätzligletscher on the northern side. In the first hours of the drainage, the water escaped abruptly, since the moulin reached the bottom of the lake. This drainage phase corresponds to the peak in the discharge curve (≈25 m³ s⁻¹) measured in the Simme Valley (blue curve in Fig. 2a and 2b). This peak discharge overwhelmed the capacities of the subglacial drainage system, which is indicated by the local ice uplift measured at all three GPS stations (Fig. 2b). As the lake level fell to the elevation of the moulin inlet, the direct connection between moulin and lake became disrupted. Subsequently, the lake connected to the moulin through a supraglacial channel which steadily incised deeper into the ice but slowed down the drainage (6–11 m³ s⁻¹). The exact transition time to this state is unknown but was within the first day of the drainage initiation.

Figure 2(a) Hydrological data recorded in the vicinity of Glacier de la Plaine Morte: discharge of the Simme River measured in the Simme Valley (blue curve; ≈4 km line of sight from the glacier terminus), level of the Trüebbach stream (gray curve; ≈1.5 km from glacier terminus), height of the water column in the lake above the pressure sensor (orange dashed), and precipitation at stations ABO and MVE (blue and purple bars; 12 km north and 10 km south of the glacier, respectively). (b) Discharge and lake level for the drainage period (same as a) and the vertical displacement of three GPS units (black lines). (c) Spectrogram of station PM05 for the same time period shown in (a). (d) Zoom of the spectrogram in (c) showing anthropogenic noise. (e) Zoom of a spectrogram of station PM32 showing moulin resonances. The white bar indicates a data gap due to station maintenance. (f) Zoom of the spectrogram in (c) showing the lake drainage. Note that data used to calculate the spectrograms are not corrected for the instruments' phase responses.

Download

Discharge magnitudes similar to those of the lake drainage period were also measured in the Simme River prior to the lake drainage (three peaks on days 213–225) and after the lake drainage (days 248–252). Most of these discharge peaks can be linked to rainfall events having a shorter duration than the lake drainage (precipitation data are provided by the Switzerland's Federal Office of Meteorology and Climatology MeteoSwiss). Since precipitation affects the entire catchment above the gauging station (more than 4 times the glacier surface area), these precipitation-related discharge events need to be interpreted with caution because part of the measured discharge at those times may be due to water flowing outside of the glacier. In general, however, the similarity of the discharge curve and the stream level height measured close to the glacier terminus suggests that Glacier de la Plaine Morte is the main contributor to the discharge measurements. In addition to the drainage of Lac des Faverges, a smaller supraglacial lake at the southern rim of the glacier (labeled SL in Fig. 1b) was observed to drain via a supraglacial canyon routing water to moulins. A field visit on day 239 revealed that the lake was draining, but the time of the drainage initiation was not witnessed.

3.2 Seismic tremors

Figure 2c shows a spectrogram for station PM05. Recent studies suggest that water routing in subglacial conduits generates seismic tremors observable in the frequency range 1–10 Hz (Bartholomaus et al., 2015; Gimbert et al., 2016). In this frequency range, however, we observe several signals of anthropogenic origin. These include a diurnal signal from Monday to Friday with sharp onset and decay times, a monochromatic signal visible as a spectral line at roughly 2 Hz starting from day 156, and most likely also the diffuse band centered around 5 Hz (Fig. 2c–d). Regarding glacier seismicity, we identify a harmonic moulin tremor with three prominent frequencies which indicate resonant modes in the water column, similar to those in Roeoesli et al. (2016) (Fig. 2d). During the lake drainage, the signal strength is increased for frequencies greater than 1 Hz, and we observe high-frequency tremors (>3 Hz) during the drainage initiation (Fig. 2c and e).

To better distinguish the seismic signal contributions, we investigate the wave field in more detail. For this purpose, we calculate 3-D particle motion polarization attributes following Koper and Hawley (2010). This approach is based on an eigen-decomposition of the spectral covariance matrix containing the power and cross spectra of a single three-component station (Vidale, 1986). One of the polarization attributes, the difference in phase between the vertical component and the principal horizontal component, ϕ_VH, allows us to distinguish between different wave types. In particular, the elliptical particle motion of a Rayleigh wave is caused by a 90^∘ phase shift between vertical and horizontal ground motion and distinguishes it from other wave types. To calculate the polarization attributes, we use the freely available toolbox hosted on the IRIS web page (http://ds.iris.edu/ds/products/noise-toolkit/, last access: November 2019) with the default parameters and processing steps (including instrument response removal; see Koper and Hawley, 2010). Figure 3 shows probability density functions of ϕ_VH for a station of each array. Consistent with an elliptical particle motion in the vertical–radial plane associated with Rayleigh wave propagation, ϕ_VH clusters around ±90^∘ in the frequency range 8.5–12 Hz for all four stations shown. Below 8.5 Hz (6 Hz for station PM33), i.e., frequencies where anthropogenic noise is evident, clustering around ±90^∘ indicative of Rayleigh waves is not present or only in narrow frequency bands (e.g., 4–5 Hz for station PM05). We do not find a difference in the polarization results prior to and during the lake drainage, though the short duration of the drainage process hinders a detailed comparison by means of a statistical representation as shown in Fig. 3.

Figure 3Probability density functions for the phase difference between the vertical component and principal horizontal component. Probabilities are calculated for the time period 22 July 2016 to 6 September 2016 (23 August 2016 for station PM33) and bins of 5^∘ width. The results are shown for one station of each array: (a) PM05, (b) PM11, (c) PM21, and (d) PM33.

Download

Since the continuous recordings below ≈8.5 Hz are contaminated by anthropogenic noise with a complex wave-type signature, we chose to analyze the frequency range 8.5–12 Hz in the context of glacial hydraulics. This frequency range is dominated by Rayleigh wave energy, which facilitates the tremor location analysis in the next section. To investigate possible correlations between discharge and seismicity for the 8.5–12 Hz range, we calculate the tremor amplitude, or median absolute ground velocity, for the vertical component of ground velocity as described in Bartholomaus et al. (2015). Figure 4 shows the resulting time series for a station of each array along with the discharge recordings. Prior to the lake drainage, variations in tremor amplitude are weak but follow the discharge curve (e.g., days 213–225). In the 4 d preceding the lake drainage, the daily melt cycle due to high temperatures is visible in both the discharge time series and the tremor amplitude curves. During and after the lake drainage, tremor amplitudes are increased and show stronger variations than prior to the lake drainage. From Fig. 4b we note that PM21's tremor amplitude correlates with discharge particularly well.

Figure 4(a) Tremor amplitude (8.5–12 Hz) time series for a station of each array (thin colored lines) and discharge (thick blue line). (b) Zoom of the gray shaded area in (a).

Download

3.3 Icequake activity

To investigate the interplay of glacial hydraulics and icequake activity on Glacier de la Plaine Morte, we build on the results of Lindner et al. (2019). This study focuses on ice-fracturing events (Walter et al., 2009) to investigate azimuthal anisotropy of seismic wave propagation, but the events are also useful to study fracturing associated with glacial hydraulics, e.g., during outburst floods (e.g., Roux et al., 2010). Lindner et al. (2019) detect icequakes by applying a short-term average and long-term average (STA–LTA) trigger (Allen, 1978) on bandpass-filtered data (10–20 Hz for A1, A2, A3; 7–15 Hz for A0), require at least three stations of an array to trigger concurrently, and disable the trigger for 3 s after a detected event to avoid overlapping event windows (for parameter details see Lindner et al., 2019).

From the event catalogs from each array, we calculate the icequake detection rate in events per hour. Figure 5 shows that icequake activity is often increased during discharge peaks, though not always. Given the correlation between tremor amplitude and discharge (Fig. 4), this implies that the tremor amplitude in turn is also correlated with the icequake rate (see blue arrows in Fig. 5). In addition, we identify times when correlation of the tremor amplitude with the icequake rate is higher than with discharge (red arrows in Fig. 5). These features correspond to icequake bursts lasting on the order of hours but less than a day. Maximum detection rates are 352, 314, 172, and 20 icequakes per hour for A0, A1, A2, and A3, respectively. For A0, this corresponds to 5.87 icequakes per minute and thus almost 18 s of disabled trigger per minute (trigger disabled for 3 s after event). This suggests that our results are a conservative estimate of icequake occurrence. We note that especially those arrays with high icequake rates (A0 and A1) suffer from icequake contaminated tremor amplitudes. This contamination might be reduced by choosing different window lengths for the calculation of the tremor amplitude, but investigating this matter is beyond the scope of this study.

Figure 5Icequake detections per hour for all arrays (gray bars; see text for details) and the discharge curve of the Simme River (blue line). The black, magenta, orange, and green lines (from top to bottom) are the tremor amplitude curves shown in Fig. 4. Note the different tremor amplitude scaling between the two panels. Blue arrows indicate times where the tremor amplitude correlates with both discharge and icequake rate. Red arrows indicate times where the tremor amplitude correlates with icequake rate only. The black dashed rectangle indicates times, where three of five A2 stations tipped over due to diminishing snow cover. The icequake rates in this interval need to be taken with caution.

Download

For further insights into the glacial hydraulics, we consider the icequake source locations as determined from plane-wave beamforming (for details see Lindner et al., 2019). Plane-wave beamforming enables us to determine the back azimuth of an incident signal with an array of sensors. Its generalization to epicentral coordinates will be introduced in Sect. 4.

Figure 6b shows the icequake detections of A1 as a function of focal time and source back azimuth. Some peaks in discharge, e.g., on day 222, are accompanied by fracturing at various source back azimuths. Since A1 is in the glacier's center and icequakes arrive from various back azimuths, this suggests that these discharge events affect large portions of the glacier. Other events, e.g., the melt cycle in the days prior to the lake drainage, are accompanied by more localized seismicity at back azimuths of 50–100^∘ only. The latter is also the case for the ≈24 h preceding the drainage initiation, where the seismicity at the back azimuth towards the main drainage moulin is increased. We also detect icequakes at this back azimuth earlier, but activity is not sustained and back azimuths do not focus on the moulin. With the onset of the lake drainage, fracturing occurs at various back azimuths with a focus on the lake basin. After the drainage, fracturing is predominantly confined to the lake basin as well.

Figure 6(a) Discharge and lake level (same as in Fig. 2). The vertical red dashed line indicates the drainage initiation. (b) Detected icequakes at A1 as a function of time and source back azimuth (white dots on black background). Icequake clustering in both time and back azimuth is visible as bright white spots. The two horizontal red dashed lines indicate the back azimuth from the array center to the main drainage moulin ±5^∘. (c) Icequakes per hour in the back-azimuth range marked with the two horizontal red dashed lines in (b).

Download

We note that STA–LTA detection thresholds might be affected by changes in the background noise level (Walter et al., 2008), resulting in biased event detections. However, since our focus is on periods with high discharge or strong melting (as in the hours prior to the drainage initiation) in which trigger sensitivity is typically decreased, we argue that our results are robust.

4.1 Matched-field processing

To locate the sustained tremor sources, we apply matched-field processing (MFP; Baggeroer et al., 1993) to our four seismic arrays. MFP measures signal coherence of a phase across an array of receivers and matches it against a synthetic wave field computed for a point source and a velocity model. By testing various source positions and velocity models for the synthetic field, a grid search finds the combination of source position and velocity model which best matches the measured coherence across the array. The result is the most likely source location and velocity model. Allowing near-field point sources, and hence circular wave fronts, MFP is a generalization of the conventional plane-wave beamforming approach used to determine the slowness and back azimuth of incoming waves (Rost and Thomas, 2002). In case the distance of a source to the array is greater than 2 to 3 times the array aperture, the circular wave front approach converges towards a plane-wave solution (Almendros et al., 1999). For MFP, this implies that far-field sources allow a back-azimuth estimate only (as is the case for plane-wave beamforming), while near-field sources can be associated with epicentral coordinates. We leverage this to locate tremor sources, some of which are in the arrays' near field as we show in the following.

The workflow for MFP is as follows (for a more detailed introduction, see, e.g., Corciulo et al., 2012). From the time-domain ground-velocity recordings of N receivers grouped to an array, the discrete Fourier transforms at some angular frequency of interest ω are calculated. The resulting complex frequency-domain values are arranged to form a column vector d(ω) of length N. From this column vector, the cross-spectral density matrix K(ω) is calculated as

where † denotes the complex conjugate transpose operation. Note that the diagonal elements of K(ω) are the autocorrelation values of the N receivers at ω, while the off-diagonal elements are cross-correlation values of receiver pairs. Both autocorrelation and cross-correlation are discrete values associated with angular frequency ω, and the latter represent average phase delays between two receivers at ω. The synthetic field at ω is calculated for each of the j=1, …, N receivers as

where i is the imaginary unit, r_j is the source–receiver distance of the jth receiver, and c is the phase velocity of the velocity model, which is constant in the case of a homogeneous ice body. Note that this representation focuses on phase information and disregards amplitude information. For j=1, …, N, the complex ${\stackrel{\mathrm{̃}}{d}}_j$ values are arranged to form the synthetic column vector $\stackrel{\mathrm{̃}}{\mathit{d}}\left(\mathit{\omega }\right)$ (equivalent to the data vector d(ω)), and phase matching is achieved via the conventional Bartlett processor (Baggeroer et al., 1993),

or via a high-resolution MFP method, i.e., the minimum variance distortionless response (MVDR) beamformer (Capon, 1969; Corciulo et al., 2012),

where K⁻¹ is the inverse of the cross-spectral density matrix K. In case the incoherent noise power is small relative to the power of the signal of interest, the MVDR processor is capable of estimating the source location and the velocity beneath the array with higher resolution than the Bartlett processor (Capon, 1969). Note that there is a trade-off between high-resolution and robustness; i.e., in contrast to the MVDR processor, the Bartlett processor might still produce meaningful results if the incoherent noise power is increased.

4.2 Single-array results

To investigate the spatial variability of tremor sources across Glacier de la Plaine Morte, we apply MFP to all four arrays individually. For this purpose, we use a sliding window of 15 min length (without overlap) over the entire data set to also resolve temporal variations. Each of these 15 min segments are processed as follows. To suppress incoherent noise, we calculate the ensemble average of K(ω) at discrete frequencies, using a 10 s sliding window with 50 % overlap (e.g., Corciulo et al., 2012). The overlap criterion yields a total set of 179 windows over which we average. In the frequency range of 8.5–12 Hz, the results from our polarization analysis (Fig. 3) suggest Rayleigh waves whose amplitude is correlated with discharge. For this reason, we calculate the MFP results in 0.1 Hz steps and average over the frequency range of 8.5–12 Hz. For the velocity c in Eq. (2), we use the local Rayleigh wave velocities which are 1600, 1800, 2200, and 1600 m s⁻¹ for arrays A0, A1, A2, and A3, respectively (Lindner et al., 2019). In the spatial domain, we apply a grid search over the entire glacier surface and its surroundings (assuming a horizontal plane) to calculate the r_j values in Eq. (2) with a spacing of 25 m in the x and y directions. Figure 7a shows the spatial clustering of the MFP results from all available time windows using the MVDR processor. The picked and shown epicenters are associated with the maximum B_MVDR value of all tested coordinates.

Figure 7(a) MFP locations (MVDR processor) over the frequency range 8.5–12 Hz assuming Rayleigh wave velocities (colored dots). Each dot is the result obtained for a 15 min window. The thick black line indicates the glacier margin in 2015, the black triangles the locations of the seismic stations, and the gray dashed circles the distance of twice the array aperture from the array center. The black “X” symbols indicate positions of moulins identified from orthophotographs (by swisstopo, SWISSIMAGE). The red “X” marks the position of the lake drainage moulin. Labels of type A0-2 refer to dominant source clusters discussed in the text. (b–e) Temporal variation in back azimuth and distance of the tremor source locations (colored dots) from (a) as seen from the array centers of A0, A1, A2, and A3. The gray line depicts the discharge curve measured in the Simme Valley for reference.

4.3 Multi-array results

To further constrain the tremor source locations, we stack the results obtained from the different arrays. Following the argumentation in the previous section, we continue to focus on Rayleigh waves and consider the MFP results obtained for the MVDR processor on the entire spatial domain used for the grid search. Figure 8 (left column) shows the results for a 15 min window on day 214 during a peak in discharge caused by precipitation. As reported in the previous section, A0 sees a persistent source to the north of the array, A1 and A2 point towards the glacier tongue, and A3 points toward the (south)west. For A3, however, a secondary lobe of the MVDR output is visible, which points to the glacier tongue as well. We combine the information from different arrays by stacking the MVDR grid-search results, which shows high MVDR values in regions where multiple arrays locate signals. The stacking allows triangulation and confirms that the main tremor source is in the region of the glacier tongue (Fig. 8, left column). We tested other time windows and found that the depicted situation is representative for the pre-drainage period which appears stable with little excursions to other source regions (see Fig. 7b–e). With the onset of the drainage, the tremor source locations change, as shown in Fig. 8 (right column). The depicted situation shows the result for a 15 min window on day 243 roughly 55 h after the drainage initiation. A0 now locates the tremor signal south of the array and A1 points towards the southeast with a secondary lobe pointing to the glacier tongue. A2 again points towards the glacier tongue with less scatter compared to Fig. 8, left column. As discussed in the context of Fig. 7b–c, A1 (secondary lobe) and A2 back azimuths are slightly increased compared to the pre-drainage period. Stacking the results from A0, A1, and A2 again (A3 has no data) shows two source regions, the glacier tongue and the main drainage moulin.

Figure 8MFP results obtained using the MVDR processor and Rayleigh wave velocities. Shown are the results for the single arrays (rows 1–4) and a stack of the arrays (last row). The left column shows the results of a 15 min window on day 214 during a peak in discharge caused by precipitation. The right column depicts the results of 15 min windows during the lake drainage on day 243. Exact times are given on top of the plots. The spacing of the ticks on the x and y axes is 500 m (see also Fig. 1).

5.1 Tremor composition

The results from our tremor analysis demonstrate that the recorded seismic wave field on timescales beyond those of discrete single events is generated by various processes. Apart from cryo-seismicity, we observe signals of anthropogenic origin. The diurnal signal occurring on working days only (Fig. 2c; also reported in Preiswerk and Walter, 2018), originates to the south of Glacier de la Plaine Morte (determined by plane-wave beamforming), likely in the Rhone valley where industry is located. The frequency range of anthropogenic noise (e.g., Anthony et al., 2015) often overlaps with the discharge–tremor band, meaning that glacioseismological data need to be analyzed carefully in glaciated regions with anthropogenic activity such as the European Alps to avoid misinterpretation. This also holds in the absence of anthropogenic noise, since our data reveal that tremors may be generated by different aspects of glacier hydraulics at the same time. We identify tremors which are dominated by energy released through ice fracturing (A0 and A1), are located at moulin locations (A0 and A3), or exhibit a characteristic frequency signature of moulin resonances (A3) and thus obscure turbulent-flow tremors. However, at A2, we argue that the recorded tremors are generated by subglacial water routing for the following reasons: (i) the tremor amplitude correlates with the discharge curve (Fig. 4) and (ii) MFP shows a persistent source in the region of the glacier tongue (Fig. 7a), from where (iii) the main glacier outlet emerges (Fig. 9). We note that subglacial water routing in turn can generate tremors both via pressure fluctuations in turbulent flow and via impact events from bed load sediment transport (Gimbert et al., 2014, 2016). Recent studies (Bartholomaus et al., 2015; Gimbert et al., 2016) typically separate the two processes by frequency at around 10 Hz; thus the frequency range associated with our results (8.5–12 Hz) may contain both processes. Even though we cannot exclude that bed load significantly contributes to total seismic power, we see evidence for water tremors being the dominant source for the following reasons. (i) The frequency ranges are controlled by various parameters (channel-to-station distance and channel apparent roughness among others) also permitting turbulent-flow tremors above 10 Hz (Gimbert et al., 2014). (ii) Ice flow of Glacier de la Plaine Morte is negligible (<1 cm d⁻¹ at A2, not shown), resulting in little sediment production by abrasion (Hallet, 1979), which we expect hinders bed load tremor generation. (iii) The A2 tremor–discharge scaling as discussed later tends to follow the drainage-regime predictions for water tremors without evidence for a hysteresis due to sediment flushing (e.g., Gimbert et al., 2016). Apart from various tremor sources, we finally note that the tremors are composed of different wave types, further increasing the complexity of the tremor signal.

Figure 9Upstream area distributions calculated from the hydraulic potential (see text for details). (a) Solution obtained for (spatially uniform) water pressures of half the ice overburden pressure. (b) Solution obtained for (spatially uniform) water pressures equaling the ice overburden pressure. The white triangles indicate the positions of the seismic stations.

5.2 Temporal evolution of the drainage system

Figure 7b–e shows that the tremor locations change over time. Since icequake tremors typically last on the order of hours (Fig. 5) but the inferred back azimuths of sources stay stable on the order of days to weeks, we attribute the dominant source locations to moulin tremors and subglacial water routing. At the end of July (when the deployment of all sensors was completed), A1 and A2 tremor sources locate towards the glacier tongue, and A3 tremor sources locate in the region to the west of the array where multiple moulins are located. In addition, we note that seismic tremors are likely generated by efficient channelized subglacial flow (Gimbert et al., 2016) and that the moulins in the vicinity of A3 seemed to evacuate meltwater without the buildup of supraglacial lakes or reservoirs. We therefore suggest that the left branch of the upstream area distributions in Fig. 9, or more general the western and northern part of the glacier, had an efficient and channelized subglacial drainage system. According to Fig. 7c, this configuration stayed stable until the end of August (day 236). At the same time, A0 saw a persistent source to the north whose origin remains elusive. A potential explanation for this source could be another moulin feeding the upstream area branch (Fig. 9) originating in the northeast of the glacier. However, neither the observations from our regular station visits nor the orthophotograph show evidence for a moulin in the area of A0. Starting on day 236, A1 points toward the southwest (Fig. 7), and we attribute this source to the drainage of the smaller supraglacial lake (SL in Fig. 1b). With the onset of the drainage of Lac des Faverges, tremors from the lake basin become dominant (A0 and A1), suggesting that the eastern part of the glacier has an efficient connection to the drainage system, as tremors are expected to originate from channelized flow (Gimbert et al., 2016). Combining this information with the theoretical pattern of subglacial water routing (Fig. 9) suggests that the seismic tremors reveal a gradual “up-glacier” (along the main branch in Fig. 9 from north to south to east) evolution of an efficient channelized drainage system as the melt season progresses. This matches both the field observations (first SL connects to the drainage system, then Lac des Faverges) and the theory of subglacial channel evolution throughout a melt season (Werder et al., 2013).

While subglacial channel evolution is typically described through the competing mechanisms of melting and ice creep (Röthlisberger, 1972), our results show that fracturing can play an important role under specific flow scenarios. We find that icequake activity in the lake basin precedes the drainage onset by several hours (Fig. 6). In combination with a lake reservoir which pressurizes the void spaces and the englacial environment, we suggest that hydrofracturing (e.g., Van Der Veen, 1998; Roberts et al., 2000) drives the drainage initiation. Since no sustained seismicity in the lake region is detected prior to that, this highlights the potential of passive seismic monitoring for early warning of glacier-dammed lake outburst floods. Apart from the lake drainage, other discharge peaks are accompanied by fracturing as well (Figs. 5 and 6). However, we note that elevated strain rate resulting from water-enhanced basal sliding may give rise to icequakes as well (Podolskiy et al., 2016).

5.3 Drainage regime

5.3.1 Theory

Water flow through ice-walled conduits is driven by the hydraulic pressure gradients along the conduits. Along the channel walls, frictional heat enlarges the channels. At the same time, ice creep closes the conduits in the case where the ice overburden pressure exceeds the water pressure in the conduit (Röthlisberger, 1972). These two counteracting processes result in a temporal evolution of conduit radius and water pressure in the conduit. Recently, Gimbert et al. (2016) suggested that pressure fluctuations due to turbulent flow in subglacial conduits can generate seismic tremors whose power scales with discharge according to the drainage regime. Gimbert et al. (2016) derive two end-member scenarios for which the relative seismic power P_rel and relative discharge Q_rel (with respect to some reference state) are related through a power law but with a different scaling exponent.

• i.

  Varying hydraulic pressure gradient and constant hydraulic radius, implying $P_{\mathrm{rel}}\propto Q_{\mathrm{rel}}^{\mathrm{14}/\mathrm{3}}$. As defined in Gimbert et al. (2016), changes in the hydraulic pressure gradient are caused by variations in the water pressure p along a conduit, i.e., $\partial p/\partial x$, where x is the distance along the channel. Such a situation is schematically depicted in Fig. 10, where, for instance, the diurnal melt cycle causes hydraulic head variations in a moulin without changes in the hydraulic radius of the conduit. At some distance from the moulin, at the glacier snout, water constantly flows at atmospheric pressure. As the hydraulic head in the moulin varies, this results in pressure-gradient changes in the subglacial conduit implying $P_{\mathrm{rel}}\propto Q_{\mathrm{rel}}^{\mathrm{14}/\mathrm{3}}$. This drainage regime is expected to dominate in filled subglacial conduits which do not adjust their hydraulic radii fast enough to accommodate discharge changes. We expect that this occurs, for example, for strong daily melt variations in the early melt season (when the capacity of the conduits is still limited) or for rapid water input due to a sudden lake drainage.

• ii.

  Varying hydraulic radius and constant hydraulic pressure gradient, implying $P_{\mathrm{rel}}\propto Q_{\mathrm{rel}}^{\mathrm{5}/\mathrm{4}}$. As the hydraulic radius of a conduit is defined as its cross-sectional area divided by the wetted perimeter, both changes in the water level of a conduit operating under atmospheric pressure and the cross section of a fully filled conduit result in variations in the hydraulic radius (Fig. 10). For instance, subglacial water routing at atmospheric pressure is predicted to be revealed by the power law $P_{\mathrm{rel}}\propto Q_{\mathrm{rel}}^{\mathrm{5}/\mathrm{4}}$ as the wetted perimeter can vary without geometrical changes in the subglacial conduits. The same scaling relationship holds for filled conduits, in case melt enlargement and creep closure of channels dominate over changes in the pressure gradient.

Gimbert et al. (2016) also derived solutions for the relative hydraulic pressure gradient S_rel and the relative hydraulic radius R_rel (again with respect to some reference state) as a function of observed P_rel and Q_rel, given as

$$\begin{array}{}\text{(5)}& {\displaystyle }{S}_{\mathrm{rel}}& {\displaystyle }=P_{\mathrm{rel}}^{(\mathrm{24}/\mathrm{41})}{Q}_{\mathrm{rel}}^{(-\mathrm{30}/\mathrm{41})},\text{(6)}& {\displaystyle }{R}_{\mathrm{rel}}& {\displaystyle }=P_{\mathrm{rel}}^{(-\mathrm{9}/\mathrm{82})}{Q}_{\mathrm{rel}}^{(\mathrm{21}/\mathrm{41})}.\end{array}$$

Figure 10Interpretation of the theory of Gimbert et al. (2016) relating seismic power to discharge. (a) Cross section through a glacier parallel to the flow direction. The hydraulic head in the moulin varies due to, for example, the daily melt cycle. The hydraulic radius of the subglacial conduit is constant. As the water routing at the glacier snout occurs at constant atmospheric conditions, a pressure gradient in the subglacial conduit is present. (b) For such a configuration of varying hydraulic pressure gradient (and constant hydraulic radius) the relative seismic power is predicted to scale with the relative discharge (relative to some reference state) to the power of 14∕3. (c) Cross section perpendicular to the flow direction. The hydraulic radius of a subglacial conduit varies through a change in water level or through changes in the cross-sectional area due to frictional melting or creep closure. The pressure gradient is assumed constant. (d) For such a configuration of varying hydraulic radius (at constant hydraulic pressure gradient), the relative seismic power is predicted to scale with the discharge to the power of 5∕4.

Download

5.3.2 Observations

At A2, we observe tremors due to subglacial water flow beneath Rätzligletscher. Knowing the source locations of subglacial tremors allows us to apply the tremor–discharge relationships to a specific area. If the source locations are not known, the tremor–discharge scalings provide an integrated view over the surroundings of the seismic measurements, whereas the locations presented in Fig. 7 allow us to investigate glacier hydraulics at a specific point, i.e., beneath Rätzligletscher. As Rätzligletscher accommodates the main outlet, we argue that discharge measured in the Simme Valley is representative for water routing at the measured A2 tremor locations. Furthermore, we expect that the number of conduits close to the outlet stays constant. Both assumptions favor the successful application of the tremor–discharge relationships.

Figure 11 shows the scaled seismic power P_rel (square of the tremor amplitude) versus the scaled discharge Q_rel (using the minimum discharge value and its associated seismic power for scaling) on a log-log plot (for details see Appendix C). In this representation, the slope equals the exponent x of $P_{\mathrm{rel}}\propto Q_{\mathrm{rel}}^x$, where the black lines indicate discharge routing accommodated by hydraulic radius adjustment ($x=\mathrm{5}/\mathrm{4}$) and the red lines discharge routing accompanied by variations in pressure gradient ($x=\mathrm{14}/\mathrm{3}$).

Figure 11Tremor–discharge scaling of station PM23. (a) Tremor amplitude (black) and discharge curve (colored) for the pre-drainage period. (b) Scaled seismic power as a function of the scaled discharge (see text for details) on a log-log plot. Color-coding corresponds to the colors in (a). Red and black lines are the drainage-regime predictions of Gimbert et al. (2016) and indicate discharge routing through variations in the hydraulic pressure gradient and variations in the hydraulic radius, respectively (see legend). (c) Distribution of slopes (and thus exponents) calculated from the log-log representation of two adjacent samples each. Black and red bars again show the expected exponents for the two drainage regimes (see legend in b). (d–f) Same as (a)–(c) but for the drainage period.

Download

In the pre-drainage period (Fig. 11a), the power–discharge representation shows a general trend towards radius-adjusting conduits. This is also revealed by the x-exponent distribution (upper right in Fig. 11a) obtained by calculating the slopes between two adjacent samples. This in turn implies that pressure-gradient adjustment occurs rarely and on short timescales only. Such a system is indicative of a well-established, channelized drainage system evacuating water efficiently without significant pressurization. We find such a configuration on Glacier de la Plaine Morte, where the source region of the tremors corresponds to the main trunk of an arborescent drainage network (indicated by the upstream area distributions). However, for the approximately 10 d preceding the drainage but in particular for the last 4 d of this time span with a pronounced diurnal melt cycle, the data suggest pressure-gradient adjustments (yellow dots). This indicates that the capacity of the conduits cannot yet accommodate the water from the melt events without pressurization.

Figure 11b shows the power–discharge scaling for the drainage period. At the drainage onset, the data points scatter along the pressure-gradient adjustment prediction (black dots). Subsequently, after a more chaotic phase associated with clockwise hysteresis, the data reveal hydraulic-radius adjustments during most of the drainage period (purple and orange dots), which is again followed by pressure-gradient adjustments at the end of the drainage (yellowish dots).

To investigate these observations in more detail, we consider the evolution of the hydraulic pressure gradient and the hydraulic radius as calculated from Eqs. (5) and (6), respectively. In Fig. 12, we compare R_rel and S_rel to the measurements of the lake level and the ice surface uplift, which also provide constraints on the drainage hydraulics. In addition, the pictures of the automatic camera provide an estimate of the time when the lake basin was empty. As already inferred from Fig. 11, the diurnal melt cycles prior to the drainage cause pressure-gradient variations while the hydraulic radius changes little. In this phase, the daily peaks of the pressure gradient occur around the time of maximum daily discharge. At the onset of the lake drainage in the evening of day 240 as the lake level starts to drop (gray dashed line in Fig. 12), the inferred pressure gradient increases and reaches its maximum when the rate of ice uplift at A2 is highest. At the same time, the hydraulic radius is described by a transient decrease. Subsequently, the pressure gradient decreases to high pre-drainage values. Concurrently, the hydraulic radius increases as the discharge increases. After the peak discharge, the hydraulic radius decreases again but remains above the pre-drainage level. Subsequent variations in the hydraulic radius and the pressure gradient stay on an elevated level. The sharp peak and drop in S_rel and R_rel on day 243, respectively, correspond to the time when we reinstalled the A2 stations directly on ice, as the snow cover was diminishing. According to the imagery, the emptying of the lake basin was finished in the night from day 246 to 247 (gray dashed line in Fig. 12). A few hours earlier, discharge starts to drop to pre-drainage values. We observe the same for the hydraulic radius. In contrast, the pressure gradient briefly increases before dropping to values lower than prior to the drainage.

Figure 12(a) Change in lake level (orange), discharge measured in the Simme Valley (blue), and vertical ice surface motion at A2 (GPS-2, black dashed). Maximum ice uplift is around 5 cm; see Fig. 2b for scale. The vertical gray dashed lines indicate the start and end times of the drainage, as determined from the lake level change and the automatic camera (as the lake level sensor was not installed at the deepest point of the basin and thus did not provide measurements until the end of the drainage), respectively. Vertical gray bars and roman numbers (I)–(IV) mark snapshots illustrated in the cartoon in Fig. 13. (b) Evolution of the relative hydraulic radius (black) and the relative pressure gradient (magenta) derived from the seismic power and the discharge curve (Eq. 6) for the same time period.

Download

5.3.3 Interpretation

From all our measurements, we deduce the following history of glacier hydraulics associated with the drainage. In the hours prior to the drainage onset, the lake reaches the drainage moulin, but the latter is not yet connected to the subglacial drainage system (situation schematically depicted in Fig. 13a). At this stage, seismic tremors are generated beneath the glacier tongue by the “background” meltwater routing where the daily melt events cause daily variations in the pressure gradient. Through hydrofracturing (Sect. 3.3), the moulin then connects to the subglacial drainage system causing a sudden water input into the drainage system. The lake discharge overwhelms the drainage system, as “an excess of water is pouring into a conduit system of low capacity” (Röthlisberger, 1972), which results in a pressurization of the subglacial environment. From our GPS measurements, it is evident that water pressures exceed the ice overburden pressure, which results in local flotation. The pressure gradient, in turn, can be approximated as the difference in pressure on either side of the tremor-generating region. Considering that water is at atmospheric pressure at the outlet of the glacier tongue, an increase in subglacial pressure due to the lake drainage would also cause an increase in pressure gradient as illustrated in Fig. 13. This is in agreement with our power–discharge-derived pressure-gradient history. Since the conduits cannot adjust their size fast enough, discharge increases only slightly as the lake level starts to drop (Fig. 12). Subsequently, the cross-sectional area (and thus the hydraulic radius) of the subglacial conduits increases due to frictional heat of pressurized flow, causing melting of the ice walls (Röthlisberger, 1972). As the conduits increase in size allowing larger discharge, water is effectively evacuated, resulting in a drop in the pressure gradient and causing the ice uplift to cease (Figs. 12 and 13c). The timescale of conduit enlargement due to melting is expected to be on the order of hours to days (Mathews, 1973).

Figure 13Illustration of the inferred history of subglacial hydraulics associated with the lake drainage. Shown is a schematic section along the major branch of the drainage system shown in Fig. 9. Blue indicates water, and red and black circles indicate seismic wave propagation which indicate discharge routing dominated by hydraulic pressure-gradient adjustments and hydraulic radius adjustments, respectively. The dominant drainage regime after Gimbert et al. (2016) is also given on the right-hand side. The times of the snapshots (I–IV) are indicated in Fig. 12. (a) Situation prior to the lake drainage. The lake reaches the drainage moulin which is not yet connected to the subglacial drainage system, but icequake activity from the direction of the lake basin is increased (indicating hydrofracturing). Tremor generation beneath the glacier tongue is caused by the “background” meltwater routing, and the pressure gradient measured between some arbitrary position along the subglacial conduit and the outlet (constant) is moderate but varies. (b) Initiation of the lake drainage. The drop in lake level causes an increase in the subglacial pressure gradient and local uplift of the ice. The capacity of the conduits is overwhelmed. (c) The subglacial conduits increase their radius by frictional melting to accommodate the lake discharge, which results in a drop in the pressure gradient to pre-drainage values. (d) At the end of the lake drainage, as discharge decreases, the subglacial conduits shrink, causing a short episode of pressure-gradient increase.

Download

As the lake steadily spills water into the moulin, the conduits adjust their size by the competing mechanisms of closure due to ice creep (Nye, 1953; Glen, 1955) and opening due to melting, without significant pressurization and radius changes. Finally, as the discharge drops at the end of the lake drainage, the conduits decrease their size. As the conduits close, another short phase of pressure buildup occurs, indicating the capacity of the conduits is decreased too quickly to maintain constant pressures (Fig. 13d). We expect that conduits tend to close due to ice creep as discharge decreases at the end of the drainage. However, we note that the contraction of conduits takes place on the order of days to weeks, especially for thin ice (less than 100 m) as encountered on Glacier de la Plaine Morte (Mathews, 1973). This suggests that our inferred closure rates of the relative hydraulic radius (Fig. 12) might be overestimated. Another explanation could be the physical collapse of parts of the conduits as discharge decreases (Mathews, 1973). Figures 5 and 6 show that fracturing is indeed pronounced at the end of the lake drainage but we cannot find evidence for strong fracturing from the direction of the glacier tongue, which is expected for mechanical failure during conduit collapse. In addition, the drop in hydraulic radius at the onset of the lake drainage remains enigmatic as we do not have a reason to believe that conduits shrink as an ice-marginal lake starts to drain. We suggest that this drop is an artifact that could be due to neglecting potential changes in channel number and position when inverting for S_rel and R_rel using Eqs. (5) and (6) or by not accounting for sheet-like flow during the ice uplift phase.