Here we present the results from a range of techniques which were used to observe and model the growth and drainage of a supraglacial lake and water‐filled crevasses at Helheim Glacier, southeast Greenland. The observations highlight an unusual pattern of filling and draining of these areas of surface water, not typically observed on the Greenland ice sheet. We first assume that the observed drainage is driven by hydrofracture, and attempt to model this behavior using a linear elastic fracture mechanics model. However, we find that this model cannot explain the observations. We therefore propose an alternate hypothesis which is supported by analysis of the basal hydrological conditions in the Helheim catchment.

Despite the observed links between lake drainages, basal water pressure, and flow speed, observations of the subglacial system following lake drainages are still limited, particularly at tidewater glaciers. A substantial volume of work has investigated the subglacial hydrology of land‐terminating glaciers through observations [e.g., Hubbard and Nienow , 1997 ; Bartholomew et al. , 2010 ; Chandler et al. , 2013 ; Cowton et al. , 2013 ] and modeling [e.g., Banwell et al. , 2013 ; Dow et al. , 2015 ]. However, the subglacial systems of tidewater glaciers are typically much more difficult to access using techniques such as boreholes, which are not suitable for the highly crevassed surface of large tidewater glaciers. Further observations are vital in order to better understand the dynamics of tidewater glaciers, as well as ice‐ocean interactions and the impacts on fjord circulation [ Straneo et al. , 2013 ].

Supraglacial lakes, which frequently cause hydrofracture [e.g., Das et al. , 2008 ; Danielson and Sharp , 2013 ; Tedesco et al. , 2013 ], are most commonly found in the south‐western region of the Greenland ice sheet (GrIS), with only 2% of lakes by number occurring in the South East [ Selmes et al. , 2011 ]. Studies and modeling work on the behavior of supraglacial lakes have therefore mainly focussed on the southwest [e.g, Box and Ski , 2007 ; Banwell et al. , 2013 ; Clason et al. , 2015 ]. Remote sensing investigations in this region have found that supraglacial lakes typically drain at progressively higher altitudes as the melt season progresses [ Sundal et al. , 2009 ; Doyle et al. , 2013 ; Morriss et al. , 2013 ], a process which has been reproduced in modeling work [ Arnold et al. , 2014 ; Clason et al. , 2015 ]. The up‐glacier progression of drainage has also been observed in the Canadian Arctic, where water‐filled crevasses close to the terminus of Belcher Glacier were observed to drain by hydrofracture much earlier in the melt season than lakes at higher elevations [ Danielson and Sharp , 2013 ].

Research has shown that hydrofracture can easily force a crevasse to penetrate through the full thickness of an ice sheet [ van der Veen , 2007 ], rapidly transporting large volumes of surface meltwater to the bed [ Das et al. , 2008 ] and leading to increases in flow speed on diurnal [e.g., Shepherd et al. , 2009 ] to seasonal [e.g., Bartholomew et al. , 2010 ] timescales. These increases in flow speed may be driven by high water pressures in the subglacial system, which reduce basal friction and lead to rapid sliding [ Iken , 1981 ; Iken and Bindschadler , 1986 ; Meier et al. , 1994 ]. This is a widely observed phenomenon on alpine [ Hubbard and Nienow , 1997 ] and land‐terminating glaciers [e.g., Bartholomew et al. , 2010 ]. At marine‐terminating margins, this mechanism has been observed at Helheim Glacier, where ice velocity lags surface meltwater production by 1 day [ Andersen et al. , 2011 ], as well as on the west coast of Greenland [ Sole et al. , 2011 ], and in Alaska [ Kamb et al. , 1994 ; Oneel et al. , 2001 ].

