Ice-Dammed Lake Drainage Evolution at Russell Glacier, West Greenland

Glaciological and hydraulic factors that control the timing and mechanisms of glacier lake outburst floods (GLOFs) remain poorly understood. This study used measurements of lake level at 15 min intervals and known lake bathymetry to calculate lake outflow during two GLOF events from the northern margin of Russell Glacier, west Greenland. We used measured ice surface elevation, interpolated subglacial topography and likely conduit geometry to inform a melt enlargement model of the outburst evolution. The model was tuned to best-fit the hydrograph rising limb and timing of peak discharge in both events; it achieved Mean Absolute Errors of <5%. About one third of the way through the rising limb, conduit melt enlargement became the dominant drainage mechanism. Lake water temperature, which strongly governed the enlargement rate, preconditioned the high peak discharge and short duration of these floods. We hypothesize that both GLOFs were triggered by ice dam flotation, and localized hydraulic jacking sustained most of their early-stage outflow, explaining the particularly rapid water egress in comparison to that recorded at other ice-marginal lakes. As ice overburden pressure relative to lake water hydraulic head diminished, flow became confined to a subglacial conduit. This study has emphasized the inter-play between ice dam thickness and lake level, drainage timing, lake water temperature and consequently rising stage lake outflow and flood evolution.


INTRODUCTION
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, 1982Clarke, , 2003Fowler, 1999Fowler, , 2009Flowers 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).

