Field methods and data processing

Fluorescent dyes Uranine and Rhodamine WT were used for the injections. The detection was performed at the gauging station of Grande Dixence SA, 1.25 km below the terminus of the glacier where proglacial discharge was measured. The dye was monitored continuously with both a flow-through and a submersible fluorometer. We present a total of 11 dye injections which were conducted into the lake-drainage moulin in the morning (0830 h local time, UTC+2), early afternoon (1400 h) and evening (2100 h) between 28 June and 18 July 2006. The diurnal variability of tracer experiments conducted using the same moulin can be large (Reference Schuler, Fischer and GudmundssonSchuler and others, 2004); to minimize the influence of such effects, only experiments conducted at the same time of day were compared directly.

We present the tracer transit speed; to characterize the breakthrough curves, an advection–dispersion model with storage (ADSM) was fitted to them (Reference Toride, Leij and van GenuchtenToride and others, 1999). (See Reference Schuler, Fischer and GudmundssonSchuler and others, 2004, for a glaciological application of this model.) The ADSM takes as input the concentration time series and the transit distance. The ADSM returns estimates of the mean tracer transit speed (not presented), the dispersion, D, the fraction of mobile water, β, and the exchange coefficient between mobile and immobile water (not presented). All the breakthrough curves considered in this paper were fitted accurately by the ADSM, and are thus fully described by these parameters. Furthermore, the fraction of returned tracer mass, M, is obtained by integrating the tracer concentration multiplied by the proglacial discharge and divided by the injected tracer mass.

Here we are interested in the time the tracer spends within the glacier, the glacial residence time. Hence it is important to estimate the proglacial residence time, i.e. the time spent in the proglacial stream. Eighteen injections were conducted in the proglacial stream at different discharge levels throughout the field campaign. To the resulting proglacial transit speeds we fitted the expression:

with a least-squares regression. The parameters were estimated to be k = 0.66 m^1−3e s ^{e −1}, Q ₀ = 6.4 m³ s⁻¹, e = 0.34 with a coefficient of determination, R ², of 0.8. The zero proglacial transit speed limit as Q → 6.4 m³ s⁻¹ is unphysical, but the proglacial discharges relevant here are all above Q = 10 m³ s⁻¹. From the fit of Equation (1) to the data, we estimated the error in the calculated proglacial residence time to be ∼30%. The measurement error on the (total) residence time was <1%. This led to a 4% error in the calculated glacial residence time, because the proglacial residence time was around 10 times shorter than the total residence time. The ADSM parameters were not corrected for possible influences of the proglacial stream.

Subglacial water pressure data were obtained from pressure transducers (Gekon 4500) installed in boreholes drilled to the bed (labelled BH1–BH4 and BH6 in Fig. 1). The lake level was monitored with a pressure transducer (Keller DX series). Lake discharge was determined from measurements of the channel cross-section and of the water flow speed (with a current meter Schiltknecht MiniWater20) and the error was estimated to be 10%. The measured discharge agreed well with the independent estimate derived from lake level, bathymetry and water input into the lake calculated with a melt model (using the method of Reference Huss, Bauder, Werder, Funk and HockHuss and others, 2007). A thermistor was immersed into the lake water directly at M4 to measure the temperature. The air temperature and precipitation were measured at an automatic weather station located off-ice at the northern glacier margin (see fig. 1 of Reference Werder, Loye and FunkWerder and others, 2009). To determine the geometry of the moulin shaft of M4, active source seismology was conducted 2 weeks after the presented tracer experiments by lowering explosives into M4. The detonations were detected on the seismic network (Reference Walter, Deichmann and FunkWalter and others, 2008) installed in the vicinity, on the glacier (personal communication from F. Walter, 2008).

Model formulation

Reference ClarkeClarke’s (2003) jökulhlaup model integrates a modified version of the Spring–Hutter equations (Reference Spring and HutterSpring and Hutter, 1982) which describe water flow through R channels. These R channels (Reference RöthlisbergerRöthlisberger, 1972) are en- or subglacial circular or semicircular channels incised into the ice. These equations describe the time evolution of water pressure, channel cross-section, water flow speed and water temperature. The Spring–Hutter equations contain several empirical relations, two of which are of interest here: firstly, the relation between pressure gradient and water flow speed and, secondly, the heat transfer relation.

