High-resolution modelling of the seasonal evolution of surface water storage on the Greenland Ice Sheet

Seasonal meltwater lakes on the Greenland Ice Sheet form when surface runoff is temporarily trapped in surface topographic depressions. The development of such lakes affects both the surface energy balance and dynamics of the ice sheet. Although areal extents, depths and lifespan of lakes can be inferred from satellite imagery, such observational studies have a limited temporal resolution. Here, we adopt a modelling-based strategy to estimate the seasonal evolution of surface water storage for the ∼ 3600 km2 Paakitsoq region of W. Greenland. We use a high-resolution timedependent surface mass balance model to calculate surface melt, a supraglacial water routing model to calculate lake filling and a prescribed water-volume-based threshold to predict rapid lake drainage events. This threshold assumes that drainage will occur through a fracture if V = Fa · H , where V is lake volume,H is the local ice thickness and Fa is the potential fracture area. The model shows good agreement between modelled lake locations and volumes and those observed in nine Landsat 7 ETM images from 2001, 2002 and 2005. We use the model to investigate the lake water volume required to trigger drainage, and the impact that varying this threshold volume has on the proportion of meltwater that is stored in surface lakes and enters the subglacial drainage system. Model performance is maximised with values of Fa between 4000 and 7500 m 2. For these thresholds, lakes transiently store< 40 % of available meltwater at the beginning of the melt season, decreasing to ∼ 5 to 10 % by the middle of the melt season; over the course of a melt season, 40 to 50 % of total meltwater production enters the subglacial drainage system through moulins at the bottom of drained lakes.


Introduction
The formation of surface meltwater lakes on the surface of the Greenland Ice Sheet (GrIS) (and other Arctic ice masses) during the melt season is a widely observed phenomenon.However, many of the lakes observed early in the melt season have disappeared or decreased in size by later in the summer (McMillan et al., 2007;Selmes et al., 2011;Liang et al., 2012;Fitzpatrick et al., 2014;Morris et al., 2013).Whilst some lakes drain slowly by overflowing the lowest part of the lake rim, allowing water to escape from the depression by surface channel incision, other lakes drain rapidly by waterdriven fracture propagation (i.e."hydrofracture") ( Van der Veen, 2007;Das et al., 2008;Doyle et al., 2013;Tedesco et al., 2013).The sudden injection of large quantities of surface meltwater to the ice sheet bed and the associated reduction in basal friction has been proposed as a mechanism to explain observations of short-term increases in summer ice velocities following such drainage events (Zwally et al., 2002;Hoffman et al., 2011;Banwell et al., 2013).The storage (and subsequent release) of meltwater in supraglacial lakes is a key control in determining both the timing and rate of delivery of water to the subglacial drainage system.This is important as it is widely accepted that the variability in magnitude and timing of meltwater input to the subglacial drainage system has a greater influence on subglacial water pressures, and thus ice motion, than simply the total volume of water input (Schoof, 2010;Bartholomew et al., 2011Bartholomew et al., , 2012;;Colgan et al., 2011).However, despite numerous studies which compare ice velocity data to lake drainage observations from satellite imagery (e.g.Bartholomew et al., 2010Bartholomew et al., , 2011Bartholomew et al., , 2012;;Hoffman et al., 2011;Sole et al., 2011Sole et al., , 2013;;Sundal et al., 2011;Cowton et al., 2013;Joughin et al., 2013), there is still uncertainty about how supraglacial lake drainage events affect Published by Copernicus Publications on behalf of the European Geosciences Union.
basal water pressures, uplift and ice velocity over seasonal and annual timescales.Accurate estimates of potential water storage and the seasonal evolution of surface lake volumes on the GrIS will be an important component of ongoing work to understand how the ice sheet will respond to future climate scenarios.