STUDY SITE AND MEASUREMENTS
A ∼1 km 2 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 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 −1 (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  Huss et al., 2007;Sugiyama et al., 2007Sugiyama et al., , 2008 Gornergletscher, 2005 1.3 144 60 10 6 Huss et al., 2007;Sugiyama et al., 2007Russell Glacier, 1987  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 timelapse 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 yearround. 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 7 m 3 over 112 days from mid-May to early September (Figure 2A), equating to an average inflow rate of 1.14 m 3 s −1 . This inflow rate is very similar to Russell et al.'s (2011) estimate of 0.6 to 1.3 m 3 s −1 , 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 3 s −1 , but in contrast only 0.86 m 3 s −1 in 2014 and 0.87 m 3 s −1 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 icedammed 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 3 s −1 during flood discharges typically of 1,000 m 3 s −1 ), 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 3 s −1 and 489 m 3 s −1 , 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.

Determination of Lake Outflow Hydrographs
In both 2010 and 2012, the ice-dammed lake water level drained below the elevation of the pressure transducer, as indicated mainly by the sensor temperature record. In order to constrain the later part of the rising limb, it was therefore necessary to reconstruct the peak and falling limb of the outflow hydrograph for both the 2010 and 2012 ice-dammed lake drainage events. Hydrograph reconstruction ( Figure 3A) assumed that the rate of change in lake level drawdown ( Figure 3B) increased linearly until a time t 1 , and then decreased linearly to zero at time t 2 , these assumptions being motivated by the roughly symmetrical shape of the hydrographs measured downstream ( Figure 2B). t 2 was chosen from the period of the flood gauged downstream FIGURE 2 | Lake water level and water temperature record for ablation seasons 2010 and 2012 (A). The gray shaded bars denote the time period depicted in Figures 4A,B. Italicized numbers refer to the maximum thickness of ice that could be liable to flotation given the water level and thus water depth and hydrostatic pressure in comparison to ice overburden pressure. Note varying y axis scales. Panel (B) depicts gauged hydrographs ∼ 23 km downstream at "Sugarloaf" and ∼ 32 km downstream at Kangerlussuaq bridge, both after data by Hasholt et al. (2013) and Mikkelsen et al. (2013).
of 11 and 20 h for the 2010 and 2012 events, respectively Mikkelsen et al., 2013). t 1 was chosen such that the total volume drained ( Figure 3C) within the time interval was consistent with that gauged downstream ( Figure 2B). We acknowledge that flood hydrograph durations typically lengthen with distance down-conduit and down-stream (e.g., Evans and Clague, 1994; so our flood duration is likely to be an over-estimate. Given the assumptions and uncertainty associated with reconstructing the flood recession phase, we will not use the falling limb to constrain our subsequent numerical modeling of the discharge evolution. However, we were motivated to reconstruct the lake water level draw-down in order to (i) estimate a final (post-drainage) lake level and a lake outflow hydrograph ( Figure 4A), and (ii) for comparison with other events at Russell Glacier and elsewhere ( Table 1).
In 2014 our pressure transducer was inoperable after mid-July and in autumn 2015 it could not be recovered due to becoming struck by rockfall, so we cannot constrain the rate of change of outflow in 2014 and 2015. Whilst we could constrain the conduit inlet elevation and the lake volume change in 2014 and 2015 via time-lapse imagery and via dGPS measurements, these together indicated lower volumes and slower rates of outflow than in 2010 and 2012 and so, without the ability to evaluate any numerical FIGURE 3 | Measured lake outflow (black filled squares) with reconstruction (gray filled circles) of the latter part of each event using rate of change of lake level (A), absolute lake level (B) and lake volume (C). Gray shaded areas are indicative of uncertainty range due to pressure transducer and dashed lines indicate uncertainty range due to reconstructions using downstream gauged data. model output, we do not consider these later two events further in this paper.
Overall, the uncertainty in our pressure-level measurements due to sensor resolution and precision is ± 0.05%. Uncertainty in the SETSM DEM elevations for the flight strip used in this study is a maximum of 4 m (Noh and Howat, 2015). Uncertainty in our subglacial bed elevation interpolations and thus ice thickness calculations (Figure 1) is of the order of 20% as a function of the radio echo sounding data collection, interpolation between transects and combination of these data with subaerial topography. Uncertainty in our lake bathymetry interpolations is likely to be ±5%. Given these uncertainties, and the propagation of them (through multiplication of area by elevation change to get volume, and through differencing of DEMs, for example), our field measurements for 11th September 2010 determined a flood volume of 3.1 ± 0.3 × 10 7 m 3 . We assumed that the flood hydrographs were near-symmetrical and our reconstructions yielded a peak discharge of 1,430 ± 150 m 3 s −1 for the 2010 event. For 11th August 2012 we measured a flood volume of 26 ± 0.3 × 10 7 m 3 and reconstructed a peak discharge of 1,050 ± 140 m 3 s −1 .

Model Equations
We model how discharge Q evolves in each flood using a simplified Nye (1976) model of a subglacial conduit (of length l c , cross-sectional area S) draining the lake. The model is "lumped"-it assumes a "short" subglacial conduit and neglects spatial variations, so variables are functions of time t only (cf. Ng, 1998;Bueler, 2014). This approximation holds because l c 1 km and the hydraulic gradient ( ) is sufficiently low so that meltwater from the conduit walls contributes negligibly toward Q, as is confirmed by the small parameter ε = l c /ρ i L ≈ 10 −3 (Ng, 1998;Fowler, 1999;Ng and Björnsson, 2003) in a formal scaling analysis of our system. Although ε does not account for extra melting along the conduit induced by lake thermal energy, our simulation of the flood events later show that this additional melt source is small and does not invalidate the approximation. Specifically we solve the differential equations: for the coupled evolution of conduit size S and lake volume V. Equation (1) expresses the competition between conduit enlargement due to wall melting (which occurs at the rate m) and viscous closure at a rate dependent on effective pressure N = ρ i gHp w , where H is ice-dam thickness next to the lake (55 m), and p w is the conduit water pressure, which is determined by the hydrostatic lake pressure at the conduit inlet; accordingly, p w = ρ w gh, where h is the lake water depth. As the lake drains, h decreases, raising the effective pressure N and enhancing closure. In Equation (2), Q IN denotes the rate of water inflow feeding the lake, which derives from snowmelt and precipitation in the subaerial catchment and any subglacial contribution (or leakage). The standard physical constants in this model are given in Table 2. The closure rate constant K 0 in Equation (1) is based on the creep rate of ice at 0 • C because we assume this part of Russell Glacier to be temperate. Temperature measurements are lacking, but our simulations below show that the closure term is entirely negligible, so cold ice (yielding smaller K 0 ) will make no difference. At any time during flood evolution, we find the lake water depth h from the lake volume V by inverting the lake bathymetry shape function known from field measurements made by Russell et al. (2011): where V is in 10 6 m 3 , z L = h + z LB is the lake surface elevation (m) and z LB = 405 m is the conduit inlet elevation above sea level. Following Nye (1976), the instantaneous discharge is calculated Heat transfer constant F 0 (= 0.205k w (2ρ w /µ w √ π)) 0.8 5,000 kg m −3/5 s −11/5 K −1 Gravity g 9.8 m s −2 Thermal conductivity of water Cuffey and Paterson, 2010) Conduit closure rate constant Water density ρ w 1,000 kg m −3 † Simulation dependent; this was the range explored.
using the Manning Equation for momentum conservation in turbulent water flow: in which the constant F 1 encapsulates the Manning roughness n' and cross-sectional shape of the conduit (see definition by Nye (1976), who symbolized it with a calligraphic-N). Since our field observations show that the conduit lies at the ice-bed interface at its two ends (at the lake and glacier snout), we assume a semicircular cross section for it; then F 1 = (2(π+2)/π) 2/3 ρ w gn' 2 . In Equation (4), is the spatial mean hydraulic gradient driving the flow: g is the glaciostatic gradient = 537 Pa m −1 from topographic data, and the fraction measures the effective pressure gradient along the conduit, which is time-varying because the lake water depth controls N, as mentioned earlier. We expect the flow to transition from closed-conduit (pressurized) to open-channel (atmospheric) some distance back from the exit portal, because thinner ice and the local stress distribution at the snout cause creep closure to become negligible there (Evatt, 2015). N exit is the effective pressure at this transition, and our definition of l c strictly refers to the length of closed-conduit flow. Neither quantity is observationally well-constrained. In our modeling we specify N exit = ρ i g × 35 m, where 35 m (the ice thickness at the transition) is taken from the regional ice thickness near the snout ( Figure 1C). To cater for uncertainty in both this assumption and the conduit's trajectory, we conduct sensitivity analysis in our simulations by varying l c from the estimated curvilinear path distance between the lake and the exit portal (700 m); e.g., l c would be higher for a more sinuous conduit, or less if openchannel flow is more extended. [Additional sensitivity analysis on N exit or the threshold thickness is possible, but not especially insightful, because changes in these can be absorbed into changes in l c ; see Equation (5)].
The use of spatial mean gradients is consistent with our "short conduit" approximation and convenient because the subglacial flood path at Russell Glacier is only approximately located. Nye (1976) and Clarke (1982) also assumed mean hydraulic gradients when applying their models, but neglected the effective pressure gradient (the fraction) from Equation (5). At Russell Glacier, the lake's proximity to the exit portal and high ratio (∼1) of its water depth to the elevation drop along the conduit means that the effective pressure gradient strongly controls when the lake level falls during an outburst. As we shall see, this effect can cause the simulated flood discharge to reach its peak and fall shortly before the lake runs empty.
Finally, we calculate m in Equation (1), which is the mean melt rate along the conduit in our lumped model. We do this by using a formula from Ng et al. [2007;Equations (A3) and (A4)]: where T L denotes the lake water temperature. This formula encapsulates energy conservation and turbulent heat transfer and is based on the pseudo-steady profile of the water temperature along the conduit, which equilibrates on a much faster timescale than flood evolution. Its derivation from Nye's (1976) original model is detailed in the Supplementary Information (SI), where we explain why it captures the heat-transfer mechanisms appropriately, and discuss a concern raised by Clarke (2003) about other formulas for m used in jökulhlaup lumped models. Appearing twice on the right-hand side of Equation (6), the thermal partitioning coefficient α (between 0 and 1) accounts for downstream heat advection, heat transfer between water and conduit wall, and dissipation of potential energy in warming the flow, all of which depend on the changing discharge. As shown in the SI, the instantaneous temperature profile along the conduit is the sum of two temperature distributions: (i) a growing exponential [∝ 1−exp(−βx/l c )] due to potential energy dissipation at rate Q along the conduit, and (ii) a decaying exponential [∝ exp(−βx/l c )] due to water entering the conduit at temperature T L and losing its heat to the walls. The dimensionless parameter β governs the growth/decay rate of these exponentials, and α is the spatial mean value of exp(−βx/l c ) over the conduit. The relationships in Equation (7) show that β ∝ l c / √ Q, and that α is a decreasing function of β such that at the limits β → 0 (a "thermally short" conduit) and β → ∞ (a "thermally long" conduit), α = 1 and α = 0, respectively. Correspondingly, on the right-hand side of Equation (6) are two separate contributions to m from potential head loss (first term) and lake thermal energy (second term). In the first term, a fraction α of Q does not cause melting as it has been used to warm the conduit flow. In our simulations, we assess the relative sizes of these terms during each flood event.

