Understanding of sudden and rapid glacier outburst floods, or “jökulhlaups” is important because glacier lakes are increasing in number and size worldwide in mountain regions (Carrivick and Tweed, 2013) and at ice sheet margins, particularly in south-west Greenland (Carrivick and Quincey, 2014). Sudden drainage of ice-marginal lakes can affect local ice dynamics (e.g., Anderson et al., 2005; Walder et al., 2006; Sugiyama et al., 2007; Riesen et al., 2010). Resultant glacier lake outburst floods (GLOFs) can cause intense and extensive downstream geomorphological change (e.g., Carrivick, 2011) and present a hazard to people and infrastructure (e.g., Carey et al., 2012; Carrivick and Tweed, 2016).

Understanding how glacier lakes suddenly drain is challenging, not least due to the different triggers and drainage mechanisms that can act and interact, but also due to a paucity of direct measurements. Potential trigger mechanisms include, but are not limited to, subaerial breaching of ice dams, overspill, ice flotation, syphoning, viscoplastic deformation of the ice dam [otherwise known as “the Glen mechanism” Glen, 1954], changes to the subglacial cavity drainage system, and volcanic activity (e.g., Tweed and Russell, 1999; Björnsson, 2002). However, understanding why and how sudden drainage occurs, remains poorly resolved both theoretically and numerically (Ng et al., 2007; Ng and Liu, 2009). There are a number of theoretical models of ice-dammed lake drainage (e.g., Nye, 1976; Spring and Hutter, 1981; Clarke, 1982, 2003; Fowler, 1999, 2009; Flowers et al., 2004; Kessler and Anderson, 2004; Kingslake and Ng, 2013; Kingslake, 2015) and where models follow the approach of Nye (1976) they often ignore flood initiation, i.e., the flood trigger, by assuming the pre-existence of a conduit whose evolution controls the simulated flood. Ng and Björnsson (2003) concluded their analysis of ice-dammed lake drainage by stressing: (i) the importance of identifying the flood trigger, (ii) a need for more monitoring of lake levels during floods and (iii) reliable measurements of ice dam thickness and lake geometry.

The primary aim of this paper is to analyse the thermomechanical evolution of a rapidly draining ice-dammed lake, taking advantage of a suite of direct field measurements. This study was motivated generally by (i) the unusually high mean lake outflow rates and (ii) the consequent rapid discharge of recent glacier outburst floods or “jökulhlaups” at Russell Glacier, west Greenland (Table 1).

A ~1 km² ice-dammed lake on the northern flank of Russell Glacier, west Greenland (Figure 1) is known to have drained repeatedly from the late 1940s until 1987 (Sugden et al., 1985; Russell and de Jong, 1988; Russell, 1989). Following 20 years of relatively stable lake level, a jökulhlaup on August 31st 2007 marked renewed ice-dammed lake drainage (Mernild et al., 2008; Mernild and Hasholt, 2009) and a new jökulhlaup cycle (Russell et al., 2011). To date, this new cycle has resulted in floods almost every year (Table 1). This lake has a history of being studied for its lake drainage mechanism(s) (e.g., Russell and de Jong, 1988; Scholz et al., 1988; Russell et al., 2011).

Ice surface elevations were surveyed in the field in 2010 and 2015 using a Leica GPS500. We used a temporary base station, positioned relative to the Kellyville International Geodetic System network continuous receiver, via post-processing of a 10-h static occupation and recording at 1 min intervals, and with a vertical precision of ± 0.01 m. From this temporary base station our rover points, or points of interest on the ice surface (Figure 1), which were obtained in real time static mode (using the geometric mean of 120 static readings), have a vertical precision of ± 0.15 m.

Further distributed ice surface data were gained from a 2 m digital elevation model (DEM) that was produced by photogrammetric processing of stereo-pairs of DigitalGlobe imagery, specifically via the Surface Extraction with Triangulated Network-based Search-space Minimization (SETSM) algorithms (Noh and Howat, 2015). Note that whilst denoted as from year 2011, this is a composite DEM, the seamless coverage being constructed from multiple image pairs from multiple flight lines from multiple dates. In order to consider any relationship between change in lake level and change in ice dam thickness, elevation differences between our 2015 dGPS field surveys and the 2011 DEM have been converted to annual rates of change in surface elevation and these show surface lowering in the vicinity of the subglacial conduit that is more than double that across the rest of the ice margin (Figure 1).