The Manning–Gauckler–Strickler formula relates water flow speed, v, to hydraulic head gradient, dH/ds:

where n _man is the Manning roughness, s is the along-channel coordinate and R _h is the hydraulic radius defined as

where S is the cross-sectional area of the channel and P _w is the wetted perimeter. The Manning roughness is one of the tuning parameters of the model. The same average roughness is used for englacial and subglacial channels. Reasonable values of n _man lie in the range 0.02 < n _man < 0.08 m^−1/3 s (e.g. Reference ClarkeClarke, 2003).

The amount of melted ice, m, at the conduit walls per unit length is proportional to the temperature difference between ice and water, ΔT:

where L is the latent heat of fusion and U _h is the heat transfer coefficient. The model uses the following empirical relation (Reference McAdamsMcAdams, 1951):

where Re is the Reynolds number, Pr is the Prandtl number, κ is the thermal conductivity of water and P _m is the ice-walled perimeter of the conduit. (Reference Bird, Stewart and LightfootBird and others, 1960, give a list of alternative heat transfer relations.)

As input the model needs the geometry of the connection channel and the thickness of the overlying ice. Channels of different sinuosities, σ, are implemented in the model by stretching the along-channel coordinate, s, by σ. The boundary conditions at the inlet are the water level and the water temperature. At the outlet the water pressure is needed as a boundary condition. Reference ClarkeClarke’s (2003) code has been modified such that the water level at the outlet (normally used to simulate a jökulhlaup which terminates underwater) can be a function of time. The fitting parameters are the initial cross-sectional area of the connection channel, c ₀, its sinuosity, σ, and its roughness, n _man . The model was fitted to the lake discharge and for some runs also to the inferred water flow speeds using, first, a grid search on the tuning parameters (n _man , σ, c ₀) and then a steepest-descent method.

Connection channel location

Only the connection channel, connecting the lake to the main drainage channel of the glacier, is subjected to a jökulhlaup-like evolution (grey arrow in Fig. 1) and hence can be simulated using Reference ClarkeClarke’s (2003) model. We consider a connection channel transit path lying in the region of the grey shaded area in Figure 1. We restrict it to this area since tracer experiments from both M2 and M1 (Fig. 1) show only a modest response to the subglacial outburst floods in 2005 and 2007 (Reference Werder, Loye and FunkWerder and others, 2009), hence connection channels lying north of M2 are not considered. The response of BH4 during the onset period can be explained by a flow direction in the lake dam area parallel to the crevasse orientation, as was also determined by tracer experiments (Reference WerderWerder, 2009). The direction of the crevasses is north–south and BH4 lies upstream of M4. BH3 did not show a response either, so connection channels south of BH3 are not considered, which would also grossly disagree with the direction of steepest descent of the hydraulic potential. Borehole BH2 also does not show a response, which is puzzling as it lies in the middle of the region of possible transit paths. Active source seismology (conducted by lowering explosives into M4) showed that the moulin/channel shaft drops very steeply (crosses on Figs 1 and 5 mark shot locations) and is orientated slightly southwards, suggesting that the connection channel lies south of BH2. A connection channel lying in the grey shaded area (Fig. 1) has a transit distance of .

Fig. 5.

A cross-section of the glacier with the assumed connection channel in grey. The crosses indicate locations of active seismic shots conducted inside M4. The boundary conditions needed at the inlet and outlet are given in parentheses.

The geometry of the channel in the vertical is depicted in Figure 5, which shows a cross-section of the glacier along the connection channel. We assume that the channel drops steeply englacially, as supported by the active source seismology, to reach the bed and then follows it. The length along the connection channel, having a transit distance , is l _con = 790 m. The vertical geometry of channels lying in the grey shaded area is very similar, and thus their geometry can be adjusted to different transit distances by stretching it in the horizontal. This leads to a range of 790 < l _con < 1120 m. To finally arrive at the channel flow path length, we assume the connection channel sinuosity to be constrained by 1 ≤ σ < 2. We absorb the above range of l _con into the sinuosity, fix l _con = 790 m and vary 1 < σ < 2.8, of which a factor of up to 1.4 can be due to different connection channel transit paths.

