Introduction
Outburst floods from ice-dammed lakes constitute a major natural hazard and can lead to extreme discharge events (Reference Clague and MathewsClague and Mathews, 1973; Reference HaeberliHaeberli, 1983; Reference Walder and CostaWalder and Costa, 1996). In densely populated mountainous areas, these lakes often represent a severe threat (Reference Haeberli, Alean, Müller and FunkHaeberli and others, 1989). Lakes impounded behind ice barriers can drain in a variety of ways. The term ‘jökulhlaup’ refers to catastrophic drainage by rapid enlargement of subsurface tunnels. This phenomenon has been investigated thoroughly both through observations and theoretical studies (Reference NyeNye, 1976; Reference Spring and HutterSpring and Hutter, 1981; Reference ClarkeClarke, 1982; Reference BjörnssonBjörnsson, 1992, Reference Björnsson1998, Reference Björnsson2003; Reference Fountain and WalderFountain and Walder, 1998; Reference AndersonAnderson and others, 2003; Reference RobertsRoberts, 2005; Reference Huss, Bauder, Werder, Funk and HockHuss and others, 2007). Ice-dammed lakes can also drain by mechanical collapse of the ice dam (Reference HaeberliHaeberli, 1983) or can overflow the ice barrier and drain over an ice-floored spillway, through a drainage channel or through an aerial side channel between the ice and the rocky edge (Reference Walder and CostaWalder and Costa, 1996; Reference Raymond and NolanRaymond and Nolan, 2000; Reference Mayer and SchulerMayer and Schuler, 2005).
In this study, we investigate the breaching of an ice dam accompanied by thermal erosion of the drainage channel. More specifically, we focus on the case where the discharge water melts and lowers the floor of the drainage channel faster than the lake level drops. In the Alps, the breaching of ice dams has led to a number of catastrophic outburst floods (Reference HaeberliHaeberli, 1983). The outburst flood from Glacier de Gietroz, Switzerland, in 1818 is one of the best-known disasters related to glacier hazards. In order to drain the large volume of water impounded by the advance of Glacier de Gietroz, the authorities decided to dig a channel through the ice. Unfortunately, unexpected rapid erosion within the channel led to a sudden emptying of 20 × 106 m3 of water and caused 40 fatalities (Reference HaeberliHaeberli, 1983). In 1985, the Bogatyr ice-dam lake in Kazakhstan was drained artificially (Reference Mochalov and StepanovMochalov and Stepanov, 1986; Reference Nurkadilov, Khegay and PopovNurkadilov and others, 1986), leading to an unexpected outburst flood of >7 × 106 m3 of water in 2 days with a peak discharge of 100 m3 s−1. More recently, other ice-dammed lakes, on Grubengletscher, Switzerland, Ghiacciaio del Belvedere, Italy, and Glacier de Rochemelon, France, have caused concern to authorities and scientists (Reference Haeberli, Kääb, Vonder Mühll and TeysseireHaeberli and others, 2001, Reference Haeberli2002). In all these cases, the possibility of digging an artificial channel in the ice was considered. Before carrying out such an operation, it is important to predict the rate of erosion of the ice and the water discharge rate that can be expected through the channel. The mechanics of the erosion process depend on numerous parameters and are still unclear. The main reason for this lack of knowledge is the lack of field measurements during the breaching of ice dams. Most of the time, the only available data relating to lake drainage through ice channels consisted of rough estimates of lake volume and peak discharge (Reference Walder and CostaWalder and Costa, 1996). Some measurements of ice erosion are available from experiments carried out in a cold room (Reference Costard, Dupeyrat, Gautier and Carey-GailhardisCostard and others, 2003; Reference Isenko, Naruse and MavlyudovIsenko and others, 2005).
In 2005, the artificial drainage of Lac de Rochemelon through an ice barrier gave us the opportunity to perform extensive field measurements during drainage. To our knowledge, this is the first time that discharge and water temperature have been monitored continuously during the drainage of a lake through an ice dam, with repeated measurements of ice erosion.
The goal of this paper is to (1) present the dataset obtained from these field measurements, (2) simulate the discharge and erosion process using a numerical model and (3) study the sensitivity of the discharge to water temperature and other parameters. This study should help provide a better assessment of the risks related to lake drainage through a channel in an ice dam.
Studied Area and Historical Background
Glacier de Rochemelon (45°13′ N, 7°4′ E) is located in the Vanoise area close to the Franco-Italian border. It spans an elevation range of 2950–3300 m, with an area of 1.6 km2 and a northwest exposure. It is located 10 km upstream of the nearest town and infrastructures. A supraglacial lake started to form about 50 years ago and grew steadily until 2004. Aerial photographs taken by the French National Geographic Institute indicate that the lake formed between 1953 and 1962. The rate of change of the surface area of this lake is easily assessed from aerial surveys carried out in 1953, 1962, 1970, 1980, 1986, 1996 and 2000. The lake was located in the upper part of the glacier, at 3218 m a.s.l., and was bounded to the southwest by a rocky ridge (Fig. 1). It was dammed by the glacier on the northeast side. Until September 2004, the lake water was naturally drained through an outlet in the rocky ridge, thereby controlling the lake level. However, the ice dam was thinning, and geodetic measurements performed between 2001 and 2004 showed that the ice dam had decreased by an average of 1.3 m a−1. Field measurements were carried out on 31 August 2004 by the Laboratoire de Glaciologie et Géophysique de l’Environnement (LGGE). They showed that the volume of the lake was about 650 000 m3 and that the freeboard below the top of the ice dam was 1.5 m, threatening several towns in the valley below. The LGGE warned the authorities in charge of public safety (Préfecture de Savoie) on 16 September 2004 (M. Vallon, unpublished report; C. Vincent and E. Le Meur, unpublished report). The freeboard decreased to 0.8 m on 17 September and to 0.2 m on 5 October. Immediate action was therefore taken to lower the lake level by 6 m before the beginning of the winter. This work was supervised by the Restauration des Terrains en Montagne (RTM) service. The company in charge of the work used a siphoning technique to drain the water in October 2004. It was subsequently decided to drain the lake entirely the following summer. The job was carried out step by step between June and August 2005. The first step was to siphon the water out of the lake using pipes crossing over the ice dam. No pumping was required. Because of the low atmospheric pressure at this altitude, it was only possible to siphon 6 m. Once the lake level had been lowered by this amount, a 1 m wide channel was dug into the ice using explosives. When the bottom of the channel had reached the lake water level, pipes were installed in the channel and the water was siphoned again. The lake level was thus reduced without any risk of overflowing the ice dam. On 25 August 2005, when the volume had been reduced to 250 000 m3, it was decided to use explosives again to dig the channel to a level below that of the lake surface, creating an ice channel through which the water could flow naturally (Fig. 2). During the discharge, extensive field measurements were performed to study the channel erosion process.
Fig. 1. Lac de Rochemelon (a) on 22 August 2004, and (b) on 15 August 2006 a year after the artificial drainage. (Photographs by M. Caplain and P. Macabies.)
Fig. 2. Channel dug through the ice during drainage. The top 6 m were dug artificially using explosives, then the water overflowed through the ice breach and dug the channel naturally. (Photograph by C. Vincent.)
Field Measurements During Drainage
Photogrammetric measurements based on aerial photos taken in 2000 provide an accurate map of this region (Fig. 3). In addition, geodetic and bathymetric measurements were carried out in 2004 and 2005 to produce a very accurate digital elevation model (DEM) of the ice dam and the bottom of the lake (Fig. 4). The accuracy of the geodetic and bathymetric measurements is 5 and 50 cm respectively. At the beginning of the drainage operation, on 25 August 2005, the volume of the lake was 250 000 ± 28 000 m3. Radar and seismic measurements were performed to determine the dam ice thickness (M. Vallon and others, unpublished report; C. Vincent and others, unpublished report). The measured ice thickness and surface topography along the channel are shown in Figure 5. The longitudinal profile of the channel bottom was measured during drainage on 25 August, 26 August, 29 August and 23 September 2005 (Fig. 5). Six ablation stakes were set up on the glacier in the vicinity of the lake to measure the contribution of ice melt to the filling of the lake during drainage. Water levels (Fig. 6) were continuously monitored using two autonomous pressure gauges. These gauges were accurate to within 0.1%, i.e. 2 cm for 20 m depth, and were set up at the bottom of the lake. The water-pressure change was corrected from simultaneous atmospheric pressure measurements. In addition, the lake water level was checked 25 times during the drainage operation using independent readings made on a leveling staff set on the bottom of the lake. After 4 days of measurements, a discrepancy of 3 cm was observed between these two independent methods. As a result, the outflow was calculated from the lake-level changes, taking into account the lake surface-area changes and the recharge rate due to ice melt (Fig. 6). The lake hypsography is known accurately from geodetic measurements performed after the lake was drained, and the surface area can easily be determined for each elevation. The recharge rate of the lake was assessed from the melt rates obtained from the ablation stakes set up in the catchment area. Note that the recharge rate never exceeded 10% of total discharge. Based on the uncertainties in the hypsography, lake-level measurements and recharge rates, the discharge accuracy is estimated to be ±0.1 m3 s−1.
Fig. 3. Map of Lac de Rochemelon. The DEM is based on photogrammetric measurements from aerial photos taken in 2000. The surface area of the lake comes from measurements made in 2004. As can be seen, the lake was bounded to the southwest by a rocky ridge and to the northeast by an ice dam. Until September 2004, the lake drained naturally through the outlet in the rocky ridge. The inset shows the area depicted in Figure 4.
Fig. 4. Detailed map of the lake, ice dam and channel dug through the ice. A, B and C refer to the cross-sections at which the water height and elevation change of the channel bottom were measured. The ablation stakes are represented by black triangles.
Fig. 5. Bottom profiles of the drainage channel at different dates (2005). The lake level corresponds to that of 24 August 2005. The origin of the horizontal distance is the breach through which the water escaped into the channel.
Fig. 6. (a) Lake temperature. (b) Lake level (continuous curve) and channel bottom elevation (dots). (c) Discharge deduced from the water-level measurements. The discharges were corrected for the ice-melt inflow estimated over the drainage basin. The time origin corresponds to the beginning of the drainage operation.
The water temperature was also monitored continuously (Fig. 6) by two autonomous thermistors, with an accuracy of 0.1°C, installed at the bottom of the lake. Three other thermistors were set up at 2.2, 7.2 and 12.2 m depth. Seven temperature measurements were made by thermistor before and after the beginning of the drainage. These showed that the differences between thermistors did not exceed the measurement uncertainty (±0.1°C). Note that the water temperature lapse rate was <0.01°C m−1. Note also that the lake was completely ice-covered until the beginning of July, with a water temperature of 0.2°C. It was partially ice-covered at the beginning of the drainage operation on 25 August and ice-free on 26 August.
Ice erosion in the channel was measured by the elevation change of the channel bottom at cross-sections A–C (Fig. 4). These measurements, as well as the water height measurements in the channel, were performed 18 times at these cross-sections during the drainage using a theodolite and a leveling staff (Fig. 6). For these measurements, the channel bottom was accessed by abseiling. Elevation changes of the channel bottom were also measured at cross-sections D, E, F and G at the beginning of drainage on 25 August, during drainage on 26 August and at the end of drainage on 29 August. Temperature measurements in ice boreholes revealed a temperature of −2°C at 12 m depth.
Analysis of the Energy Fluxes Dissipated into the Channel
As water escapes through a breach in an ice dam, the energy dissipated by the flow leads to melting of the ice. However, the mechanisms involved in this ice-dam erosion process are unclear. To study the efficiency of erosion, we first analyze the energy loss through the channel. The melting of the channel floor is related to the energy lost by the flowing water through energy conservation. Consequently, some of the potential energy and thermal energy of the water is transferred to the ice for melting. The question is how much. For numerical modeling purposes, Reference Raymond and NolanRaymond and Nolan (2000) and Reference Mayer and SchulerMayer and Schuler (2005) assume that the potential energy and thermal energy are transferred to the ice floor of the channel uniformly over its length according to the average slope and rate of water cooling. Thus, they assumed that temperature decreases linearly to freezing point at the exit of the channel. However, this assumption is not based on observations and is subject to question. In our study of Lac de Rochemelon drainage, the lake temperature was measured continuously. Also, the elevation of the bottom of the channel was measured three times at cross-sections D–E, E–F and F–G to determine the total ice ablation in the channel (Fig. 5). From these measurements, it is possible to compare the available energy with the energy loss resulting in ice melt (Fig. 7). According to the conservation of energy over a longitudinal section, the output energy equals the input energy corrected for the change of internal energy. This can be expressed in the following form:
where E t is the thermal energy of water, E p is the gravitational potential energy, E k is the kinetic energy of water, E m is the energy used for channel ablation (melting or mechanical abrasion), E a is the heat exchange with the atmosphere and E i is the dissipation by thermal conduction of heat through the ice walls beyond the melted slab of ice. Subscripts 1 and 2 respectively refer to the entrance and exit of the longitudinal section over which the energy balance is computed. Energy exchanges with the atmosphere, E a, which can correspond to a loss or gain of energy according to the temperature of the atmosphere, are neglected (Reference Walder and CostaWalder and Costa, 1996). Heat conduction through the walls of the channel beyond the ice-melting boundary, E i, is also assumed negligible (Reference Walder and CostaWalder and Costa, 1996). By expressing the difference of gravitational potential energy ΔE p, Equation (1) can be rewritten:
This equation has been applied to six 15 m longitudinal sections along the channel (Figs 5 and 7) for different time periods, making it possible to compare the energy loss resulting in ice ablation and the available energy flux at the end points of the section. The entry of the first section is the breach through which the water escaped from the lake into the channel. Consequently, the first term of Equation (2) is reduced to E t1 + ΔE p because of null kinetic energy at the entrance. The left-hand side of the equation includes the two sources of energy that will be used to melt the ice at the bottom of the channel. The energy used for channel ablation, E m, is calculated from our erosion measurements in the channel (see below). The terms E t2 and E k2 depend on water temperature and water flow at the exit of the section. They remain unknown, but their sum E t2 + E k2 is obtained from the other terms. They represent the remaining energy after the loss due to melting and constitute the input energy for the following section.
Fig. 7. Relationship between the energy used to melt the ice at the bottom of the channel and the available energy from both water temperature and potential energy. The energy quantities have been calculated over six 15 m longitudinal sections (A–B, B–C, C–D, D-E, E-F, F–G). The corresponding measurement dates (26/08 is 26 August, etc.) are reported.
The thermal energy E t (J) is obtained from the mass of water discharged, water temperature and specific heat of water, which can be written as:
where Δt (s) is the duration of the measurements, Q (m3 s−1) is the discharge over time Δt, ρ w (kg m−3) is the water density, c w (J kg−1 °C−1) is the specific heat of water and θ w (°C) is the water temperature.
The change in gravitational potential energy ΔE p (J) is deduced from the channel slope. It is given by:
where g is acceleration due to gravity (N kg−1) and Δh is the elevation difference between the two ends of the longitudinal section along the channel.
The energy loss E m (J) is obtained from the average melt rate measured along the longitudinal section and the latent heat of fusion. It is given by:
where m (m s−1) is the average melt rate measured at the bottom of the channel in the longitudinal section, S (m2) is the bottom surface area of the longitudinal section, ρ i (kg m−3) is the density of ice and L f (J kg−1) is the latent heat of fusion.
This implies that the erosion is due to melting and that mechanical abrasion is neglected. The importance of mechanical abrasion as a process for enlarging ice walls is unknown (Reference ClarkeClarke, 2003) and here we implicitly assume that the heat loss by mechanical abrasion is included in the heat flux due to melting.
These values have been calculated for six 15 m longitudinal sections and for different time periods, depending on available data (Fig. 5). At each section, the energy used upstream for melting was subtracted from the available energy. As a result, a clear relationship can be determined between energy loss and available energy (Fig. 7). For each 15 m longitudinal section, the energy loss is 7.3% of the available energy. This rate seems to be similar whatever the section or time period. Consequently, on average, the energy loss per meter along the channel is only 0.5% of the available energy. Since the channel is 110 m long, the energy used for melting is 55% of the available energy. A similar calculation over the entire period (4 days) and over the entire length of the channel, using the mean discharge, mean temperature and total melt volume, confirms this result. In addition, note that the contribution of potential energy is low, representing at most 7% of the available energy. From these results, we conclude that about half of the available energy is used for breach erosion.
Numerical Modeling
We now assess the ability of a numerical model to reconstruct the ice erosion and water discharge observed through the channel during drainage. However, the main objective is to study the sensitivity of ice erosion and discharge to water temperature, channel width and lake surface area and volume. For hazard analysis, it is essential to predict the ice erosion rate and the maximum expected discharge. The following numerical model is based on ideas proposed by several studies (Reference ClarkeClarke, 1982; Reference Walder and CostaWalder and Costa, 1996; Reference Raymond and NolanRaymond and Nolan, 2000; Reference Costard, Dupeyrat, Gautier and Carey-GailhardisCostard and others, 2003) and provides an estimate of drainage through an ice channel versus time once water overflows the ice dam.
First, the discharge through the breach can be estimated from well-known relationships in hydraulics. For a breach with a rectangular cross-section, assuming that critical flow is reached (Froude number >1), the discharge Q (m3 s−1) is given by Reference LencastreLencastre (2005):
where l (m) is the width of the breach, h (m) is the lake level relative to the bottom of the breach and μ is the discharge coefficient. The initial geometrical parameters of the breach are provided as initial conditions.
The water depth D (m) in the breach is given by Reference LencastreLencastre (2005):
Thus, the mean velocity u (m s−1) of water is obtained from the discharge (Equation (6)) and wetted cross-sectional area Dl:
The melt rate m (m s−1)at the bottom of the channel can be written as:
where q c (W m−2) is the heat flux at the water/ice interface, θ f (°C) is the temperature of fusion and θ i (°C) is the temperature of the ice.
The heat flux transferred to the ice is expressed as a function of the heat transfer coefficient h t (W m−2 °C−1) and the water temperature θ w (°C) (Reference Costard, Dupeyrat, Gautier and Carey-GailhardisCostard and others, 2003):
The heat transfer coefficient h t characterizes the heat transfer from the turbulent flow to the ice wall. It is given by:
where Nu is the Nusselt number, k w is the thermal conductivity of water (W m−1 °C−1) and L c1 (m) is a characteristic length. Reference Walder and CostaWalder and Costa (1996) suggest that L c1 can be considered equal to the part of the wetted perimeter on which ice melting occurs. Assuming vertical erosion, we take L c1 = l.
The Nusselt number is calculated from an empirical relationship using the Prandtl and Reynolds numbers, Pr and Re respectively. It reads:
where A, α and β are empirical coefficients. The Prandtl number Pr is the ratio of the dynamic viscosity to the thermal diffusivity of the water. For a water temperature of 0°C, Pr is fixed at 13.2. The Reynolds number is given by:
where η w is the kinematic viscosity (m2 s−1). L c2 (m) is equal to 4R h, where R h (m) is the hydraulic radius (Reference Walder and CostaWalder and Costa, 1996).
The heat flux resulting in ice melting (Equation (9)) is calculated from Equation (10) using the heat transfer coefficient. This coefficient is calculated from Equations (11–13). The ice melt rate m makes it possible to calculate the new elevation of the bottom of the breach. Finally, the change in lake level, dz (m s−1), is obtained from the flow balance and the surface area of the lake:
where A l (m2) is the surface area of the lake and Q i (m3 s−1) is the recharge rate (inflow) of the lake. The hypsography of the lake is well known from bathymetric measurements, and the surface area can be determined easily for each elevation. As explained above, the recharge rate of the lake was assessed from the ablation stakes set up in the drainage basin. The changes in lake and breach elevations result in a new value of h which is then introduced in Equation (6) for the next step. The time-step in our numerical experiments is 60 s. The values of the constant parameters are shown in Table 1. The coefficient α has been set to a value of 0.333 based on previous studies (Reference Lunardini, Zisson and YenLunardini and others, 1986; Reference Costard, Dupeyrat, Gautier and Carey-GailhardisCostard and others, 2003). If we used a value of 0.4 for α, i.e. the value chosen by Reference McAdamsMcAdams (1951), the results would be very similar. The coefficients μ, A and β were deduced from our measurements. The discharge coefficient value μ results from our discharge and channel geometry measurements. For this purpose, the channel elevation and width measurements, performed 17 times during the drainage once the breach had a regular shape, have been used. The discharge coefficient was then inferred from the relationship between the discharge Q and the value of l(2g)1/2 h 3/2 (Equation (6)). From these measurements, the discharge coefficient has been calculated to be 0.18. Note that the measurements carried out over the first 12 hours of drainage have not been taken into account because the shape of the breach was not yet regular. For a concrete reservoir with vertical walls, the discharge coefficient ranges from 0.4 to 0.55 (Reference LencastreLencastre, 2005). Our data show that the discharge coefficient for ice is much lower than these theoretical values.
Table 1. Physical constants and model input parameters
The coefficients A and β used to calculate the Nusselt number were determined empirically in order to obtain the best agreement between calculated discharge and melting and the corresponding measurements. Coefficients A and β used by Reference McAdamsMcAdams (1951) and Reference Lunardini, Zisson and YenLunardini and others (1986) lead to an underestimation of melting and discharge. The best agreement was found with A = 0.332 and β = 0.74, and the Nusselt number then reads:
The value of α = 0.4 used by Reference McAdamsMcAdams (1951) would lead to A = 0.24 and β = 0.75 and the results would be very similar. Since the glacier is cold, part of the energy is used to warm the ice to melting temperature by diffusion of heat through the melting slab (see Equation (9)). The extra flux that diffuses beyond the melted ice (through the channel walls) has been assumed negligible (Reference Walder and CostaWalder and Costa, 1996). The difference between ice temperature and fusion temperature leads to a melting decrease of 3%. Heat loss to the atmosphere is assumed negligible. The melted ice contributes <0.5% to the net flow and is thus also negligible. The drainage started on 25 August at 1900 h. Simulated and measured data are compared in Figure 8. At the beginning of drainage, the shape of the channel was not regular since it had been dug with explosives. For this reason, the numerical modeling experiment starts 5 hours after the beginning of the drainage. The irregular shape of the breach over the first 12 hours could explain the difference between simulated and measured data. Twelve hours after the beginning of drainage, the channel was about 1.3 m wide and very regular. The discharge variations observed in the channel seem to be related to water temperature changes (Fig. 6) due to the diurnal cycle. The water warming observed on 26 and 28 August between 1000 and 1800 h led to an increase in discharge, although this effect is not clear on 27 August.
Fig. 8. (a) Lake-level and channel-bottom elevation change and (b) channel discharge from measurements and numerical modeling. The time origin corresponds to the beginning of the drainage operation.
As stated above, coefficients used for the Nusselt number have been adjusted to obtain the best agreement between simulations and measurements. Consequently, the above results do not prove anything about the ability of this numerical model tool to predict ice erosion and discharge through an ice channel. However, it is probably a good tool to study the sensitivity of ice erosion and water discharge to the water temperature, channel geometry and lake surface area and volume.
Discussion
The drainage discharge depends essentially on a complex combination of water temperature, channel geometry, lake hypsography and recharge rate, with strong feedbacks. A water temperature rise leads to an increase in breach melting and discharge. In addition, the rate at which the lake level drops depends on the discharge and the surface area of the lake. High discharge can occur when the discharge melts and lowers the breach floor faster than the lake level drops. This can occur when the lake surface area is large. Conversely, if the surface area is small, the lake elevation drops faster than the breach elevation. In this case, the discharge will decrease along with the ice melt.
Our case study is relatively simple, as the channel geometry is very regular with time and the recharge rate is very low. Consequently, Lac de Rochemelon drainage depends essentially on two input variables: water temperature and lake surface area. Some relevant features are noted when comparing measurement data and modeling results. Figure 8a shows a strong change in the rates at which the lake level and channel bottom are lowered, about 53 hours after the beginning of drainage. During the first period, from the beginning of the drainage until 0200 h on 28 August, the lake level decreased at a nearly constant rate of about 9.5 cm h−1. Given the decreasing lake area, we would have expected an acceleration of this rate. In fact, the decrease of water discharge due to a decrease in melt rate more than offsets this effect. After this first period, the lowering of the lake level suddenly decreased to 6.2 cm h−1 until 0500 h on 30 September, at the end of drainage. On the other hand, the discharge fluctuations show a different pattern. Except for the large fluctuations in the first 24 hours of drainage, the discharge decreased steadily from 1900 h on 26 August to 1100 h on 28 August. Note the subsequent strong change in the rate of change of the discharge, 9 hours after the change in lake lowering rate. During this intermediate period, the lake level continued to drop faster than the channel bottom despite the decrease in lake lowering rate. Consequently, the discharge continued to decrease. This corresponds to low water temperatures, <0.8°C. When the discharge reached about 0.35 m3 s−1,at 1100 h on 28 August, the lowering rates of the lake and channel were similar and the discharge was almost constant over the next 15 hours. The model successfully simulates the timing of these changes in lowering rate and discharge. During the last period, the discharge increased following the warmer temperature of 29 August. The modeling results are generally consistent with the measured data except for this last increase in discharge which is poorly simulated. This discrepancy could be related to the uncertainty in the lake hypsography since discharge is increasingly sensitive to hypsography as the lake empties.
More generally, the drainage discharge is very sensitive to the water temperature. From numerical modeling experiments, it can be shown that a temperature increase of 0.1°C could lead to a difference of 0.15 m3 s−1 in the peak discharge and to a difference of 2.5 cm h−1 in the lake lowering rate. The drainage discharge is also very sensitive to the channel width. If the channel width decreased by 0.1 m, the peak discharge would increase by 0.1 m3 s−1. However, the channel width remained almost constant during drainage, and the ice walls in the channel were more or less vertical. The reason for this constant width remains unclear. Melting would be expected to occur on all the ice in contact with the water and should therefore cause an overhang to develop. In reality, however, melting appears to be concentrated at the bottom of the channel. In addition, similar channel widths have been observed in other natural drainage events of supraglacial lakes in the Alps (e.g. Ghiacciaio della Croce Rossa, Italy, in 2004 (personal communication from L. Mercalli, 2007) and Glacier du Baounet, France, in 2005).
In September 2004, the crucial question was posed as to what would happen if the lake overflowed the ice barrier. To answer this question, we have performed numerical experiments with the lake full and with different water temperatures. For this experiment, the initial water depth at the breach was fixed at 10 cm. The results are reported in Figure 9. For all experiments, the drainage starts very slowly and then rises exponentially for the highest temperatures. With a water temperature of 1°C, the peak discharge will reach 1.5 m3 s−1 and drainage will take about 10 days.
Fig. 9. Numerical modeling results for a full lake at the beginning of drainage for different water temperatures (+4°C, +2°C and +1 °C). (a) Lake level (continuous curves) and channel bottom elevation (dashed curves). (b) Discharges calculated with different channel widths (thin curves: 0.5 m; thick curves: 1.0 m; dashed curves: 1.5 m).
Given that the measured water temperature was +1°C on 30 September 2004, it is likely that this discharge would have reached this value if the lake had overflowed the ice barrier at this time. However, as observed in August 2005, the temperature can rise by almost 0.5°C on a warm day. With water temperature of 1.5°C, the peak discharge will reach 2.8 m3 s−1 and drainage will take 5 days. Although not consistent with previous observations, a channel width reduced to 0.5 m (and a temperature of 1.5°C) will lead to a discharge of 6 m3 s−1 and drainage will take 2 days. As a result, it can be concluded that natural drainage of this lake at the beginning of October 2004 probably led to a discharge of about 1.5 m3 s−1 and is not likely to have exceeded 6 m3 s−1.
As seen in Figure 9, further experiments with warmer water show a strong sensitivity to water temperature and to the width of the channel. The question arises as to whether these temperatures can be reached in this glacier-dammed lake. During summer, when the lake surface is ice-free, the water is warmed under the influence of solar radiation and air temperature. In addition, the warming is enhanced by thermal convection. As the density of fresh water is maximum when temperature reaches +4°C, the warmer surface water sinks to the lake bottom, flows toward the ice front and cools. The density therefore decreases and the water rises again to the lake surface. Thus it is believed that the water temperature cannot exceed +4°C when the lake is dammed by an ice barrier or when the water is in contact with ice. Temperature data from glacier-dammed lakes in other studies show temperatures ranging from 0 to +4°C, supporting this assumption. The water temperature was 0.7°C for the marginal lake of Black Rapids Glacier, Alaska, USA, in June 1993 (Reference Raymond and NolanRaymond and Nolan, 2000), 2°C for periglacial lake No. 5 at Grubengletscher, Switzerland, during the summer of 1994 (Reference Haeberli, Kääb, Vonder Mühll and TeysseireHaeberli and others, 2001), 2–4°C for Lac d’Arsine, France, in 1971 (Reference TournierTournier, 1971), 0.5–1.1°C for periglacial Gornersee lake, Switzerland, in 2007 (personal communication from M. Werder, 2008), 2°C for supraglacial Grindelwald lake, Switzerland, in 2008 (personal communication from M. Werder, 2008), 2–4°C for Lake George, Alaska (Reference Walder and CostaWalder and Costa, 1996), 0.2–1.2°C for Hidden Creek Lake, Alaska (Reference AndersonAnderson and others, 2003), and 3–4°C for the Bogatyr ice-dam lake, Kazakhstan (Reference Nurkadilov, Khegay and PopovNurkadilov and others, 1986). Using the maximum temperature of +4°C, the drainage of Lac de Rochemelon would have led to a peak discharge of 11 m3 s−1 (for a width of 1.3 m).
Finally, further experiments have been performed to highlight the sensitivity of drainage to the lake surface area (Reference Raymond and NolanRaymond and Nolan, 2000). These show that with the same volume, a temperature of 1°C and a channel width of 1.3 m, the peak discharge of Lac de Rochemelon would have reached 4.2 m3 s−1 if the surface area of the lake had been twice as large as it was in 2004 (compared to 1.5 m3 s−1).
Conclusions
The potential danger of water bodies that form at the surface of glaciers can generally be determined at an early stage and monitored. The main problem is to estimate the risks related to uncontrolled overflow.
Outburst floods can result from either supraglacial flooding, subglacially routed flooding or a mechanical collapse of the ice dam. In the case of Rochemelon, the outburst from the lake could have been triggered by overtopping and consequent erosion of the ice dam, subglacial water drainage and rapid enlargement of the ice channel close to bedrock or mechanical collapse of the ice dam.
In all cases reported in the Alps (Grubengletscher; Belvedere; Rochemelon), the possibility of digging an aerial channel in the ice has been considered in order to drain water and avoid any risk of outburst flooding. This method appears to be the easiest and least expensive way to solve the problem. Consequently, the question of predicting the rate of erosion of the ice and the water discharge through an aerial channel is crucial.
The artificial drainage of Lac de Rochemelon provided a very good opportunity to investigate the breaching of an ice dam accompanied by thermal erosion of the drainage channel. Using extensive field measurements, the available gravitational and thermal energy and the energy used for melting have been compared. The results show that only half the available energy has been used for breach erosion.
In addition, ice-dam breaching was analyzed using a numerical model calibrated with the numerous measurements made during the drainage operation. Overall, the model results are consistent with the measurements made during drainage. As expected, the numerical model experiments reveal a strong sensitivity to the water temperature and the width of the channel. The model was used to simulate discharge and ice erosion with the lake full, as was the case in September 2004 when it was decided to drain the lake artificially to avoid any risk of an outburst flood. According to the model, natural drainage would have led to a peak discharge of 1.5 ± 1.5 m3 s−1, probably representing no danger to life and property downstream.
Note that the surface area of the lake is a key parameter when assessing the risk of outburst flooding. For a given water volume, the danger increases with the surface area of the lake.
Acknowledgements
The field campaigns were funded partly by the ‘Service de Restauration des Terrains en Montagne (RTM) de Haute Savoie’ in France. We thank the Sens Alpin company and RTM for pressure-gauge measurements in the lake. We thank B. Laily, P. Macabies and C. Maréchal for advice and for taking part in collecting the field measurements. We thank M. Vallon whose comments and advice improved the quality of this study. We are very grateful to M.T. Gudmundsson, G.K.C. Clarke and an anonymous reviewer whose comments greatly improved the manuscript.