Glacier bed topography, indicated by contour lines in Figure 1, was obtained by combining 1 km resolution ice thickness data (Bamber et al., 2013), a regional InSAR-derived ice surface (Palmer et al., 2011) and local IceBridge flightline data (Leuschen et al., 2017) and was gridded at 0.1 km cell size. This bed topography shows that the Russell Glacier margin abuts the ice-dammed lake on an adverse gradient bed slope (Figure 1). Recognition of the fact that the glacier bed is inclined downwards in a cross-flow direction is important for understanding ice margin dynamics, the configuration and evolution of the subglacial drainage network and subglacial water pressure, all of which we refer to later.

A seamless digital elevation model was constructed by combining the 0.1 km interval subglacial bed contours with the 5 m proglacial terrain elevation contours and interpolating the remaining ice-dammed lake floor using the ANUDEM algorithm, which is within the ArcGIS “Topo to Raster” tool. This seamless digital elevation model enabled us to make spatially distributed calculations of hydrostatic pressure and ice overburden pressure for scenarios of varying water depths and varying ice surfaces or ice thickness.

The ice-dammed lake drains (almost every year; Table 1) via a subglacial conduit with stable inlet and exit positions (Figure 1). From these positions, we derived an approximate straight planform length of 650 m and a conduit gradient of ~0.015 m m⁻¹ (Figure 1). The hydraulic gradient of this conduit is low, due to the relatively shallow ice surface and bed topography slopes (Figure 1). Field observations and photographs of the conduit inlet and exit portals show that the conduit is elliptical in cross-section and has a diameter (at both inlet and exit portals) of ~ 5 m vertically and ~ 15 m horizontally; direct measurements have never been possible because both portals are unsafe to access. Lake basin bathymetry (Figure 1) is well-known from differential Global Positioning System (dGPS) surveys conducted when the lake basin was almost completely drained (Russell et al., 2011).

In anticipation of ice-dammed lake drainage, a HOBO U20 pressure transducer, with range from 0 to 76.5 m and an accuracy of ± 0.05%, was installed in mid-May 2010 (as indicated in Figure 1) when the lake level was at 432.5 m a.s.l., which is 27.5 m above the conduit inlet portal. The transducer was weighted to ensure it remained on the lake bed. The HOBO recorded average lake water pressure and lake water temperature at 15 min intervals. Lake drainage events were recorded by the pressure transducer on 11th September 2010, 11th August 2012, by time-lapse cameras on 3rd August 2014, and by the same cameras accompanied by dGPS measurements on the 28th July 2015. The 2010 and 2012 lake water pressure records were corrected for atmospheric air pressure recorded 2 km away from the ice-dammed lake with an automatic weather station (AWS) using Campbell Scientific sensors operating continuously year-round. The AWS data showed that seasonal weather up to and including the drainage events was not unusual, thereby assisting in ruling out any external trigger to the drainage events. A simple calculation of the melt season up to the day of lake drainage, made with the AWS data, gives 116, 75, 46, and 60 positive degree days for 2010, 2012, 2014, and 2015, respectively.

In 2010, the ice-dammed lake filled by 1.1 ± 0.1 × 10⁷ m³ over 112 days from mid-May to early September (Figure 2A), equating to an average inflow rate of 1.14 m³ s⁻¹. This inflow rate is very similar to Russell et al.'s (2011) estimate of 0.6 to 1.3 m³ s⁻¹, which was derived from supraglacial stream gauging only. Calculations performed on the 2012 lake level data suggest a lake refill rate of 2.4 m³ s⁻¹, but in contrast only 0.86 m³ s⁻¹ in 2014 and 0.87 m³ s⁻¹ in 2015. The 2012 refill rate is presumably high because of the extreme ice surface melt experienced over the Greenland Ice Sheet (GrIS) (Nghiem et al., 2012) during this year. Comparison between these lake refill rates and stream inflow rates demonstrates that there are no substantial subglacial inflows to the ice-dammed lake, that subglacial leakage from the ice-dammed lake during lake filling is negligible, and thus that the seasonal and event water balance(s) can be considered simple and closed.

These refill rates and our interpretations of water inflows to the ice-dammed lake are important because the few previous studies that have obtained direct measurements of ice-dammed lake drainage have found it necessary to deal with complications in the event water balance. These complications have primarily been: (i) significant upstream inputs to the lake (e.g., at Kennicott Glacier: Anderson et al. (2003); at Grímsvötn: Flowers et al. (2004); at Merzbacher lake where Ng et al. (2007) reconstruct inflows of up to 100 m³ s⁻¹ during flood discharges typically of 1,000 m³ s⁻¹), and (ii) release of significant volumes of water from storage within the glacier over and above that draining from the lake; e.g., Huss et al.'s (2007) calculations imply that up to 40% of a jökulhlaup from Gornersee came from non-lake water within Gornergletscher. In comparison, the Russell Glacier ice-dammed lake drainage events have a mean rate of water release that is up to two orders of magnitude greater than other (non-volcanically triggered) measured sudden lake level falls (Table 1). Indeed the 2010 and 2012 events at Russell Glacier had mean outflows of 426 m³ s⁻¹ and 489 m³ s⁻¹, respectively (Table 1), which is equivalent to ~10 and ~11 Olympic swimming pools draining per minute, respectively. These mean outflows exceed the normal summer ice melt discharge observed (but not gauged) through the exit portal by two orders of magnitude. They also have a relatively simple water balance; (i) subglacial water inputs to the ice-dammed lake throughout the summer were negligible as evidenced by lake refill rates that agree well with gauged supraglacial meltwater runoff, as described above, and (ii) inflow to the lake was volumetrically negligible during the lake drainage events since both events lasted just hours rather than days.