Hydraulic potential in Helheim catchment was calculated from the Greenland Mapping Project (GIMP) surface DEM [] and the IceBridge BedMachine Greenland Version 2 bed DEM [], both at a spatial resolution of 150 m. Based upon the work of], the hydraulic potentialwas calculated usingwhereandare the surface and bed elevations, andis a fraction of overburden pressure, typically set to 1 based upon the assumption that the entire catchment is at overburden pressure. We relaxed this assumption in order to test how variations in this fraction affected the hydraulic potential within the catchment []. We use values ofbetween 0.5 and 1.2 in order to test the likely range of basal effective pressures.

We use a constant background stress field in all years. This appears to be reasonable based upon the sensitivity tests, which show that the model is insensitive to the background stresses. However, in order to test this we also ran the model on stress fields from images pairs collected in November 2010, March 2011, and July 2012. The results showed variations of no more than 1 day between years, consistent with the sensitivity tests. This suggests that the background stress field does not have a significant impact on the results from year to year.

Two potential errors are introduced by the flow routing assumptions used here. First, we do not account for flow delay and retention in the catchment, potentially leading to an assumption of early hydrofracture. However, as can be seen from Table 2 and Figure 1 , the catchments are relatively small, average flow paths are up to 5.8 km for area L and much smaller for other areas, and the catchments largely become snow free early in the melt season, leading to minimal retention. The maximum flow path length is 12 km for area W 1. Based upon a conservative assumption that water travels at an average of 0.1 m s −1 within the catchment, the longest flow path would only take around 36 h to drain the surface water areas. Therefore, we expect this to have a minimal impact on the time of hydrofracture. Second, by spreading the water volume evenly over the maximum water area we underestimate the maximum water depth in the early stages of filling, potentially leading to an assumption of later hydrofracture than in reality. This is therefore a conservative assumption for the purposes of this model.

The combined sensitivity test allows all parameters to vary at random within the ranges defined above. This provides an estimate of the overall uncertainty in the model and is used to define the error bars in result plots. For the combined test, the results show an IQR of 14 days. The 5th and 95th percentiles show that 90% of the model results lie within ± 15 days of the modeled date of hydrofracture. Given the wide ranges within which the parameters are allowed to vary, it is unlikely that the error of the model is outside this range.

Field measurements and observations of crevasse widths forced by hydrofracture are limited. A maximum crevasse width of 0.4 m was measured by Doyle et al. [ 2013 ], following hydrofracture of a lake on the west coast of Greenland. Further research by Krawczynski et al. [ 2009 ] suggests that these widths may be up to 1–2 m for very deep cracks (>1500 m) or for high longitudinal stress gradients. In order to capture this variation, we tested a range of values between 0.1 and 2 m, taking a base value of 1 m.

Values of K I C between 100 and 400 kPa m 1/2 were used by van der Veen [ 1998 ] based upon previous laboratory testing of glacier ice [see van der Veen , 1998 , and references therein]. We extend this range to cover the values used by Mottram and Benn [ 2009 ] who tested values as low as 10 kPa m 1/2 , which may be more appropriate for weaker ice close to the terminus of Helheim. The range of values for K I C used here therefore cover the range 10–400 kPa m 1/2 , with 150 kPa taken as a base value consistent with previous work.

Varying the tensile stress has a very similar impact to varying the crevasse spacing, in that the size of the first term in equation 1 is increased or decreased relative to the other terms in the equation. It therefore has a similar impact on the modeled results; crevasses penetrate deeper during spin‐up, but the IQR of 1 day shows a very minor impact on the day of hydrofracture.

The parameter space for the tensile stress was defined by allowing for errors in the selection of the value of A ; a temperature dependent flow parameter used in the creep relation (equation 4 ). Previous work has assumed an ice temperature of −5°C [e.g., Clason et al. , 2012 ], which gives A a value of 9.3 × 10 −25 s −1 Pa −3 [ Cuffey and Paterson , 2010 ]. We took the limits of our parameter space as 0°C and −30°C which give values of 2.4 × 10 −24 s −1 Pa −3 and 3.7 × 10 −26 s −1 Pa −3 for A , respectively. The effect of the choice of this parameter is to alter the mean stress in the catchment from a base value of 303 kPa at −5°C to lower and upper bounds of 220 kPa and 889 kPa for 0°C and −30°C, respectively. We therefore took these as the limits of our parameter space in the Monte Carlo sensitivity testing.