Methods
Most studies of lake extent and dynamics have been based on a combination of direct surface observation and/or the use of multi-spectral satellite imagery (e.g.Box and Ski, 2007;Sneed and Hamilton, 2007;Sundal et al., 2009;Selmes et al., 2011;Liang et al., 2012).Whilst such studies have allowed the extent and the seasonal evolution of surface lakes to be ascertained for the times when such imagery is available, they are somewhat limited by the fact that lakes are transient features; not all lakes will be observed on any given image, and an image may not capture the maximum extent of any given lake.Recent studies using automated detection algorithms (e.g.Liang et al., 2012;Fitzpatrick et al., 2014;Morris et al., 2013) have alleviated these problems to a large extent, but such studies still face some limitations due to the availability of imagery of sufficient quality, and are perhaps best suited to tracking changes in total lake area over wide regions, rather than tracking the behaviour of individual lakes.
Model-based studies provide a valuable complement to observational studies.Although models require calibration and validation, they are not subject to issues such as the spatial and temporal resolution of available imagery, or gaps in coverage due to cloudy weather conditions.Models also allow the sensitivity of the predicted results to parameter and/or boundary conditions to be investigated, and can also be used for forecasting purposes.In this study, we adopt a model-based approach following Banwell et al. (2012b) and subsequently Banwell et al. (2013), by linking a highresolution surface mass balance (SMB) model (originally developed for Arctic glaciers (Rye et al., 2010), and subsequently used to calculate melt over the Paakitsoq region of the GrIS; Banwell et al., 2012a), to a new high-resolution surface routing and lake filling (SRLF) model.This SRLF model simulates supraglacial water flow by assuming Darcian flow for snow-covered portions of the surface, and open channel flow for bare ice, in order to calculate the filling rates of surface lakes over the course of a melt season.This in turn is linked to a water-volume-based threshold model of surface lake drainage (SLD) (Clason et al., 2012;Banwell et al., 2013).
Although high-resolution modelling of supraglacial hydrology is in its infancy, several previous studies have adopted a similar overall strategy for modelling the filling of supraglacial lakes by linking some form of surface melt model to a model for the movement of water over the ice surface.Luthje et al. (2006) link an ice melt model which allows for the absorption of solar radiation in lakes to a simple overland flow model which assumes that water velocity is proportional to the gradient of the surface slope, but do not model the difference between ice-and snow-covered surfaces; they use a tuning parameter to control the impact of slope on water flow to fit their model results to observations.Leeson et al. (2012) use daily runoff totals from the MAR regional climate model (Fettweis et al., 2011) as their water input, and then route water over the ice surface using the same combination of Darcian versus open-channel flow as we use.However, neither of these studies allow for surface lakes to drain through the ice sheet.
In this study, we run our combined glaciohydrological model (which we call "G-Hyd") for a ∼ 3600 km 2 area of the Paakitsoq region of West Greenland, just north of Jakobshavn Isbrae (Fig. 1).We use remotely sensed imagery to assess the performance of the model by first comparing modelled lake locations with observed lake locations in nine Landsat images from 2001, 2002 and 2005, and second by comparing predicted water depths in modelled lakes with water depths calculated from the imagery.
Our primary aims are threefold: (i) to evaluate the model performance over a larger area of the ice sheet than that assessed by Banwell et al. (2012b); (ii) to investigate the impact of altering the water volume threshold for lake drainage on model performance and the behaviour of the supraglacial hydrological system; and (iii) to consider the possible implications of the modelled behaviour on the subglacial hydrology of the ice sheet.