Flow speed estimate

An estimate must be made of the tracer flow speed in the connection channel in order to be able to compare the water flow speeds from the model with the experimental data. We divide the total residence time, , i.e. the time the tracer took in experiment i to reach the detection station, into the sum of three parts: the proglacial residence time, , the main channel residence time, _, and the connection channel residence time, :

where was measured and can be determined from Equation (1). is unknown but must be greater than 0 to assure a finite flow speed in the main channel. To get an upper bound on , we assume that does not change much due to the additional lake water influx or due to changing meltwater discharge conditions from day to day, as suggested by the constant amplitude and mean of the proglacial discharge (Fig. 4a). Hence an upper bound for is the shortest measured glacial residence time

where experiments conducted in the morning, early afternoon and evening are treated separately, i.e. a different t _shortest is obtained for each time of day. All of them achieved the shortest glacial residence time on 7 July, the day after the onset period. The t _shortest for morning, early afternoon and evening correspond to main channel transit speeds of 0.65, 0.86 and 0.79 m s^{− 1}, respectively.

The inequality on ,

together with Equation (6) leads to upper and lower bounds for :

These bounds on rest only on the assumption that does not exceed t _shortest. To translate this into flow speed bounds which can be compared to Reference ClarkeClarke’s (2003) model, an assumption on the connection channel flow path length needs to be made, as discussed in the previous subsection:

The values of will be referred to as inferred flow speeds or flow speed bounds. Note that the flow speed bounds depend on the sinuosity, σ.

Model set-up

The model is run with the channel geometry discussed above (Fig. 5). The channel cross-section is set to circular for the englacial part and to semicircular for the subglacial part. To reduce the number of free parameters, we use an average Manning roughness and do not distinguish between the roughness of the ice walls and the glacier bed. The boundary conditions at the channel inlet are the measured water temperature (Fig. 4b) and lake level, which was constant during the onset period. At the channel outlet we use the water pressure measured in borehole BH1 and also data from the same borehole but from the previous year, 2005, as the borehole was not so well connected in 2006 (Fig. 4d).

Four different model runs (Mod1–Mod4) are used to investigate the influence of different boundary conditions, physics and fitting procedures. We tune model Mod1 to the lake discharge only and Mod2–Mod4 are tuned, in addition, to the inferred flow speeds. The tuning parameters are the Manning roughness, n _man , the initial channel size, c ₀, and for Mod2–Mod4 also the sinuosity, σ.

It is found that the modelled flow speeds in Mod1 are higher than the inferred flow speeds. If the sinuosity is increased, the inferred flow speed also increases, which is simulated in Mod2 and Mod3. Another possibility is to lower the modelled flow speeds. For pressurized channel flow, the discharge is equal to the flow speed times the channel cross-sectional area. Hence, to decrease the flow speed for a given discharge, the cross-sectional area needs to be larger. This can be achieved by increasing melt of the channel walls (i.e. the heat transfer) or by reducing creep closure. Creep closure can be neglected in these settings as the high water pressure prevents it. This leaves the heat transfer, which is increased by a factor of 2 in Mod4. Table 1 gives a summary of the four model runs performed.

Table 1.

Summary of the settings and results (n _man , σ and c ₀) for model runs Mod1–Mod4

Observations