The results of the sensitivity testing show that the model is insensitive to the crevasse spacing. The IQR of the parameters tested is 1 day, which is equal to the model time step. Closer inspection of the model results shows that increased crevasse spacing leads to an increase in the initial crevasse depth forced by the background stress field, but once runoff is added hydrofracture occurs at much the same time.

For the sensitivity tests performed here, a ratio of crevasse spacing to depth of 1:9 was taken as a lower bound estimate of the spacing at Helheim, giving S a value of 0.1. As crevasse spacing increases , therefore S = 1 was taken as an upper bound on the range of S . In the absence of more detailed information, we also used S = 1 as the base value for other model runs, but as shown by the sensitivity results this had a minimal impact on the results.

The effects of crevasse spacing were tested by] through modification of the tensile stress term in equation 1 such that the first term becomeswhereis the ratio of the crevasse spacing to the sum of the crevasse depth and the spacing, such thatas the spacing decreases.) is then an empirical function ofwhich varies between 0.5 and 1.12. For the full definition of these factors, see].

The sensitivity tests show that the model is highly sensitive to the runoff, with the interquartile range (IQR) of the difference between modeled and observed results equal to 9 days, the second largest IQR of the parameters tested. The extended parameter range shows that for α f <0.2 hydrofracture frequently fails to occur, therefore indicating a lower limit on the possible catchment size.

The largest source of uncertainty in the runoff arises from the catchment delineation. Smaller catchments lead to later hydrofracture, up to a threshold where there is insufficient water for hydrofracture to occur. However, while small catchments may explain later hydrofracture, they cannot explain why water is only found in discrete areas on the surface of the glacier. Based upon these limitations we set the lower bound of the runoff as a 50% decrease in catchment area, and the upper bound as a doubling in catchment size. Further to this, we tested an extended parameter range for runoff which represents a 90% reduction in the catchment size. Such small catchments are difficult to reconcile with observations; however, it is informative to demonstrate the model behavior under these conditions. A multiplication factor α f , consistent across catchments, was used to alter the runoff in the model using values between 0.5 and 2, and a base value of 1, with the extended case using a lower limit of 0.1.

Two types of test were performed. In the first, a single parameter was varied while all others were kept at their base values, thus allowing the sensitivity of the model to each individual parameter to be tested. In the second type of test, all parameters were varied simultaneously, allowing the overall uncertainty in the modeled time of drainage to be tested. In each test the model was run for 5000 different randomly selected parameter combinations. The range and base values of each parameter are discussed below and summarized in Table 3 . The results are discussed and presented as the difference in days from the base case, and the interquartile range (IQR) is used as a metric to compare the relative sensitivity of the different parameters.

Sensitivity tests were carried out on an idealized model setup in order to isolate the sensitivity of the model from spatial variability in the model forcings. The setup was forced with a uniform background tensile stress and an idealized runoff profile, illustrated in Figure S2 . Five key parameters were tested, which are discussed individually below. The sensitivity tests were run using a Monte Carlo method, where the likely range of each parameter was defined and a different randomly selected value from within this range was used for each model run. The probability distribution within each range was assumed to be uniform; however, as all of the parameter ranges are skewed to some extent, an equal number of samples were selected from above and below the base values.

There is some uncertainty in the selection of the crevasse width. Based upon the work of Krawczynski et al. [ 2009 ], a width of 1 m is a conservative estimate for the areas we are studying, where the ice is around 1000 m deep and under moderate tensile stress gradients. However, in order to address the uncertainty, sensitivity tests were carried out on a range of values of the crevasse width. Varying the crevasse length has a similar impact to the crevasse width as it effectively increases the volume of the crevasse.