The G-Hyd model
The SMB model consists of three coupled components: (i) an energy balance component that calculates the energy exchange between the glacier surface and the atmosphere using measured meteorological variables; (ii) an accumulation routine to calculate winter accumulation, and possible summer snowfall; and (iii) a subsurface component, simulating changes in temperature, density and water content in the snow, firn and upper ice layers, and hence refreezing and net runoff of water (Rye et al., 2010;Banwell et al., 2012a).
The SMB model is forced by meteorological variables from the JAR1 GC-Net station, supplemented by data from the Swiss Camp GC-Net station when JAR1 data were unavailable (Fig. 1) (Steffen and Box, 2001;Banwell et al., 2012a).The main forcing variables are the incoming global short-wave radiation, air temperature, relative humidity and wind speed at a height of 2 m above the ice surface.Incoming long-wave radiation data were not available for the Paakitsoq region, so these data were calculated using parameterisations (Banwell et al., 2012a, following Konzelmann et al., 1994).The accumulation routine uses measured hourly precipitation from the Asiaq Greenland Survey station 437, ∼ 4 km west of the ice margin at an elevation of 190 m a.s.l.(Fig. 1).Precipitation is distributed over the ice sheet using an elevation-dependent precipitation gradient.Snowfall is calculated from this using a threshold temperature for snowfall of 2 • C (Rye et al., 2010;Banwell et al., 2012a).Output from the SMB model consists of hourly totals of available water in each grid cell, equivalent to surface melting on the ice surface for snow-free cells, or to percolating water entering the base of the snow/firn layer (effectively the total amount of melting, less any re-freezing within the snow/firn) for snow-covered grid cells.For the subsequent discussion we refer to this water as "runoff"; our previous work has shown that around 6 % of total melt and rainwater re-freezes within the snow pack in the study region (Banwell et al., 2012a).Full details of the SMB model, and the calibration and validation schemes for Paakitsoq, are given in Banwell et al. (2012a).
The SRLF model links an algorithm for calculating flow accumulation (upstream area) values over a digital elevation model (DEM) (Arnold, 2010) with a supraglacial water flow algorithm previously used on valley glaciers (Arnold et al., 1998) to calculate input hydrographs for all the depressions (and hence potential lakes) in the DEM.The supraglacial water flow algorithm assumes that water flows by Darcian flow in a saturated layer at the base of the snow pack in the case of snow-covered cells with greater than 0.275 m w.e. of snow, or as open-channel flow over shallower snow or bare ice surfaces.The snow depth for the transition from Darcian flow to open-channel flow was derived by calibration experiments reported in Banwell et al. (2012b), and allows for the formation of channels within shallow snow or slush, as are commonly observed on the ice sheet surface.From the DEM, the surface slope and water flow direction in each cell can be calculated; together with the inferred flow characteristics, this then allows the time taken for any "parcel" of water to cross each grid cell to be calculated.By integrating these times down the surface slope, a total delay time from "source" (initial runoff production) to "sink" (i.e. a lake/moulin, or exit from the model domain) for each hourly quantity of runoff in each grid cell can be calculated, allowing water input hydrographs to each sink cell to be derived.
The lake-filling algorithm assumes that any closed depression in the DEM acts as a sink for surface water, and may contain a lake.A catchment area which feeds meltwater into each depression can be calculated from the surface slopes and aspect of the DEM cells.The DEM also allows the lake hypsometry to be calculated, which, together with the input hydrographs from the supraglacial water flow component, allows lake depth and areal extent at any given time to be calculated.If a lake fills to its maximum possible extent (controlled by the elevation of the topographic "lip" which defines the depression), any further water inputs overflow into the downstream catchment.In this way, water can potentially flow from its source on the ice sheet in a series of "cascades" through multiple full lakes.Additionally, low points at the DEM boundary can act as sinks, and allow water to leave the model domain.Water could potentially enter our DEM domain at the margins which could be a possible source of error.However, this is minimised because the eastern (upstream) edge of our domain is at > 1500 m elevation where very little runoff occurs, and water flow is typically directed towards the margins at the northern and southern edges of our domain, as the ice surface beyond our domain is generally lower (Jakobshavn Isbrae is just to the south, Eqip Sermia to the north) than the ice surface within our domain.The western edge of our domain is ice-free terrain, and the model allows for runoff to exit freely at the edge of the ice sheet.
Full details of the model, and the calibration and validation of the assumed flow characteristics, are given in Banwell et al. (2012b), who applied this model to calculate water www.the-cryosphere.net/8/1149/2014/The Cryosphere, 8, 1149-1160, 2014 volumes in an instrumented lake within a 100 km 2 area of Paakitsoq in 2011 with a high degree of accuracy.
The SLD model uses a water-volume-based threshold to trigger lake drainage events.Fracture propagation and consequent lake drainage is assumed to occur if a lake reaches the volume needed to fill an inferred fracture extending from the ice surface to the bed (Clason et al., 2010;Banwell et al., 2013).The inferred width and length (i.e. the surface area) of the fracture can be varied in the model, and is assumed constant across the model domain.The threshold water volume needed is calculated by multiplying the fracture area by the local ice thickness.Thus, lakes over thicker ice have to reach a larger water volume in order to drain.We explore the model behaviour for inferred fracture areas (F a ) from 500 m 2 to 10 000 m 2 , and also with an infinite volume threshold, which prevents lake drainage events.If a lake reaches the volume threshold for drainage the lake empties, and any further water inputs are assumed to enter the subglacial drainage system directly through a moulin at the lowest part of the depression which is assumed to stay open for the remainder of the melt season (Banwell et al., 2013).
Following Banwell et al. (2013), we assume that water can only enter the subglacial drainage system via moulins which form at the bottom of drained lakes.Whilst this is undoubtedly a simplification, several studies have shown that moulins which occur away from depressions on the ice surface do not capture large volumes of water compared with those within depressions (e.g.Catania and Neumann, 2010;Hoffman et al., 2011).Over our study area, given this assumption, we calculate an overall moulin density, given the number of depressions, of 0.2 km −2 , which is very similar to the estimates given by Zwally et al. (2002) (0.2 km −2 ) and Colgan and Steffen (2009) (0-0.89km −2 ) for the Paakitsoq region.
For each of our study years (2001, 2002 and 2005) we begin with the DEM "empty" of water.The amount of runoff stored over the winter on the GrIS in refrozen lakes is largely unknown.However, in a study based on MODIS satellite imagery to the south of our region, Johansson et al. (2013) found that 78 % to 88 % of lakes they observed at elevations < 2500 m drained rather than refroze at the end of the melt seasons in 2007, 2008 and 2009.Given this, and given that refreezing is more likely at higher elevations where total runoff production is much smaller and lakes are smaller (Johansson et al., 2013), beginning with zero storage is the simplest and most justifiable initial DEM condition.Leeson et al. (2012) also make this assumption.

DEM processing
We use the GIMP product (Howat et al., 2014) for the highresolution surface DEM required by our modelling strategy, at the posted 90 m spatial resolution.Initial inspection of the data showed too many small (often single DEM-celled) depressions within the DEM.To combat this, we smoothed the data with a 2 × 2 cell median filter to remove much of the noise, and then applied an 11-cell radius Gaussian filter to remove the "terracing" effect that is a consequence of the 1m vertical resolution of the original data.
We are not aware of other research which has used this DEM for high-resolution surface modelling of melt and/or lake extent; the study by Leeson et al. (2012) used a DEM derived from 1996 InSAR data for their study region (Palmer et al., 2011).In some ways, then, this study acts as a type of validation of the small-scale (sub ∼ 1 km horizontal, ∼ 1-10 m vertical) resolution of the GIMP DEM.
Ice thickness data are needed to calculate the lake drainage threshold volume for each lake.We calculate ice thickness using the surface DEM and the 750 m resolution bed DEM described in Plummer et al. (2008), resampled using bilinear interpolation to 90 m resolution.