Numerical runs are made with different combinations of Manning roughness n' and (closed) conduit length l_c to simulate each year's hydrograph for comparison with the observed hydrograph. Since the equations are autonomous and the lake highstand too loosely located in time for us to put the “t = 0” precisely on the observational time line, in each comparison we always first align the hydrographs in time by sliding the simulated one against the observed one until the Mean Absolute Error (MAE) in Q between them is minimized for their overlap period. The MAE is then reported as a percentage of the mean observed discharge for the period. This procedure uses only measured discharge data on the rising limb and ignores reconstructed data (as these are more uncertain), so it is not aimed at fitting the flood peak.

For both 2010 and 2012, the map of MAE across the parameter space of n' vs. l_c (Figures 4A,B) implies a weak influence of conduit length on the success of fit. MAEs <10% are found along the valley axes (Figure 4C). Within a plausible range of l_c (400 to 1000 m), neither map shows a minimum point in MAE that identifies the best combination of n' and l_c. The weak influence of l_c found here means that it cannot be constrained based on the fit on the rising limb alone.

For each choice of conduit length, the minimum MAE identifies an optimal n' allowing best fit of the rising limb of the hydrograph (white curves in Figures 4A,B); we call the corresponding numerical run an “optimal simulation”—for that l_c. For the l_c range studied here, optimal values of n', ≈ 0.03 to 0.045 m^−1/3 s, fall within the range encountered in jökulhlaup studies [~0.01 to ~0.1 m^−1/3 s Ng, 1998; Fowler, 1999; Clarke, 2003; Werder and Funk, 2009; Kingslake, 2013]. Figure 5 shows the modeled lake level and outflow hydrograph histories in the optimal simulations for l_c = 700 m, and Table 3 lists the corresponding n'. The simulated lake levels match the observed histories well for most of the flood duration (Figures 5A,B). In both years, the model captures the shape of observed discharge on the rising limb successfully (MAE = 3.0% in 2010 and = 5.6% in 2012). The simulated peaks occur slightly later than expected from our reconstructions and overshoot the reconstructed peak discharge by 100 to 200 m³ s⁻¹. For both years, these matches to the reconstructed peak discharge are improved if we assume a shorter section of closed-conduit flow, e.g., l_c = 500 m (blue curves, Figures 5C,D). Then the simulated peaks are lowered to near the reconstructed peaks, while the MAEs for the rising limb remain acceptable (4.3% in 2010 and 4.0% in 2012, Figure 4C). These findings suggest that there might have been a considerable (≈ 200 m long) section of open-channel flow behind the exit portal in both the 2010 and 2012 events; however, we emphasize that this inference is very tentative because our peak discharge values are reconstructed, with large uncertainties of ±140 to 150 m³ s⁻¹. On the other hand, this exercise of varying l_c shows that a matching procedure that best-fits the peak discharge as well as the rising limb can constrain l_c as well as n'—at least for the parameter region studied here. An accurate and complete measured hydrograph would be needed for this purpose. Note that the choices of l_c in these experiments are illustrative: we are not implying that l_c is the same in different floods. However, our modeling assumes l_c to be constant in each simulation/flood, and this may not be an accurate description.

The modeled floods terminate with Q falling to zero abruptly, in a manner noted by Clarke (2003), when the lake level falls below the conduit inlet elevation (Figure 5). While this describes the situation at the lake, remaining water in the conduit takes time to evacuate, and an expansion wave causes Q at the exit portal to decrease more gradually (Fowler and Ng, 1996). For each year, the flood volume released by the lake is nearly invariant across our simulations because the model uses the same initial lake level (h_HS), and the cumulative inflow into the lake (which is drained in the flood) is very small. The flood volume at the exit portal is slightly higher due to the added conduit meltwater, but only by ≈ 0.2%, so the short-conduit approximation used in the model is self-consistent.