In order to calculate the water depth b required in equation 1 , some assumptions were made about the geometry of the lake or crevasse being filled. We used the geometry illustrated in Figure 2 , where the water was allowed to fill a crevasse and any remaining water pooled on the surface. Similar to Clason et al. [ 2015 ], we used a crevasse width w c of 1 m, assumed to be uniform with depth, and we set the length equal to the grid size, in this case 40 m. The crevasse depth d was initialized with a depth of 1 m. During the first time step the model was spun‐up, allowing crevasses to penetrate to the depth driven by the background tensile stress.

Daily runoff data were gathered from the MARv3.5.2 model forced with ERA‐Interim Reanalysis data [ Fettweis et al. , 2013 ]. An example of the runoff data averaged over the catchment of W 3 is illustrated in Figure S2 in the supporting information . The runoff was routed to surface water areas identified from remote sensing data using a D‐Infinity flow routing algorithm [ Tarboton , 1997 ]. The flow was routed over a DEM which was pit filled everywhere except for locations where water was observed on the surface; thus, sinks could only form in these areas. Without this assumption we see a more uniform distribution of water on the ice surface and no cases of hydrofracture. Any runoff flowing into sinks within a surface water area was added to the volume of water within that area. This volume was then distributed evenly across grid squares within each surface water area. The catchments and approximated flow path lengths are summarized in Table 2 .

The resulting von Mises equivalent stress was used to represent the tensile stress R x x in equation 1 . The von Mises stress has been shown to be reliable for predicting the failure of glacier ice by Vaughan [ 1993 ] and has been widely used since [e.g., Hubbard and Hubbard , 2000 ; Clason et al. , 2012 ; Albrecht and Levermann , 2014 ]. For the purposes of the model, the stresses can be assumed to be constant through the depth of the ice [ van der Veen , 2007 ].

We ran the model on a 40 m resolution grid within a 24 km × 28 km domain, covering the terminus of Helheim Glacier. Within the domain, equation 1 was evaluated at each grid square with a time step of 1 day. The model was initialized with zero water depth at the start of each year. The model was spun‐up in the first time step, allowing crevasses to penetrate to the depth resulting from the background stress field.

When the tip stress, K I , reaches a critical fracture stress, K I C , the fracture begins to propagate downward. Equation 1 is then solved iteratively with d increased until K I < K I C or the crevasse has reached the bed. Following Clason et al. [ 2015 ], we used a K I C value of 150 kPa m 1/2 . The model was forced with surface stresses R x x , derived through feature tracking of synthetic aperture radar (SAR) imagery, and water depth b , determined using modeled runoff and an approximated crevasse geometry, both described in more detail below.

The high water line of the lake was identified in the imagery as a change in the ice surface from smooth, white ice within the lake basin, which we interpret as having been submerged, to rougher, darker ice, which characterizes the surrounding area and does not appear to have been submerged. The boundary roughly follows a line of constant elevation, and we therefore interpret this as the predrainage height of the lake in 2007, thus allowing the lake volume and depth to be derived.

Aerial photographs collected when the water level in lake L was relatively low were used to produce a digital elevation model (DEM) of the lake basin from which lake volume was derived. The photographs were captured on 24 July 2007 using a fully calibrated, aircraft‐mounted Wild RC‐10 Aviphot vertical aerial camera system and digitized with a high‐precision scanner in order to maintain radiometric and geometric fidelity. Photogrammetric processing was carried out in the SocetSET Photogrammetry Suite v.5.6 using ground control data extracted from a temporally coincident airborne lidar DEM, as described in James et al. [ 2006 ]. The photogrammetric adjustment yielded an RMSE in the adjusted ground control of 2.4 m in X and Y and 0.5 m in Z , which provides a good estimate of any systematic errors. An initial DEM of the lake and surrounding area was collected manually at a grid spacing of 50 m taking advantage of SocetSET's interactive 3‐D editing capabilities. This low‐resolution surface was used as a “seed” surface to constrain the automated terrain extraction in SocetSET's NGATE module. The resulting 5 m resolution DEM was then manually edited where required to produce the final DEM. Unlike most airborne lidar, which cannot penetrate even clear water, features on the submerged lake bottom were easily visible to the aerial camera through the clear water and thus were included in the final DEM without relying on interpolation.