Analysis of lakes observed in satellite imagery
The modelled lake areas and volumes are compared with observed lake areas and volumes derived from nine Landsat 7 ETM+ images acquired in 2001 (7 July, 14 July, 1 August, and 8 August), 2002 (30 May, 17 June, 2 August and 3 September) and 2005 (16 June).Digital numbers were converted to reflectance values using standard methods (Chander et al., 2009).Following Box and Ski (2007), pixels with a Band 1 to Band 3 (blue to red) reflectance ratio above a certain threshold are considered to contain water.The threshold value used was chosen by comparing those pixels above the threshold by eye with true-colour projections of the imagery, in order to match as well as possible the apparent visible extent of water in the images, but to avoid as far as possible "false positive" identification of other blue pixels as containing water.Lakes are then identified as contiguous sets of water-containing pixels.The water depth in the pixels identified as containing water is then calculated using the method of Sneed and Hamilton (2007).This procedure uses the band 2 and band 4 reflectance data, but it also requires knowledge of the albedo of the bottom of each lake.We estimate the lake-bottom albedo on a lake-by-lake basis by creating a mask of each lake identified in the images, then dilating this mask by one pixel.The original lake mask is then removed, leaving a ring of pixels around the lake outline.The average reflectance of these pixels is assumed to be representative of the pixels within that lake.The volume of each lake is then calculated from the water depth of each pixel within the lake.For a fuller description of the application of the Box and Ski (2007) and Sneed and Hamilton (2007) method to the Paakitsoq region, see Banwell et al. (2014).

Model evaluation
We evaluate two aspects of model performance; the first focuses on the accuracy of the DEM in predicting potential lake locations; the second examines the predicted volumes of modelled lakes on the appropriate dates versus volumes calculated from the satellite imagery.Assessing the agreement between observed lakes in the visible imagery and DEM depressions is difficult for several reasons.Conventional classification methods based on contingency tables are biased because there will be many more locations that do not contain lakes or depressions than those that do, and also because whilst a lake can only form in a depression, it is possible for a depression to exist in the DEM but for it not to contain an observed lake; insufficient water may have accumulated in the depression for it to be observed, or the lake may have drained.Thus, an apparent disagreement between visible imagery and the DEM need not be due to an error.Given this, we assess co-location between observed lakes and depressions within the DEM by calculating the centroid of each lake identified within the visible imagery, and then determining the numbers of lake centroids within depression boundaries identified in the processed DEM versus the number of lake centroids outside depression boundaries.These values were converted into z scores by calculating the mean and standard deviation of the ratio for 1000 sets of the appropriate number of randomly distributed centroid locations within the study area.
We assess the modelled water depths by comparing them with the calculated observed water depths in the visible imagery at the appropriate dates using conventional regression techniques, and using the Nash-Sutcliffe model efficiency measure.The coefficients are calculated by taking the modelled water depth for all pixels containing water on the date of each image, and comparing them with the water depths as calculated from the Landsat image.Each image is treated as a distinct data set, but we amalgamate the resulting pairs of water depths (modelled vs. observed) from each date to calculate our overall coefficients of agreement to maximise the sample size.We also use the "weighted R 2 " value (Krause et al., 2005), calculated as where b is the slope of the regression relationship.Weighting R 2 in this way quantifies over-or under-prediction by the model as well as the overall level of agreement between modelled and observed depths (Krause et al., 2005).

Modelled lake location
Overall agreement between observed lake locations in the Landsat imagery and depressions within the DEM is good, and is summarised in Table 1; depression outlines are also shown in Fig. 1 for visual comparison with the image obtained on 7 July 2001.The flow accumulation algorithm of the SLRF model (Arnold, 2010) identified 644 depressions within the DEM, ranging in size from 1 pixel (8100 m 2 ) to 791 pixels (6.4 km 2 ).Using a blue/red reflectance ratio threshold of 2, a total of 505 separate lakes were identified (taking account of lakes which appear in more than one image) in the satellite imagery, with a maximum size of 332 pixels (2.7 km 2 ), of which 252 had centroids located within depressions.Using a blue/red ratio threshold of 3, a total of 229 lakes were identified, of which 179 had centroids within depressions, a success rate of 78 %.This compares favourably with the 66 % success rate achieved by Leeson et al. (2012).The largest lake for this threshold had an area of 277 pixels (2.2 km 2 ).As shown by the corresponding z scores for these numbers (Table 1), the chance of the coincidence in location between observed lakes and DEM depressions being due to random variation is vanishingly small.