Flood Simulation: Set Up
Equations (1) and (2) are integrated numerically forward in time by explicit finite-difference (Euler) method until h = 0. As the model neglects the flood-initiation process (e.g., Fowler, 1999;Kingslake and Ng, 2013), we specify initial conditions for S and V at a time after flood initiation. A novelty of this study is that these are known coincidently from the "highstand" lake depth h HS in our recorded lake-level histories (Figure 2A). Highstand marks the point after flood initiation when Q is still small but has grown to balance Q IN momentarily (Ng et al., 2007); its timing is not easy to pinpoint due to fluctuations in Q IN and pressuretransducer noise. However, h HS yields accurate initial conditions for V(t = 0) via Equation (3) and S(t = 0) via Equations (4) and (5) (where we set Q = Q IN and N = ρ i gH − ρ w gh HS ). Table 3 lists h HS for the two outburst floods. Without lake-level monitoring data, other ways of constraining the initial conditions would be necessary, e.g., using the observed flood volume together with observed flood peak discharge, as discussed by Ng et al. (2007). We prescribe Q IN as the mean rate of lake-volume rise over the preceding 2 months ( Table 3). We set T L to the mean value of the measured lake-water temperature during each flood, which varied by 0.5 • C in 2010 (between 2.7 and 3.2 • C) and by 0.7 • C in 2012 (between 4.2 and 4.9 • C). Although we did not measure vertical temperature profiles in the lake, the limited (and smooth) variations suggest the lake water to be relatively well-mixed, not strongly stratified. In our model runs, we therefore assume that water temperature measured by the pressure transducer is representative of the water temperature at the inlet portal. Its temporal mean is adopted because the measured variations could be due to circulating lake water with spatially-non-uniform temperature, in which case their size (up to ±0.35 • C) would indicate the potential error of using the temperature at the transducer as the inlet temperature. While we lack observations to support/refute this assumption, sensitivity analysis later ( Figure 6A) shows that such error would alter the simulated peak discharge by several hundred m 3 s −1 if other parameters are held constant.