In order to quantify the error in automatic classification, the area of lake L was manually digitized from 30 m resolution Landsat 7 and 8 images. These images were acquired from the U.S. Geological Survey using the LandsatLook Viewer, for the same day as the automatically classified MODIS images. From these images, 28 high‐quality, coincident image pairs were found between 2010 and 2014. Comparison of the areas from the automatic and manual classifications gave a root‐mean‐square error (RMSE) of 0.08 km 2 and an R 2 of 0.73, showing that the automatic classification performed reliably.

Water area was classified within pixel windows using an automatic algorithm. The method used here took the mean reflectance of the pixel window, and any pixels with reflectance values below a threshold of this mean were taken to contain water. The thresholds and window sizes used here are shown in Table 1 . This technique has been widely used in previous work [e.g., Box and Ski , 2007 ; Selmes et al. , 2011 ] and is reliable where there is a strong difference in reflectance between the ice and the water surface. For area L , the window contained an area of dark‐colored ice which led to a strong bias on the mean values from this window. Therefore, for this window only, we took the mean value from an adjacent control window which fell entirely on the light‐colored ice.

The lake, L , and areas containing water‐filled crevasses W 1, W 2, and W 3 (Figure 1 ) recur annually in the same positions. We therefore defined a number of pixel windows within the MODIS images centered on these areas, which we were able to use for all years. The sizes of the pixel windows are shown in Table 1 . If any pixels within the windows did not meet the filtering criteria, the entire window was discarded so as to avoid any contamination of the window. This strict filtering left between 62 and 104 high‐quality MODIS window images per year.

Atmospherically corrected, 250 m resolution MODIS Terra MOD09 Level 2 Surface Reflectance imagery [ MODIS Land Science Team , 2015 ] was used to automatically classify the presence of surface water at four locations. Images were acquired from the Level 1 and Atmosphere Archive and Distribution System Distributed Active Archive Center for the period 20 May to 30 September for all years from 2007 to 2014. Images were filtered using the “250 m Resolution Surface Reflectance Band Quality Description” and the “1 km Resolution Data State,” both produced during processing and supplied with imagery. Data were strictly filtered to exclude any pixels which were not classified as “highest quality” after processing. Pixels not identified as highest quality include those with an acute solar angle, noisy detector, or those that contain cloud or fall in the shadow of cloud.

This is a major simplification of a highly complex system, but it is illustrative of the rapid switching which can occur with changes in basal effective pressure. We also note that values of k will vary across the catchment, rather than the uniform values used here; therefore, the calculated flow routing indicates that in different years, and even within the same year, we would expect to see the flow routed in different ways within the catchment.

Two flow accumulation maps, calculated from the hydraulic potential, are also shown in Figure 6 . These are calculated using values for overburden fractions of 1.00 and 1.14. The full range of values 0.5 < k < 1.2 are presented in Movie S1 . Figure 6 and the Movie S1 both highlight the dramatic changes in flow switching which occur at different values of the overburden fraction k . At low values of k , the steep gradient of the hydraulic potential drives water into the central parts of the glacier. The water areas are located within a few hundred meters of the flow paths but may not be directly connected to the subglacial flow routing. As k increases, the hydraulic potential gradient shallows. Flow paths migrate toward the locations of surface water, with a switch connecting areas W 2 and W 3 around k = 0.97. Above k = 1, all areas of surface water appear to be connected but are not in a downstream order until a major switch between areas W 1 and W 2 around k = 1.13. A number of rapid switches in flow routing occur when k > 1; as the gradient of the hydraulic potential gets shallower, rapid flow switching appears more likely to occur.