Modelled water depths
Visual examination of the thresholded satellite imagery showed that a blue/red threshold of 2 tended to identify pale blue pixels some distance above the transient snow line as lakes.Sometimes these occurred in quasi-linear features that could perhaps be slush fields in shallow valleys on the ice surface; others, however, were isolated, very small clusters of pixels with no obvious topographic control.Thus, for the remainder of our analysis, we adopted a blue/red ratio threshold value of 3 to classify lake pixels, and from this determine water depths, and lake areas and volumes.Calculated water depth is not affected by this threshold value, however, and given that lower thresholds only increase the number of very shallow pixels around the edge of lakes (as well as the total number of apparently water-containing pixels with no obvious topographic control, especially at higher elevations), calculated lake volume is insensitive to this threshold.
Table 2 shows four measures of model performance as an estimator of observed water depths for the nine images across the three years.four images spanning early July to early August in 2001, and late May to early September in 2002, show similar results; the highest R 2 values occur for F a values of 2500 m 2 to 4000 m 2 , but the regression slope shows that the model generally under-predicts water depths for these F a values.As F a increases, the R 2 decreases marginally, but the weighted R 2 value and the Nash-Sutcliffe statistic increase until a F a of 7500 m 2 , after which they decrease again.The regression slope increases as F a increases, however, reaching a maximum some way above 1 for an infinite volume threshold, which effectively prevents any drainage, and leads to modelled water depths greatly exceeding their observed values.The year 2005, with just one image, from 16 June, shows different behaviour.Here, increasing F a from 1000 m 2 to 2500 m 2 results in an improvement in model performance, but as F a is increased more, there is no further change in model performance.Similar behaviour also occurs for the two early-season images in 2002 -no significant change in model performance occurs for F a higher than 2500 m 2 .

System behaviour for different drainage thresholds
Figure 2a shows the modelled evolution of total water storage in lakes in comparison with total modelled runoff for the three years, for F a = 5000 m 2 .The total volume of runoff produced in the three years is quite similar (3.6 − 3.8 × 10 9 m 3 ), but the seasonal evolution of runoff is quite different.Runoff begins earliest in 2005, stops, and then begins again, increasing rapidly to ∼ 30 June, when a cool period slows down the rate of increase until around 20 July, when runoff production increases again.The year 2002 shows a 1 Figure 2. a) Modelled evolution of total water storage (as a proportion of total ru 2 in comparison with total modelled runoff for the three years discussed in the 3 5000 m 2 .b) Cumulative lake volume (Cum.Lake Vol.), the cumulative amou 4 water (see text) (Resid.water), and the cumulative amount of runoff whi 5 subglacial drainage system via drained lakes (Subg.runoff), all as proport 6 cumulative runoff for 2001, for F a = 1000 m 2 ; 5000 m 2 and 10000 m 2 .7 (b) Cumulative lake volume (Cum.Lake Vol.), the cumulative amount of residual water (see text) (Resid.water), and the cumulative amount of runoff which enters the subglacial drainage system via drained lakes (Subg.runoff), all as proportions of total cumulative runoff for 2001, for F a = 1000 m 2 ; 5000 m 2 and 10 000 m 2 .
later, slower start in melt, but then a generally rapid rise in cumulative runoff for much of the season, to reach the highest overall total of the three years.In 2001, runoff starts even later, and rises only slowly until around 10 June, when the rate begins to increase to a similar level to 2002 and 2005.
Total lake volumes at the end of the melt season are also similar between the three years; around 5 % of the total runoff produced is stored in lakes at the end of the season.All three years show a peak in storage early in the season of around 35-40 % of runoff, with an initial rapid fall in 2001 and 2002 to around 15 % by 30 May, followed by a slow decrease (but with considerable short-term variability) throughout the rest of the melt season.The year 2005 behaves differently; volume as a proportion of runoff increases rapidly, but then decreases sharply once the short, early-season runoff period finishes.Once runoff begins again around 18 May, lake volume as a proportion of overall runoff increases until around 30 May, when it begins a similar downward trajectory to that observed in 2001 and 2002.
The cumulative lake volume, the cumulative volume of runoff which enters the subglacial drainage system via moulins in drained lakes, and the cumulative volume of the runoff which is not stored in lakes, and which does not enter the subglacial drainage system through the moulins in drained lakes for 2001 for F a values of 1000 m 2 , 5000 m 2 and 10 000 m 2 are shown in Fig. 2b.The latter portion could simply run off the surface of the ice sheet, or it could enter the en-or subglacial system via moulins that are not contained with lake basins, or via crevasse fields; it could also be stored within the ice sheets in crevasses which do not penetrate to the bed or other englacial voids.As stated above, we do not model these processes, and we dub this water "residual water" in the rest of the paper.Altering F a has a large impact on the relative proportions of the runoff which remains in storage on the ice, enters the subglacial drainage system, or remains as residual water.Very early in the melt season, all three F a values exhibit similar behaviour; around 35 % of runoff is stored in lakes, and around 60 % is residual water (the remainder is stored "in transit", within slow flow within the snowpack).
For F a = 1000 m 2 , the first drainage events occur within ∼ 3 days of the onset of runoff (15 May), leading to a rapid increase in the proportion of runoff entering the subglacial system to around 15 % by 18 May (∼ 1.8 × 10 6 m 3 out of a total cumulative runoff to date of 6.8 × 10 6 m 3 ), which then rises increasingly slowly to around 30 % by 10 June (4.2× 10 7 m 3 out of 1.5 × 10 8 m 3 ).Residual water remains at around 60 % (9.0 × 10 7 m 3 ) of cumulative runoff until this time.After 10 June, the proportion of runoff flowing as subglacial runoff begins to increase, steadily at first but then at a decreasing rate, reaching a final value of around 62 % of runoff (2.3 × 10 9 m 3 out of 3.8 × 10 9 m 3 ).Residual water as a proportion of total runoff begins to decrease after 10 June, quickly at first but at a gradually declining rate, reaching final proportions of ∼ 35 % of runoff (1.3 × 10 9 m 3 ).The proportion of water stored in supraglacial lakes decreases at a diminishing rate after the onset of drainage at around 18 May.For F a of 5000 m 2 and 10 000 m 2 , residual water continues to rise until around 30 May and 9 June, reaching peak values of ∼ 80 % (3.6 × 10 7 m 3 out of 4.5 × 10 7 m 3 ) and ∼ 85 % (9.0 × 10 7 m 3 out of 1.1 × 10 8 m 3 ), respectively.Subglacial drainage begins on 23 May and 12 June, respectively, leading to a decrease in the proportion of residual water, and an increase in the proportion of runoff entering the subglacial system.This increases steadily at first, and then at a decreasing rate, reaching end of summer values of ∼ 40 % (1.5 × 10 9 m 3 ) and ∼ 28 % (1.0 × 10 9 m 3 ) for 5000 m 2 and 10 000 m 2 , respectively.Storage within lakes shows a general gradual decline, but with episodes when the proportion varies by 3-5 % (∼ 3-9 × 10 6 m 3 ) over periods of a few days (particularly from ∼ 10 June to ∼ 20 July), with sudden drops which mark the drainage of individual large lakes.These individual drainage events are also visible in the time series of subglacial water volume.Final stored proportions are ∼ 5 % (1.5 × 10 8 m 3 ), and ∼ 8 % (2.6 × 10 8 m 3 ) of total runoff.Total residual water decreases, reaching final values of ∼ 54 % (2.0 × 10 9 m 3 ) and ∼ 62 % (2.4 × 10 9 m 3 ).
The upglacier progression of lake drainage events during 2001 for F a values of 1000 m 2 , 5000 m 2 and 10 000 m 2 is shown in Fig. 3a.Fewer events occur with higher F a values, but in all three cases there is a very clear upglacier progression of lake drainage events.The gradient of the curves around which drainage events cluster is steeper for higher values of F a ; the need for a higher volume of water leads to later lake drainage at any given elevation, as more time is needed for runoff to accumulate in the lake basins.Interestingly, clear clustering of drainage events occurs; for F a = 1000 m 2 , there are two distinct clusters of events around 28/29 June, for lakes with elevations of 900-1000 m, and a second around 14 July at ∼ 1100 m.Given the smaller total numbers of drainage events for higher F a values, clustering is less apparent but there is still a cluster of drainage events around 2 July for F a = 5000 m 2 for lakes with elevations around 950 m, and about 5 days later for F a = 10 000 m 2 .There is also a small cluster of events for F a = 5000 m around 24/25 July, at around 1100 m. Figure 3b shows the spatial distribution of lake drainage events for 2001 for F a = 5000 m 2 .Again, the upglacier progression of drainage is very clear; three pairs of lakes within ∼ 5 km of each other drain within 24 h of each other, and several other sets of neighbouring lakes drain within 2 days of each other (e.g. the pair of green lakes in the lower centre of the region, and three of the cluster of five yellow lakes just down-glacier from them).
Animations of the season-long filling and draining of surface lakes for 2001 for F a = 2500 m 2 , 5000 m 2 and 7500 m 2 are available in the Supplement (Fig. S1a, b and c).