Our sensitivity experiments reveal an acute influence of lake water temperature on the flood hydrograph (Figure 6A), which is stronger than that of other factors within their conceivable range of variations (Figures 6B–E). The same qualitative finding was reported by Werder et al. (2010) in their modeling study of a supraglacial lake that drained through englacial/subglacial pathways. Warming the lake by a couple of degrees Celsius raises the peak discharge and shortens the flood substantially, whereas T_L near 0°C stretches the hydrograph to a much lower peak over a week or more (Figure 6A). This contrast is similar to that between volcanically-triggered (rapidly-rising) jökulhlaups and the more common (slower-rising) jökulhlaups from Grímsvötn (Björnsson, 2002; Roberts, 2005). We conclude that T_L of several degrees Celsius, coupled with the short ice dam, are the most major factors controlling the peak magnitude, the timing of peak magnitude, and the short flood durations (<1 day) in 2010 and 2012. These same factors are likely responsible for similar outburst characteristics in other years at Russell Glacier (Table 1).

The sensitivity results in Figures 6B–E provide additional insights for the simulation of the Russell-Glacier floods in the contexts of hydrograph matching and prediction. The simulated hydrograph is largely insensitive to the lake inflow rate Q_IN (Figure 6E), so our using a Q_IN-estimate from the months before each flood to approximate Q_IN during the flood is reasonable. This approximation may serve well in other purely predictive runs, although dedicated sensitivity tests should be performed in each case. The peak region of the simulated hydrograph is sensitive to the highstand lake depth h_HS (Figure 6B), with the peak discharge varying by roughly ±50 m³ s⁻¹ per meter error in h_HS (a similar result is found in 2012 and 2010, for both l_c = 500 and 700 m). Thus, lake-level monitoring to determine h_HS with sub-meter accuracy is necessary for the model to match/predict the flood peak correctly to within this uncertainty range—our pressure transducer easily meets this requirement. The results for l_c and n' (Figures 6C,D) inform the question of constraining these parameters through matching the observed flood hydrograph, which we encountered earlier. Whereas the peak of the simulated hydrograph and the curved shape of its rising limb are both sensitive to n', only the peak is sensitive to l_c. This finding explains why the MAE of fit on the rising limb depends weakly on l_c (Figure 4), why the procedure of minimizing this particular measure of error constrains n' but not l_c, and why matching the peak discharge also (if this is known reliably) should constrain l_c as well as n'. These interpretations from Figure 6 are based around a specific set of model parameters; model sensitivity far from this region of the parameter space may be very different.

Our numerical simulations yield new insights into the thermomechanical characteristics of the 2010 and 2012 floods. When these insights are considered with the information from our other datasets, we are able to consider the entirety of each sudden drainage event, not just the rising limb, and also the most likely trigger of drainage.

Understanding the controls on the timing, magnitude and thermomechanical dynamics of GLOFs is complex. This study reports a rare set of comprehensive field measurements from an ice-dammed lake in west Greenland and has used them to evaluate the likely lake drainage trigger and to quantitatively assess the mechanisms of drainage evolution.

Following hypothesized flotation, we propose that hydraulic jacking sustained rapid subglacial water egress during the first third of the rising limb of both floods. Our numerical modeling of the 2010 and 2012 outbursts at Russell Glacier shows that after flood initiation, the rate of melt enlargement of the flood-conduit walls was controlled strongly by lake water temperature. Lake temperature and falling lake water level both sensitively affected flood hydrographic evolution, including peak discharge. Viscous closure of the subglacial conduit was negligible due to the short duration of each flood and because the overlying ice is relatively thin. Classical studies of jökulhlaup systems where the overburden ice thickness/pressure is considerably higher (e.g., Grímsvötn) show that viscous closure overtaking melting can cause the turnaround of discharge at the flood peaks, but this is not case for the two Russell Glacier floods. Instead, their turnarounds are due to diminishing hydraulic gradient, preconditioned by the high ratio of lake-level change to elevation drop along the conduit.

Our analysis has emphasized inter-play between ice dam thickness and lake level, drainage timing, lake water temperature and consequently rising stage lake outflow and flood evolution. We have also shown how different thermodynamical factors behind the total outflow, and their relative importance, evolved during each of the drainage events. Our quantification of the control of lake temperature on lake outflow reinforces the notion that warming of lake water by rising air temperatures can affect glacier lake outburst flood (GLOF) timing and flood magnitude (Ng et al., 2007).

JC conceived of the study, led all data collection and analysis and the manuscript write up. FT assisted with the original concept and the writing and completely edited the manuscript. FN made all of the mathematical formulation, all of the modeling/simulations, and wrote the discussion of modeling results and the SI. RH assisted with preliminary modeling. All the other authors assisted by collecting data in the field or by providing complimentary datasets.