The model consistently predicts the drainage of the lake L within 2 days of the observations, with the only outliers being 2009, which has been discussed previously, and 2010, which the model also slightly underestimates. The RMSE of the difference between the modeled results and observations for L is 4 days, showing that the model performs reliably for this area. However, the results for the lower water areas show much less consistency; RMSEs are 24, 13, and 15 days for areas W 1, W 2, and W 3, respectively. For W 1 all observed dates of hydrofracture lie above the 95th percentile of the results run with the standard parameter range. Extending the runoff parameter range captures all of the observed drainage times, but this implies that both runoff and crevasse width must be at the extremities of their parameter ranges to match the observations in a number of years. The results for W 2 and W 3 are slightly more consistent, but still two out of eight results are above the 95th percentile of the standard parameter range for W 2 and five out of eight for W 3. The strong relationship between modeled and observed results for L , which is not present in areas W 1, W 2, or W 3, is clearly illustrated in Figure 5 . While the difference in an individual year could be attributed to uncertainty in parameter choice, there is no consistency in this difference by area or by year. The lack of a relationship in the results would therefore require a different parameter selection for each year and area in order to see better agreement between modeled and observed results.

Modeled versus observed day of hydrofracture by area. Box plots are used to represent the range of the modeled hydrofracture results from sensitivity testing, where the whiskers represent the 5th and 95th percentiles, and the box represents the first, second, and third quartiles. Crosses represent the 95th percentile of the extended runoff sensitivity test. The black line represents a one‐to‐one relationship between the model and observations.

Figure 4 shows the 5 m resolution DEM produced from aerial photography. From the DEM, the maximum volume of L is estimated to be ∼9.7 × 10 6 m 3 in 2007. The moulin down which the lake L drained in 2007 is easily identifiable to the southeast of the lake. It is notable that in 2007 the lake is split by an ice divide, which is crossed by a narrow channel; while this would still allow the majority of water to drain from the lake, it may slow the rate of drainage.

The one exception to the down‐glacier progression of drainage is 2009, where we see slow growth of L with the maximum area occurring after W 1, W 2, and W 3 have drained. The slow filling rate cannot be explained by lower runoff, as we see no significant difference in runoff volume when compared to other years (Figure S2 ). We therefore identify two possible explanations for this: (i) surface flow routing is different in 2009, and runoff drains through a different connection to the bed rather than collecting in the lake, or (ii) hydrofracture creates a constricted connection to the bed allowing drainage at a rate less than the input of surface runoff. Both would result in suppressed lake growth but similar behavior to other years in W 1, W 2, and W 3.

Results from remote sensing and modeling of surface water areas. Filled diamonds represent water areas derived from satellite imagery, where the maximum and minimum water areas have been picked and used to scale diamonds in the vertical. Box plots are used to represent the range of the modeled hydrofracture results, where the whiskers represent the 5th and 95th percentiles, and the box represents the first, second, and third quartiles. Crosses represent the 95th percentile of the extended runoff sensitivity test. Colors of diamonds and box plots correspond to colors of water areas used in Figure 1

A summary of results from all years is shown in Figure 3 , and full results from the automatic classification of surface water for all years are included in the supporting information (Figures S3–S10 ). The filling and draining patterns of all water areas can be clearly identified, and variations in area are substantially larger than the errors associated with the method. While there is considerable annual variability in the maximum area of water, a number of clear patterns emerge. The most obvious is the consistent pattern of the drainage of L preceding the filling and draining of W 1, W 2, and W 3, all downstream and at lower elevations. Typically, L drains between the 20 and 30 June, while W 1, W 2, and W 3 generally drain in early July. The maximum area of L ranges between 0.25 and 0.52 km 2 , with the maximum area of 0.52 ± 0.08 km 2 observed in 2011. W 1, W 2, and W 3 have areas ranging between 0 and 0.48 km 2 , with W 2 and W 3 usually larger than W 1.