Discussion
There are three controls on the volume of any given lake; the size and shape of the depression on the ice surface in which it forms; the volume of meltwater which has flowed into the depression by a particular time, and whether a drainage event has occurred (which previous studies have shown may be either complete or partial).Our finding of a significant coincidence in location between observed lakes and depressions in the DEM suggests that the DEM topography is accurate at the horizontal spatial scales used here, and also at vertical scales of 1-10 m.The good match between modelled and observed lake depths also supports this.The fact that the model performance early in the melt season does not change with F a values above 2500 m 2 suggests that the melt model is supplying accurate estimates of water inputs to the lakes, but that early in the melt season, the water volume drainage threshold plays little role in determining observed lake volumes as very few have collected sufficient meltwater to reach a drainage threshold; their volumes are effectively melt-limited.Later in the melt season, however, the drainage volume threshold plays an increasingly important role in determining which lakes drain, and when.Too low a threshold allows lakes to drain too quickly and at too small a volume; too large a threshold means that lakes drain too late, and reach too high a volume.Our findings suggest that F a values of around 4000 m 2 to 7500 m 2 times the local ice thickness produce the best match, although there is no particular value which emerges as the overall best-fit value, as the different measures of model performance suggest different optimum values.However, the fact that this behaviour occurs in both years with good seasonal image availability suggests that the water volume drainage threshold effectively simulates at least some of the physical processes that control lake drainage events.
The temporal complexity of behaviour exhibited by the model (Figs. 2 and 3) results from the interaction between the upglacier progression of runoff over the course of a melt season, the time taken for lakes to fill to the critical volume needed to trigger a drainage event (and the proportion of lakes which drain), and the distribution of lakes over the ice surface.Figure 4 shows the cumulative maximum water storage on the ice sheet surface with elevation.At lower elevations (below ∼ 700 m), water storage is small; lakes are relatively few in number, and small in volume.Between 700 m and ∼ 1200 m, water storage increases markedly, with the impact of several very large lakes clearly visible.Above ∼ 1200 m, the rate of storage increase with elevation declines, with the exception of one very large lake at ∼ 1450 m.
Thus, very early in the melt season, runoff is confined to the lowest elevations on the ice sheet, and largely consists of residual water which probably runs off supraglacially.Lakes close to the margin begin to fill with water.They are relatively small, however, and even for low values of F a , many do not drain; those which do reach their critical volume reach this volume, and therefore drain, within a few days.Drainage of these lakes therefore begins to limit the total proportion of runoff available as residual water for the runs with low F a values; at higher F a values, these lakes fill and then overtop their margins, giving a greater portion of residual water.For low F a values, the proportion of water draining subglacially begins to increase early in the melt season, but at higher F a values, water continues to accumulate on the surface, or forms residual water.
As the season progresses, the proportion of runoff stored within lakes decreases, as total lake volume is limited by the size of the depressions themselves, or by water entering the subglacial drainage system through drained lakes for low F a value runs.At some point, however, drainage events begin to occur at higher values of F a for the larger lakes found higher on the ice sheet.Thus, the proportion of subglacial runoff begins to increase in these runs; storage in lakes continues to decrease, but with periods when storage increases as large lakes at higher elevations begin to fill, and with "step" decreases marking the drainage of these larger lakes at their volume thresholds.Residual water begins to decrease as more runoff enters the subglacial system through lakes that have drained.The most active period of lake filling and draining occurs between around 10 June and the end of July; after this time, lake drainage events become fewer, but still occur into mid-August.Generally, most lakes which could drain -i.e.where the potential maximum water volume (determined by the topography) exceeds the volume threshold (determined by F a and the local ice thickness) -have drained by the middle of August; water volumes in the lakes at higher elevations are limited by the small amounts of runoff which occur at high elevation, so drainage events at higher elevations are rare.
The clear upglacier progression of lake drainage events (Fig. 4) is driven by two factors; the upglacier progression of melting at the onset of the melt season (linked to the general trend for lower runoff at higher elevations), which leads to later and smaller inputs of water to lakes at higher elevations, and the overall trend for thicker ice in interior parts of the ice sheet, which leads to higher water volume thresholds for higher-elevation lakes.The clustering of lake drainage events is interesting, however.Fitzpatrick et al. (2014) have observed clustering of drainage events in the Russell Glacier region and argue that this suggests that some form of synoptic trigger mechanism, based on the seismic and velocity response of the ice sheet to an "initial" event, could trigger additional, nearby drainage events, and that a lake-volumebased drainage trigger seems unlikely.However, our model results suggest that a lake volume threshold can lead to clustering of drainage events; lakes in a given area will receive broadly similar meltwater inputs (due to the relative spatial uniformity of runoff across a relatively flat ice sheet surface with quite uniform albedo and surface energy flux at scales of ∼ 1 km to ∼ 10 km), and will also have broadly similar ice depths beneath them, and so will have similar volume thresholds for drainage.Lakes with unusually early (or late) drainage dates (given their elevation) occur where the local ice thickness is lower (or higher) than the average for that elevation, or where the supraglacial catchment feeding the lake is unusually large (or small).
Of course, our results do not rule out the "triggering" of drainage events by other nearby events, but they show that clustering of drainage events does not rule out a volumebased trigger for the drainage of individual lakes either.In many ways, some combination of both effects would seem most likely; a lake would most probably need some threshold volume of water in order to create a hydrofracture to the bed, but the "moment" of drainage could be triggered by flexural stresses associated with a nearby drainage event once the lake was "primed" with sufficient water.The clear upglacier progression of drainage observed in all of the remote sensing-based studies, and the need for at least one lake in a region to drain "spontaneously" in order to trigger other drainage events, also seems to support some form of water volume threshold as a condition for lake drainage.