The tracer experiments show a pronounced evolution of the drainage path of M4 during the onset period (Fig. 2). The increase in glacial transit speed on the first day of the onset period is mainly due to the increased pressure difference driving the water. The further increase on the second day is caused by the enlargement of the channel. Before the lake drained, the channel was inefficient but then developed quickly, as can be seen from the reduction of dispersion, D, and the increased fraction of mobile water, β. However, contradictory to this are the large diurnal water-level fluctuations in BH6 (located only 50 m from M4), which are commonly assumed to be caused by an efficient drainage system. The trace right at the beginning of the onset period shows a greatly reduced fraction of returned tracer mass, M. This might be caused by water being pushed out laterally from the channel into storage areas. On the second day of the onset period the local storage areas may have been filled and the fraction of returned tracer mass increases again to previous values. The data from BH4 (Fig. 3b) also indicate the release of stored water: its water pressure drops 8 hours after the pressure in BH6 drops, which could be due to water flowing out of storage keeping the pressure high. This 8 hour time lag in the signal of BH4 measurements is also present before and after the drainage, indicating that storage-release processes also play a role during normal discharge regimes (cf. Reference Hock, Iken and WanglerHock and others, 1999). The time lag also indicates that BH4 lies upstream of BH6, because otherwise the water pressure in BH6 would lag the pressure in BH4. Hence, water flow does not always follow the direction of the gradient of the calculated hydraulic potential (Fig. 1; also Reference WerderWerder, 2009).

The bounds on the flow speed in the connection channel (Fig. 6) depend on the sinuosity, σ, of the channel and range from a maximum speed of 1.4 m s⁻¹ for σ = 1 to 2.7 m s⁻¹ for σ = 1.86. The inferred speeds on the first day of the onset period are slow (<0.6 m s⁻¹). Combined with the high dispersion, this suggests that the drainage system was not very efficient, as could be the case for a braided flow path or flow through a network of englacial cracks (Reference Fountain, Jacobel, Schlichting and JanssonFountain and others, 2005). The lower bounds on the flow speeds are always very low due to the poor constraint on the upper bound of , but from the comparison to the model, it seems that the upper bounds on the flow speed are more relevant.

Comparison measurements and model

Even though the lake drained by overspilling into a moulin on its shore, the situation was similar to a subglacial outburst flood during the onset period. The main difference is that during a subglacial outburst, the water pressure is higher at the channel inlet which is then located at the lake bottom. However, this difference makes the used model (Reference ClarkeClarke, 2003) no less applicable, as the upper boundary condition can be adjusted accordingly and thus we can test it against our measurements.

Mod1 (Fig. 6a, b) shows that a good fit to the lake discharge can be achieved without matching the flow speeds inferred from the tracer experiments. Hence, discharge measurements alone are not enough to test the Spring–Hutter equations applied to a jökulhlaup. In order to match the inferred flow speeds, either the sinuosity/length of the connection channel (Mod2, Mod3) or the heat transfer (Mod4) need to be increased. However, on the first day of the onset period, the model cannot fit the inferred flow speeds while still matching the lake discharge and the inferred flow speeds of the second day. Thus, there is a qualitative difference in the drainage path between the first and the second day of the onset period. This is supported by the other observations discussed above, in particular the decrease in tracer dispersion. The inability to model the flow speeds on the first day of the onset period may be due to the following reasons: the water does not initially flow through R channels (e.g. braided channels or englacial cracks); the Spring–Hutter equations are not valid for R channels at these low flow speeds and discharges; the geometry changes profoundly between the first and second day (e.g. the sinuosity decreases considerably). Any of these reasons indicates that an accurate model of the very beginning of a jökulhlaup needs to go beyond the Spring–Hutter equations. So far only one jökulhlaup model does this and includes also sheet flow (Reference Flowers, Björnsson, Pálsson and ClarkeFlowers and others, 2004), which probably could capture flow through a more distributed system with suitably chosen parameters.