4 Discussion

The down‐glacier progression of drainage observed in our remote sensing data is unusual, and contrary to the results of Sundal et al. [2009] and Danielson and Sharp [2013] amongst others, who identified the order of drainage progressing upstream as the seasonal melt extent spread inland. The observations are particularly significant as the pattern is seen to occur in the majority of years between 2007 and 2014.

The observations could simply be explained by water draining over the surface of the glacier between water areas. However, we can clearly identify a moulin following the drainage of lake L from both Figure 4 and handheld photography collected in other years; this indicates that water is draining to either the subglacial or the englacial system rather than over the surface. While englacial drainage is a possibility, the highly stressed, heavily crevassed conditions at Helheim make it unlikely that water would flow for more than a few kilometers without encountering existing fractures or weaknesses allowing access to the bed. Additionally, the work of Andersen et al. [2011] suggests efficient drainage between the surface and the bed, implying that water is not retained in a complex englacial system.

We attempted to explain the order of drainage with variations in surface melt and catchment size using a LEFM model. The model correctly predicts the date of drainage of the lake L within a few days for the majority of years (RMSE =4 days); however, the results for the lower areas show much less consistency with the observations (RMSE =13 − 24 days). Sensitivity testing of the model shows that changes in two parameters, the runoff and the crevasse width, could account for the difference between observed and modeled results.

While the change in crevasse width could account for the difference between modeled and observed results in an individual year, it cannot account for the variability in the difference from year to year. Uncertainty in the crevasse width can largely be attributed to uncertainty in the shear modulus of ice [Krawczynski et al., 2009, Figure 1]. Any error in the estimation of the shear modulus would introduce a systematic bias to the results which would lead to a consistent overestimation or underestimation of the results. However, in order to explain the interannual variation in the difference between observed and modeled results, the shear modulus of the ice would have to vary by an order of magnitude from year to year. Variations in the tensile stress also affect the crevasse width, but on the order of 0.1 m per year [Krawczynski et al., 2009, Figure 1], which could not account for the interannual variability between observed and modeled results. It therefore seems unlikely that the uncertainty in crevasse width could explain the interannual variability of the difference between modeled and observed results.

Therefore, the only parameter which could be responsible for the interannual variation in hydrofracture appears to be the runoff. However, for runoff to account for the difference in the modeled and observed drainage times, catchment sizes would have to fluctuate in area by ±20–80% from year to year. The resulting changes in catchment size would lead to dramatically different patterns of surface water in different years. However, the observations show no such variation in the pattern of surface water, and, as drainage basins are typically tied to bedrock topography, they would not be expected to vary significantly [Karlstrom and Yang, 2016]. The variation in catchment size necessary to explain this variability is therefore difficult to reconcile with observations at Helheim, particularly the consistent locations and distribution of surface water. This evidence strongly suggests that hydrofracture is not the cause of the observed pattern of filling and draining at Helheim.

We therefore propose an alternate explanation which does not rely upon hydrofracture. We suggest that the down‐glacier order of filling and draining can be explained by a high‐pressure wave propagating down glacier following the lake drainage, controlling surface water levels as it passes. Transient high‐pressure waves such as this have been theorized and observed in association with jökulhlaups [Walder and Driedger, 1995; Tweed and Russell, 1999], producing pressures sufficient to flood the surface in areas of weakness [Tweed and Russell, 1999; Jóhannesson, 2009], and also during surges [e.g., Kamb et al., 1985; Fowler et al., 2001]. The propagation speed of the inferred pressure wave at Helheim is an order of magnitude lower than those observed in jökulhlaups; but drainage rates are similar to those observed in Antarctic subglacial lakes, where lower pressure gradients are thought to limit the rate of drainage [Fricker et al., 2007]. The system may therefore be more analogous to the sequential filling and draining of subglacial lakes observed in Antarctica [Fricker et al., 2016]. At the terminus of Helheim the ice is shallower and much more heavily crevassed than in Antarctica; therefore, where pressures greater than overburden occur, rather than raising the surface of the ice, this theory suggests that the water can penetrate to the surface and pool in crevasses.