RESULTS AND ANALYSIS
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 3 s −1 . 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 3 s −1 . 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 3 s −1 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 FIGURE 5 | Best-fit simulated lake level histories (A,B) and flood hydrographs (C,D) for the 2010 and 2012 outbursts when l c = 700 m (black curves) and l c = 500 m (blue). Solid squares show measured data and gray circles reconstructed data. In all panels, the time axis t places the first measured data point used in the hydrograph-fitting exercise at t = 0. In (C,D), the mean absolute error (MAE) of misfit between model and measured data and the optimal Manning roughness n' are indicated.
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.

DISCUSSION
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.

Thermomechanical Contributions
Throughout each simulation in Figure 5, melt enlargement dominates conduit evolution, with viscous closure being negligible in comparison. Indeed, we find that removing or doubling the closure term in Equation (1) merely changes the simulated peak discharges by ≈1 m 3 s −1 . This is partly due to low overburden pressure from a thin ice dam, and partly because other factors precondition an intense melt rate, as will be discussed next. Consequently, the simulated turnaround of discharge at the flood peaks is not caused by closure overtaking melting, but by rapidly decreasing potential energy and hydraulic gradient as the lake depth h falls toward zero. In our model, including the effective pressure gradient (the final fraction) in Equation (5) is crucial for capturing this behavior at Russell Glacier. Accordingly, our simulations show that the conduit continues to enlarge at the times of peak discharge (Figures 7A,B) while the hydraulic gradient causes the melt rate to reach a maximum and then decrease (Figures 7C,D).
In contrast, Nye (1976) envisaged for the 1972 jökulhlaup from Grímsvötn that sudden dominance of an accelerating  viscous closure rate over melt rate (in Equation 1) is what caused flood discharge to attain its peak value in that system. We do not think that this difference in behavior between Russell Glacier and Grímsvötn indicates fundamentally different floodthermomechanical processes in operation in these systems, or that any process has been overlooked or misinterpreted in either system. Instead, it reflects different relative magnitudes of melting and closure. Notably, much greater thickness of the ice dam at Grímsvötn (≈220 m) implies that creep closure exerts a stronger influence on flood evolution of a comparable size to that of melting. On the other hand, the much smaller surface area of the lake at Russell Glacier facilitates its rapid lowering (toward emptying) such that a short flood duration limits the time over which closure reduces the conduit cross-section. Our results here thus highlight the importance of local site factors in GLOF dynamics. Through mathematical analysis of the ratio of the closure timescale to the flood-duration timescale for jökulhlaup systems generally, it should be possible to explain what delineates the two types of behavior (Ng and Björnsson (2003) explored an early theory, neglecting lake thermal energy, for lakes whose depth is a small fraction of the conduit elevation drop).
The finding that melt enlargement dominates outburst flood evolution at Russell Glacier prompts us to examine the contributions to m and S, and this leads us to appreciate how strongly lake water temperature (T L ) determines the floods' size and duration, and the role of a short subglacial conduit in this effect. Equation (6) identifies (i) energy dissipation from potential head loss in the flow (Q ) and (ii) lake thermal energy as the two contributions; these contributions have independent additive effects on conduit wall melting, through their influences on water temperature (see Supplementary Information). Figure 7 shows that the contribution to melt enlargement is dominated by lake water temperature. The reasons for this are two-fold. Firstly, the maximum possible warming of the water from potential head loss, g[h HS + ∆z]/c w = 0.15 • C, is itself an order of magnitude smaller than T L in both floods (several degrees Celsius; Table 2). The actual warming experienced by the water (when it reaches the exit portal) is still less than the maximum possible, because the conduit is "thermally short, " so the temperature addition from "Q " (contribution i) is still increasing spatially toward the maximum. This is reflected by a high value of α ( 0.6) throughout flood evolution (Figure 8), and the attendant low value of (1 −α) multiplying into Q in Equation (6). Secondly, the same factor makes the lake-thermal contribution (contribution ii) efficient because the conduit is short and thus the water temperature has not decreased much from T L by the time it reaches the exit portal. Hence the mean temperature driving heat transfer behind melting along the conduit is high (a large fraction of T L ), and this is reflected by the high α-value in the last term of Equation (6). Despite the resultant high heat transfer, there is a limited decay in water temperature along the conduit because of the short transit time of the flow. In these evaluations, it is important that Equation (6) captures the conduit heat transfer appropriately. This is shown to be the case in the SI, where we also discuss Equation (6) alongside other heat-transfer formulas that had been critiqued by Clarke (2003). As noted before, we have observed that the inlet and exit portals of the conduit are elliptical and with 5 and 15 m vertical and horizontal width dimensions, respectively. The outlet portal survives from one season to the next with "normal" ice meltwater discharge. It seems unlikely that a channel at the ice margin that had grown to 100 m 2 , as suggested by our modeling, would completely collapse between each field season. This is because the ice is thin and whilst we have looked (and surveyed with dGPS) there is no evidence that a large channel was formed, or that one collapsed, such as surface depressions on the ice surface. Nonetheless, whilst some (lateral) advection of the conduit probably occurs between floods, due to ice motion northwards, we suppose that the conduit is surviving from one season to the next. Indeed, in the field we observe "normal" meltwater exiting the outlet portal throughout each spring-summerautumn season. We therefore consider it likely that the conduit seals between floods only in the vicinity of the inlet portal.
One challenge is that at Russell Glacier, and indeed most likely at many other sites of ice-marginal lake drainage, segments of conduits that are at atmospheric pressure for some distance up-glacier from the conduit exit portal are best represented by conditions of open conduit flow. Whilst Nye (1976) did not assume that the conduit was all at ice overburden pressure in his model; i.e., that the effective pressure N = 0, he did assume that dN/dx = 0, where x is the along-conduit dimension. That assumption does not hold for drainages of most ice-marginal lakes where ice dam thickness is just a few tens of meters.