On the second day of the onset period the model performs well if either the sinuosity/channel length is increased or the heat transfer from the water to the channel wall is increased. A more sinuous or longer channel (Mod2, Mod3) leads to a reduced roughness and higher flow speeds to fit the discharge. However, because the inferred flow speed also depends on the sinuosity of the channel, the model results match up with the measurements. For Mod4, the heat transfer was increased, as it has been noted that the outlet water temperature is overestimated by the Spring–Hutter equations (Reference JóhannessonJóhannesson, 2002; Reference ClarkeClarke, 2003); however, Figure 7c shows that this does not significantly reduce the water temperature. The use of a two-fold enhanced heat transfer coefficient (Equation (4)) leads to a very high Manning roughness (n _man = 0.105 m^−1/3 s) to keep flow speeds low in the larger channel. A high resistance to water flow can be caused by a channel which is low and broad. Reference Hooke, Laumann and KohlerHooke and others (1990) introduced low and broad channels with the shape of a circular segment characterized by its central angle, θ. Such a Hooke channel has a circular cross-section for θ = 360° and a semicircular one for θ = 180°. The hydraulic radius, R _h (Equation (3)), of a Hooke channel becomes smaller with decreasing θ for a fixed cross-sectional area, S. The hydraulic radius in terms of S and θ is given by . A smaller hydraulic radius reduces water flow speed (Equation (2)). Hence, the high Manning roughness of Mod4 can be reduced to an effective value of n _man = 0.049 m^−1/3 s (like Mod2) by the smaller hydraulic radius of a Hooke channel with a central angle of θ = 16.5° (this geometry also leads to higher closure rates, but these are not important in these settings). This channel shape is very low and broad but could occur during the onset of a jökulhlaup if the drainage path evolved from a braided one, as suggested above.

Models Mod1, Mod2 and Mod4 produce virtually identical discharge hydrographs and are distinguished only by their water flow speeds. Hence, the comparison to measured flow speeds is needed for further model validation and leaves Mod2 and Mod4 as options. Unfortunately, it is impossible to discriminate between these two using our experimental data. This is because the flow speeds inferred from the tracer experiments are themselves dependent on the sinuosity. The water temperature at the outlet is almost identical (Fig. 7c) for Mod2 and Mod4, indicating that it is difficult to distinguish high sinuosity from enhanced heat transfer solely by temperature measurements. However, in situations where there are no direct flow speed measurements, low water temperatures at the outlet can indicate high sinuosity or enhanced heat transfer since Mod1 produces higher outlet temperatures than either Mod2 or Mod4. A model of the whole subglacial drainage system would not need the borehole water pressure data as a boundary condition, thus allowing a comparison for validation. However, many more assumptions would be needed to set up such a model and thus it is unlikely to yield more accurate results.

The presented interpretations are based on a model which was deliberately kept simple; nonetheless, there are a few potential shortcomings. An argument was presented to constrain the location of the connection channel to the grey area in Figure 1. However, it is also plausible, from the hydraulic potential, that the connection channel could pass between M1 and M2. Also, the connection channel is likely to have tributaries, influencing its discharge, initial size and pressure. This would lead to increased tracer flow speeds, because in R channels the flow speed is an increasing function of discharge. Hence, the model, not taking into account tributaries, should underestimate rather than overestimate water flow speeds. The water pressure boundary condition at the connection channel outlet is insufficiently known. Figure 6e and f show that the use of the previous year’s subglacial water pressure data with large diurnal fluctuations as a boundary condition leads to variations in the modelled speeds and discharge that are greater than those observed in the measurements. In all likelihood, the diurnal fluctuations found in the pressure data recorded during the onset period were too small, whereas those from the previous year’s data were too great (Fig. 4d). The additional water coming from the lake would keep the subglacial water pressure higher during the night, as was observed, for example, during the 2005 outburst on 11 and 12 June (Reference Werder, Loye and FunkWerder and others, 2009), when the nightly pressure minima became less pronounced at lake discharges comparable to those encountered during the presented measurements. This higher subglacial water pressure at night would suppress the large modelled variations. This sensitivity on the pressure boundary condition at the connection channel outlet shows that the rest of the subglacial drainage system, dictating the subglacial water pressure, can have a great influence on the trigger and initial jökulhlaup evolution. Hence, to simulate the onset of a lake drainage, the prevailing conditions in the subglacial drainage system need to be taken into account. The model does not include storage–release processes, for which there is clear evidence from tracer experiments and borehole water pressure data (see also Reference Huss, Bauder, Werder, Funk and HockHuss and others, 2007). They can influence the pressure and discharge conditions considerably, in particular during the onset. Water diverted into storage should lower the tracer flow speed by the fraction of diverted water. The inferred flow speeds obtained from the tracer measurements on the first day are about two to three times less than the modelled speeds. This could be caused by forcing one-half to two-thirds of the water into storage during this time. On the first day, tracer return mass diminished by half, suggesting that the low flow speed could be caused for the most part by storage processes.