Modeling studies have indicated that pressures above overburden can be sustained for between 4 days and 4 weeks and typically occur between late June and early July on the west coast of Greenland [Banwell et al., 2013; Dow et al., 2015]. Such pressures have been observed previously on the GrIS [Cowton et al., 2013; Meierbachtol et al., 2013] but rarely for periods of more than a few days. We suggest that the early season lake drainage transfers a large volume of meltwater to the bed at a rate sufficient to overcome the capacity of the existing subglacial system, thus forcing the propagation of a high‐pressure front down glacier as the system capacity increases. The observations and timescales are consist with the sustained high pressures observed by Banwell et al. [2013] in their modeling results. While Banwell et al. [2013] use a channel‐only model, the theorized response is also comparable to the hybrid channel and sheet model of Hewitt [2013], who also observe the subglacial system being overwhelmed at the beginning of the melt season. We do not have sufficient information to attempt to infer the morphology of the subglacial system present here; however, the models of Banwell et al. [2013] and Hewitt [2013] both suggest that the drainage of the lake L may trigger a change in the subglacial system which propagates down glacier.

Previous work on jökulhlaups has identified basal water flooding the surface; however, in Landsat imagery from 2015, the water visible in crevasses at Helheim has spectral characteristics more similar to isolated supraglacial lakes than to the turbid water seen in plumes and marginal water bodies. This suggests that while high basal water pressures may control water levels, the water visible on the surface is more likely surface melt which is prevented from draining through the crevasses by high basal water pressure. However, we note that Andersen et al. [2011] did observe turbid water upwelling into an open relict‐moulin structure just behind the calving front at Helheim, clearly showing that pressures sufficient for basal water to reach the surface can and do occur. The LEFM modeling work indicates a way by which a hydraulic connection between the surface and the bed could be created in the necessary areas during the early melt season. Alternatively, basal water pressures may be sufficient to fill or open existing fractures, as appears to have happened in the 1996 Grimsvötn jökulhlaup [Jóhannesson, 2009].

Further support for this hypothesis is provided by the hydraulic potential and flow routing results. Figure 6 shows that water‐filled crevasses appear in points of convergence and depressions in the hydraulic potential surface, particularly at higher values of k, which would be expected following a lake drainage. We acknowledge that the hydraulic potential at the bed is strongly influenced by surface slope, and it is therefore difficult to distinguish which is causing the water to collect in these areas. However, as we have highlighted previously, the surface of Helheim is heavily crevassed and we see no evidence for significant flow over the surface of the glacier which would cause water to collect in these areas. The positions of the water areas are therefore consistent with where water might be expected to collect if it were forced by conditions at the bed.

Figure 6 also shows dramatic variations in flow routing within the catchment when the basal water pressure is at different fractions of overburden, consistent with the work of Lindbäck et al. [2015]. As values of the overburden fraction k vary spatially and temporally within the catchment, the differences identified in both the order of filling and maximum extent of water areas can be explained by spatial and temporal variations of k. The flow switching in the catchment which connects all areas in a downstream order occurs within a few kilometers of the lake, within the region where uplift has been observed following previous lake drainages [Das et al., 2008; Doyle et al., 2013; Tedesco et al., 2013]. It is therefore quite probable that basal effective pressures reach the necessary values in this region for water to cross the catchment and connect to other areas of surface water downstream. In future, it may be possible to identify the flow routing regime and overburden pressures from more detailed study of these water areas, which could be used to complement other techniques such as boreholes and dye tracing. However, for the present we simply take this as evidence to explain the variability in timing of the filling and draining of the water‐filled crevasses.