Conclusions
In this study, we have applied a linked surface mass balance/supraglacial water flow/lake drainage model to a 3600 km 2 region of the Greenland Ice Sheet (GrIS).
Our research suggests that the current generation of highresolution DEMs available for the GrIS are sufficiently accurate to enable the prediction of potential supraglacial lake positions and extents.In particular, it seems unlikely that the data used to derive the GIMP DEM systematically includes "full" lake surface elevations within it, as the depressions within the DEM produce realistic water depths when compared with remotely sensed imagery.This is almost certainly a result of the wide range of data sources used in the compilation of the DEM (Howat et al., 2014), and suggests that this broad-based approach is very effective at capturing small-scale topographic variations on the ice sheet surface, certainly at elevations up to around 1500 m a.s.l.
The overall pattern of behaviour produced by the model matches that found in other observation-based studies, and matches well with observed water depths in our study area.This suggests that our surface mass balance (SMB) model and surface routing and lake filling (SRLF) model performs sufficiently well to provide realistic water inputs to lakes over a wider area of the ice sheet than has been investigated in previous studies (Banwell et al., 2012a(Banwell et al., , b, 2013)).
When linked to a water volume threshold-based surface lake drainage (SLD) model, the combined glaciohydrological model (G-Hyd) is capable of producing good agreement between modelled lake extents and observed lake extents from remotely sensed imagery during three melt seasons (2001, 2005 and 2005), and reproduces the overall upglacier trend in lake filling and drainage commonly observed.The model results suggest that the water volume threshold needed to trigger lake drainage is somewhat higher than that suggested or used by previous studies that have used volumebased thresholds (e.g.Clason et al., 2012;Banwell et al., 2013).We find the best model performance for inferred fracture areas of between approximately 4000 m 2 and 7500 m 2 times the local ice thickness.However, no clear "best fit" value emerges, as different measures of model performance produce different optimal values.
The model results indicate that whilst a linear relationship between the water volume threshold and ice thickness is undoubtedly a simplification, the volume threshold needed www.the-cryosphere.net/8/1149/2014/The Cryosphere, 8, 1149-1160, 2014 for drainage must vary with elevation (and by inference, ice thickness).Early model experiments with fixed drainage threshold across the model domain showed that low thresholds were needed to allow for drainage of the small lakes near the ice margin, but these low thresholds prevented lakes higher on the ice sheet reaching the areas and depths commonly observed in satellite imagery, as these lakes would drain much too quickly.By contrast, a high threshold (which allowed these lakes to reach observed areas and depths) effectively prevented small lakes at lower elevations from draining.
The model also produces spatial and temporal clustering of drainage events, as has been observed; unlike the proposed synoptic triggering of Fitzpatrick et al. (2014), the clustering produced by this study is purely the result of nearby lakes experiencing similar filling rates (due to similar local runoff totals), and having similar water volume thresholds due to similar ice thicknesses in the area of the lakes.This study does not rule out some form of synoptic triggering, but it suggests that clustering can occur simply due to some overall influence of ice thickness on the water volume needed to force drainage to the bed of the ice sheet.
The water volume threshold acts as a primary control on the overall proportion of runoff which is stored on the ice surface versus that which enters the subglacial drainage system.Low volume thresholds lead to smaller amounts of supraglacial storage in lakes, and larger amounts of runoff entering the subglacial system compared with larger volume thresholds.For the best fit values of inferred fracture areas of between 4000 m 2 and 7500 m 2 , lakes transiently store < 40 % of runoff at the beginning of the melt season, but this decreases to ∼ 5 to 10 % by the middle of the melt season.40 and 50 % enters the subglacial drainage system through the bottom of drained lakes, and a similar amount remains on the surface as supraglacial water, which could enter other forms of supra-or englacial storage (e.g.waterfilled crevasses) or which could enter the subglacial drainage system via moulins outside of lake basins, or crevasses which penetrate to the bed.Near the margin, this water could simply run off supraglacially.
Numerous recent studies have highlighted the complexities in the relationship between supraglacial (and subglacial) hydrology, and surface velocity for areas of the GrIS (e.g.Bartholomew et al., 2010Bartholomew et al., , 2011Bartholomew et al., , 2012;;Hoffman et al., 2011;Sole et al., 2011Sole et al., , 2013;;Sundal et al., 2011;Cowton et al., 2013;Joughin et al., 2013).In particular, the recent study by Joughin et al. (2013) highlights the fact that the timing of lake drainage events controls when much of the observed seasonal speedup occurs; they also argue that it is the combination of surface and bed slopes which act together to determine water flow directions, and the consequent patterns of velocity.The G-Hyd model reported here, used in conjunction with recent high-quality, high spatial resolution surface topography data, is capable of predicting the locations, volumes and timings of water inputs to the subglacial drainage system with a good degree of accuracy.In future, models of this type will be able to provide the key surface water input estimates which will subsequently permit more effective modelling of the impact of subglacial hydrology on ice dynamics.
The Supplement related to this article is available online at doi:10.5194/tc-8-1149-2014-supplement.

1Figure 1 .Figure 1 .
Figure 1.Location map of the study area.Detailed map shows Landsat true-colour image for 2 7 July 2001.Coordinates are UTM, zone 22W.Red lines show maximum possible lake 3 extents calculated from the DEM using the method of Arnold (2010).Black rectangle in inset 4 map shows study area location within Greenland.5 6

Figure 2 .
Figure 2. (a)Modelled evolution of total water storage (as a proportion of total runoff) in lakes in comparison with total modelled runoff for the three years discussed in the text, for F a = 5000 m 2 .(b) Cumulative lake volume (Cum.Lake Vol.), the cumulative amount of residual water (see text) (Resid.water), and the cumulative amount of runoff which enters the subglacial drainage system via drained lakes (Subg.runoff), all as proportions of total cumulative runoff for 2001, for F a = 1000 m 2 ; 5000 m 2 and 10 000 m 2 .

Figure 4 .Figure 4 .
Figure 4. Cumulative maximum potential surface water storage in surface depressions by 2 elevation for the Paakitsoq region.3

Table 1 .
Co-location statistics for observed lakes and DEM depressions.

Table 2 .
Model performance for different fracture area thresholds (F a , m 2 ).