Interpretation of Flood Trigger and Influence on Drainage Evolution
Analysis of ice thickness, bed topography, lake bathymetry and knowledge of lake levels demonstrates that the threshold for icedam flotation being reached was possible and indeed very likely. Specifically, flotation can occur if hydrostatic pressure, ρ w gh, where h is water depth as calculated by maximum water level minus the elevation of the conduit inlet, exceeds ice overburden pressure, ρ i gH where H is ice thickness (e.g., Sturm and Benson, 1985; Tweed, 2000). With the conduit inlet at 404.5 m a.s.l. and given water depth evolution through the ablation season (Figure 2), hydrostatic pressure would have been reached if the glacier ice-dam was ≤45.05 m thick in 2010 and ≤38.46 m thick in 2012. This means that the ice in the vicinity of the conduit inlet had the capability to float in both 2010 and 2012 (Figure 9). Quantifiably, if the lake water surface elevation reached 440 m, as in the year 2012, then lake water depth at the conduit inlet would exceed 27 m and ice with thickness of <30 m, as colored in Figure 9, would be susceptible to flotation.
Our conceptual model for the entirety of both drainage events therefore hypothesizes flotation as the trigger (Figure 10). Notwithstanding our spatial analysis of ice overburden pressure vs. hydraulic head, this hypothesis requires further investigation. We suggest that the two floods were triggered with the lake level lower in 2012 compared to 2010 because the icedam had thinned, for which there is evidence from oblique field photographs, eye-witness reports from the Greenlandic mountain guide Adam Lyberth [pers. comm.], and comparison in this study of ice surface elevation made via ground-based dGPS surveys in 2015 and from a digital elevation model constructed from stereo images acquired in 2011.
Hydraulic jacking due to lake water pressure/head in the vicinity of the conduit, simultaneous with thermal erosion, could have facilitated extremely rapid water egress during the first third of the rising limb of each flood (Figure 10). This contribution of hydraulic jacking might explain the particularly rapid (very high mean outflow rate) water egress at Russell Glacier, in comparison to that recorded for other ice-marginal lakes (Table 1). Thus, rather than a model of a single R-channel, a coupled sheet-conduit model such as that of Flowers et al.
FIGURE 10 | Schematic to illustrate evolution of ice-dammed lake water level and subglacial conduit to explain lake drainage mechanisms. Arrows denote likely subglacial conduit pathway. Contours denote bed elevation.
(2004) could be more appropriate, as has recently been applied by Einarsson et al. (2017) to a rapidly-rising jökulhlaup in Iceland.
As lake water level fell, ice overburden pressure caused remaining outflow to be progressively confined to the conduit. Thereupon melt enlargement became the dominant mechanism permitting discharge increase up to peak outflow (Figures 5C,D,  10). The sharply-inclined falling limb was produced as a function of the low hydraulic gradient, which itself was a function of the high ratio (∼1) of change in lake water depth to the elevation drop along the conduit.

CONCLUSIONS
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).

AUTHOR CONTRIBUTIONS
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.