Simulating the growth of supraglacial lakes at the western margin of the Greenland ice sheet

We present a new method of modelling the growth of supraglacial lakes at the western margin of the Greenland ice sheet, based on routing runoff estimated by a regional climate model across a digital elevation model (DEM) of the ice sheet surface. Using data acquired during the 2003 melt season, we demonstrate that the model is 19 times more likely to correctly predict the presence (or absence) of lakes than it is to make incorrect predictions, within an elevation range of 1100 to 1700 metres above sea level (m a.s.l.), when compared with MODIS satellite imagery. Of the 66 % of observed lake locations which the model correctly reproduces, the simulated lake onset day is found to be correlated with that observed with a Pearson correlation coefficient of 0.76. Our model accurately simulates maximum cumulative lake area with only a 1.5 % overestimate. However, because our model does not simulate processes leading to lake stagnation or decay, such as refreezing or drainage, at present we do not simulate absolute daily lake area. We find that the maximum potential lake-covered ice sheet area is limited by topography to 6.4 %. We estimate that this corresponds to a volume of 1.49 km3, 12 % of the runoff produced in 2003. This can be taken as an upper bound given uncertainty in the DEM. This study has proved a good first step towards capturing the variability of supraglacial lake evolution with a numerical model. These initial results are promising and suggest that the model is a useful tool for use in analysing the behaviour of supraglacial lakes on the Greenland ice sheet in the present day and potentially beyond.


Introduction
Supraglacial lakes (SGLs) have been observed to form during the summer melt season across much of the ablation zone of the Greenland ice sheet (GrIS) (McMillan et al., 2007;Sundal et al., 2009;Sneed and Hamilton, 2007;Selmes et al., 2011).Observations show that they form in the same locations each year (Echelmeyer et al., 1991), in depressions that are controlled by the underlying bedrock topography and by spatial variations in the degree of basal ice lubrication (Gudmundsson, 2003).Runoff, a combination of melted surface snow and ice, and wet precipitation, flows over the ice and through firn to areas of lower hydraulic potential.Prior to being routed off the ice sheet and into the surrounding ocean, this water may collect in SGLs.Model studies have suggested that the area covered by SGLs is controlled by surface topography (Luthje et al., 2006).
Supraglacial lakes impact upon the mass balance of the ice sheet in several different ways; melt ponds have been shown to reduce the albedo of large areas of ice (Perovich et al., 2002), thereby promoting additional melting.They are also temporary water storage sites, which can modify the rate at which runoff leaves the ice sheet.At later stages of the melt season, a proportion of SGLs drain rapidly (Selmes et al., 2011) as a result of hydro-fracture, presumably once a critical lake volume has been reached (Krawczynski et al., 2009;van der Veen, 2007), while others drain continuously at slow rates.Rapid SGL drainage, which can occur over periods as short as a few hours, has been observed to contribute to short-term increases in the flow velocity of the ice sheet (Das et al., 2008).This rapid drainage may also precede seasonal increases in the rate of ice flow (Shepherd et al., 2009).It is believed that draining lakes provide a mechanism by which this occurs through opening up conduits to the base of the ice sheet and delivering large amounts of water to lubricate flow (Joughin et al., 2008;Zwally et al., 2002).This leads to a third mechanism through which ice sheet mass balance may be altered, because changes in the flow of ice can lead to changes in ablation, through modified hypsometry.Observations and model studies have shown that an increase in meltwater supply to the ice bed results in a more efficient subglacial drainage system, which can offset this speedup (Schoof, 2010;Sundal et al., 2011).However, short-term spikes in water supply, such as from the drainage of SGLs, can lead to short-term, high-magnitude velocity increases even after this efficient drainage system has been established (Sole et al., 2011).It has been suggested that in a warming climate, these short-term velocity variations will propagate further inland due to the increased abundance of meltwater (Sole et al., 2011;Bartholomew et al., 2011), which could result in a net acceleration.
Supraglacial lakes on the surface of the GrIS have been surveyed using a variety of satellite imagery (e.g.Selmes et al., 2011;Sneed and Hamilton, 2007).These satellite data are, however, sparsely distributed in time relative to the period over which lakes are typically present.Knowledge of SGL lake evolution over shorter time periods would be of benefit for pinpointing the timing of and conditions required for lake drainage.
Here, we present a new model, which simulates seasonal SGL evolution at the western margin of the GrIS.The model combines a digital elevation model (DEM) with primitive hydraulic flow equations to route and pond surface water runoff estimated using the MAR (Modèle Atmosphérique Régional) regional climate model (RCM).We discuss the model performance in so far as it is able to predict the location of and describe the growth of SGLs in 2003.Finally, we extend the study to incorporate later melt seasons (2005 to 2007) to assess the model's capacity to represent inter-annual variability in SGL lake evolution.

Study area and data
The model domain is restricted to a 6753-km 2 region of western Greenland for which fine spatial resolution elevation data are available (Palmer et al., 2011), lying above 1100 m above sea level (m a.s.l.) (Fig. 1).This incorporates much of the area of the Russell glacier catchment.In this region, the ice sheet elevation ranges from 1100 to 1752 m a.s.l.This entire region experiences temperatures above freezing during the ablation season, leading to abundant melting and runoff.
We use a DEM derived from interferometric synthetic aperture radar (InSAR) data acquired by the European Remote Sensing satellites in 1996 (Palmer et al., 2011) as the basis of our scheme for routing water across the ice sheet surface.The DEM uses the Universal Transverse Mercator (UTM) projection system and is posted at a grid resolution of 100 m.We use this as our reference projection and translate all other datasets into the same coordinate system.We geolocate the DEM using the UTM and a DEM covering the whole of Greenland (Bamber et al., 2001) to ensure accurate spatial referencing.Based on a comparison with airborne laser altimeter data, it has been estimated that DEMs formed from repeat-pass InSAR data achieve a relative accuracy of between 2.5 and 10.0 m, depending upon the length of the perpendicular baseline of the InSAR data (Joughin et al., 1996).The DEM used here was formed with a perpendicular baseline of 120 m, and we estimate the relative accuracy to be 10.0 m.Although radar elevation measurements correspond to horizons at depth relative to the ice sheet surface (due to penetration of the microwave signal into the snowpack), a study of InSAR-derived topography has demonstrated that these horizons follow the ice sheet surface (Rignot et al., 2001).This is because topography is strongly correlated with basal conditions that are transmitted through the ice (Gudmundsson, 2003), and so we are able to use this product to approximate surface topography.Large topographic gradients exist in the DEM at the data coverage margin, which arise as an artefact of smoothing performed during DEM derivation.We minimise the impact of these gradients on our model by removing all cells which exhibit gradients that differ by more than one standard deviation from the local mean.This removes a margin of ∼ 1 to 3 cells around the ice edge only.
We assume that the ice sheet DEM is "empty" at the start of the melt season, i.e. that local depressions in the surface contain no water.We know from observations that lakes which have not drained freeze over at the end of the melt season (Sundal et al., 2009;Selmes et al., 2011).If lakes do not freeze completely, then we can assume that there is an ice "lid" covering the lake.The DEM used in this study was formed from InSAR data collected in the winter of 1996.Rignot et al. (2001) showed that for ice, specifically the surface of Jakobshavn Isbrae, which is 193 km further north than the centre of our study area, the InSAR penetration can be as shallow as 0 m, typically 1 m (±2 m).For firn this depth was found to be closer to 10 m, however, much of the DEM (89 %) lies below the permanent snowline (∼ 1600 m a.s.l. in MAR).If the ice lid is thick (1 to 3 m), the surface of the lake and the surrounding ice in the DEM will be slightly below their absolute values, but the gradient between them will be the same.If the ice lid is thin, the radar is reflected at the lid/lake boundary; radar cannot penetrate water, and the gradient between lake/surrounding ice will be artificially shallow by up to 3 m.
We use estimates of surface runoff derived from a daily product of the MAR RCM (Fettweis, 2007;Fettweis et al., 2011) run at 25-km resolution for Greenland and forced at its boundaries by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis.For this study, we principally use runoff estimates produced for the year 2003, as this is the year for which the most complete time series of satellite observations of lakes in this area is available for comparison.MAR is a dynamic downscaling RCM, which takes a limited geographical area at fine spatial and temporal resolution and embeds it in a global General Circulation Model (GCM) on a coarse grid/timescale.Atmospheric conditions such as temperature and pressure at the boundary of the RCM are prescribed every 6 model hours by the GCM, and additional physics are employed to downscale these data to a finer spatial and temporal resolution.The runoff product provided by the MAR model has been calculated using a snow and ice melt model which accounts for retention, percolation and refreezing of meltwater (Lefebre et al., 2003).Although runoff from MAR has not been explicitly validated due to lack of observations, the annual MARsimulated surface mass balance (precipitation-runoff) has been successfully compared to independent observations acquired over the period 1990 to 2011 along a transect located in the study area (Tedesco et al., 2011).In addition, Fettweis et al. (2011) have shown that the daily melt extent simulated by MAR compares well with that derived from satellite data over the period 1979 to 2009.The margin of the ice sheet experiences high climatic variability, primarily due to the steep elevation gradient but also due to the heterogeneity of surface albedo relative to the interior (Wientjes, 2011;van den Broeke et al., 2008).Studies of the impact of spatial resolution on the ability of MAR to resolve surface mass balance (SMB) gradients, such as temperature and runoff, close to the ice sheet margin have shown that the 25-km resolution product is best used further inland (Franco et al., 2012).The lowest extent of our study area lies ∼ 25 km, i.e. one MAR grid box, away from the margin and represents the lowest point from which the MAR data can be said to accurately resolve SMB.The MAR data are re-projected and oversampled to the model domain before being supplied to the surface for routing.
We use observations of supraglacial lake area derived from MODIS satellite imagery to assess the skill of the model in predicting their locations and area.Observed changes in the area of 492 SGLs were determined using an automated classification of cloud-free MODIS images acquired on 28 separate days during the 2003 summer melt season and 13 separate days during the melt seasons of 2005 to 2007 (Sundal et al., 2009).The spatial resolution of the MODIS instrument is 250 m, and SGLs of smaller area than this are therefore not resolved.It is estimated that total SGL area is underestimated by 12 % when using the relatively coarseresolution MODIS data, as compared to a manual classification of fine-resolution (15 m) ASTER imagery (Sundal et al., 2009).There is also an error associated with the miscategorisation of ice-covered lakes as bare ice using this automated classification method (Sundal, 2009).The data used in this study have been manually corrected for this error.To facilitate the comparison between simulated and observed lakes, we re-projected the distribution of observed lakes on each day into the model domain.
MODIS observations are not available for 1996, the year during which the InSAR data from which the DEM was formed were acquired.Since surface topography is strongly controlled by basal conditions (Gudmundsson, 2003), we assume a low degree of inter-annual variability in the surface topography and use this DEM to represent the ice sheet surface in the period 2003 to 2007.This assumption is supported by the findings of Echelmeyer et al. (1991) and Selmes et al. (2011), who observed that SGLs occur in the same locations year on year.

Method
We employ a fully transient 2-D hydrology model based on Manning's equation for open channel flow (Manning, 1891) and Darcy's law for flow through a porous medium (Chow, 1988) to simulate the evolution of SGLs on the surface of the GrIS.The study area is represented on a model grid of square cells of 100-m length.The model first calculates flow direction on the surface and subsequently calculates water displacement between cells for a given time step.Lake depth is calculated using an iterative algorithm at the end of each time step.In the standard version of the model, the time step is 90 s.
The direction of water flow within each model grid cell is calculated at the beginning of each time step, and is defined A. A. Leeson et al.: Simulating the growth of supraglacial lakes as being towards the neighbouring cell whose free surface (ice surface plus lake plus free water) elevation is lowest with respect to the reference cell.If the cell itself is the lowestlying, then it is designated a "sink" to be potentially filled with water.
Water displacement is calculated on a cell-wise basis.In order to estimate the change of water volume in any given cell, flux out of the cell is combined with the rate of runoff production within the cell and the resulting equation is integrated over a given time step.Volumetric flow rate, or flux, Q, can be calculated using Darcy's law for flow through a porous medium (Chow, 1988) where snow is present in the cell: Here, k is the permeability of the snow determined using Shimizu's equation (Shimizu, 1969) and snow density calculated by MAR.
A is the cross-sectional area of flow (Fig. 2) and µ is the viscosity of water, taken here to be 1.763 × 10 −3 Pa s.Where the cell is bare ice, Q can be calculated using Manning's equation for open channel flow (Manning, 1891): Here, v is the cross-sectional average flow velocity, n is "Manning's n", an experimentally derived roughness coefficient accounting for channel friction, and R h is the hydraulic radius of the "channel", which can be approximated by depth of flow for wide, shallow channels.In this study, n takes the value 0.011, the midpoint of a range of n values (0.01 to 0.012) derived experimentally for ice by Lotter et al. (1932).In both cases, S is the free surface slope, which is calculated by dividing the hydraulic head by the path length of flow, P (Fig. 2).The hydraulic head is calculated by: where z i and z j are the ice/lake surface elevation and d i and d j are the depth of runoff in cells i and j , respectively.Path length is L for flow between cells which share a border and √ 2L 2 for cells with a diagonal path.Rate of change of depth of water in the cell ( ḋ -we use dot notation to denote the time derivative) can be found by stating that the volumetric flow rate is equal to the rate of change of depth (d) due to flow out of the cell multiplied by the surface area of the cell L 2 : There is also likely to be a change in depth due to runoff production within the cell.Runoff produced within the cell through melting or precipitation is provided by MAR in mm Fig. 2. Two-dimensional schematic of model processes including a sink cell (j ).(Top) state of model at time t = 0, a layer of MARsimulated runoff (blue) covers each cell.(Bottom) state of model at time t = t (one standard time step); the sink cell (j ) now has a lake surface (teal), given by Eq. ( 7), and cells i and k contain a combination of runoff from contributing cells and that internally generated by MAR. per day.The sum of these contributors to rate of change in depth (in seconds) is therefore given by: Because Q is proportional to surface slope, which is described as a function of two independent variables, d i and d j , an additional equation is required in order to integrate Eq. ( 5) with respect to time.Equation ( 5) can be rewritten for d j since ḋFLUX This gives a system of coupled differential equations, which in this study are integrated with respect to a given time step using a fourth-order Runge-Kutta approximation.
It is possible that one cell receives water from many contributing cells.The new depth of water in any given cell, n, at time t is given by: Equation ( 6) adds the summation of water passed into the cell from its neighbours (if any) and the amount generated by MAR within the cell to the depth of water at time 0 and then subtracts the displaced water due to flow.A condition is imposed where if the calculated amount of displaced water is more than the sum of the depth of water in the original cell at time 0 and the water generated within the cell by MAR, only the original cell content plus MAR runoff is displaced.
At the end of each time step, after the flow from each cell has been calculated, an iterative accumulation algorithm brings water together into lakes.In a similar manner to the flow direction algorithm, each cell is considered in turn and if a cell is lowest-lying with respect to its nearest neighbours, it is designated a sink cell.If a cell is found to be a sink and the height of the cell plus water is less than the height of its lowest-lying nearest neighbour, all the water in the cell is incorporated into the lake/ice surface DEM.If the height of the cell plus water is greater than the height of its lowestlying neighbour, the difference between the heights of the two cells, plus 1 mm, is incorporated.The algorithm loops round the cells repeatedly until all potential sinks are filled.Once all possible lake water has been accumulated, the lake depth is calculated by subtracting the baseline ice surface elevation from the new lake surface, which is equal to the height of the ice surface plus that of the lake.This process is illustrated in Fig. 2. At time t = 0 the three cells all have water present, cell j is a sink and cells i and k flow into j .The new depth at time t is given by Eq. ( 7), but because cell j is a sink and the new depth is less than the free surface height of its nearest neighbour, all of this water is incorporated into the new lake surface.
We do not model the rapid drainage of lakes, although slow lateral seepage and overspill are inherently included in this method.Knowledge of the conditions required to promote rapid drainage is limited (Selmes et al., 2011) and so it is not appropriate to attempt to parameterise it at this time.
We assessed the model sensitivity to runoff amount, because the runoff produced by MAR has not been explicitly validated.To do this, we ran the model with runoff equal to fixed fractions of that estimated by MAR.We compared the model results for each of these sensitivity experiments to observations of lake location, area and filling rates determined from the MODIS satellite dataset.

Results
We assess the skill of the model in predicting lake location on a lake-by-lake and on a cell-wise basis across our study area.Figure 3a shows a comparison between observed (red) and simulated (blue) lake distribution, with coincident lake area highlighted in purple, in our study area.There are simulated and observed lakes that are in close proximity, but that are not coincident.There are also a number of incidences of one large observed lake coinciding with several smaller modelled lakes.However, in general we see good agreement between the model and observations in terms of locating lakes; we find that we correctly simulate the location of 66 % of observed lakes greater than 0.125 km 2 in our study area.
In order to quantify model skill in predicting the presence or absence of lake on a cell basis, we use the "odds ratio", the Heidke skill score (HSS), the Peirce skill score (PSS) and also the Threat score.We refer to Stephenson (2000) for a comprehensive description of each method.The number of simulated lake cell hits (predicted and observed), false alarms (predicted but not observed), misses (observed but not forecast) and correct rejections (neither observed nor forecast) for 2003, the year during which the sampling of the satellite observations is most dense, are presented in the form of a contingency table (Table 1).Using this contingency table we calculate the odds ratio to be 18.93, indicating that using our model we are ∼ 19 times more likely to predict a cell to be part of a lake or not than to make an incorrect prediction.The odds ratio can be positively skewed by a low event-non-event ratio.These data presented here apparently exhibit this rare event bias, given that there are many more cells which do not form part of a lake than there are cells which do.However, we consider a correct rejection (correctly predicting a cell not to contain a lake) to be a positive result of successful water routing and ponding in this case.The Heidke skill score is found to be 0.32, indicating that the model performs 32 % better than a random assignation of lake/non-lake cells over the study area.The Peirce skill score is found to be 0.33 and the Threat score is found to be 0.21, both of which indicate a useful degree of skill.
We investigated the extent to which the model run at 90s resolution was able to reproduce the temporal evolution of SGLs, both in terms of onset day and subsequent growth.The www.the-cryosphere.net/6/1077/2012/The Cryosphere, 6, 1077-1086, 2012  simulated onset days of the 66 % of observed lakes successfully placed by the model (159 of 242 lakes) are found to be correlated with the observed onset days of these lakes with a Pearson correlation coefficient of 0.76 (Fig. 4).A linear fit through these data suggests that observed lakes appear up to 5 days earlier than modelled lakes in the early part of the melt season.Later in the melt season, some lakes are simulated before they are observed.We calculate the mean lead in the observations to be 4 days when considering all 159 lakes.Lakes initially form at lower elevations and subsequently form progressively further inland over the course of the melt season (Fig. 5).Lakes first appear between 1100 to 1300 m a.s.l., then between 1300 to 1500 m a.s.l., and lastly between 1500 to 1700 m a.s.l.This pattern can also be seen in the MODIS dataset presented here (Fig. 5).Simulated lake initiation is contemporaneous with observations in the upper two elevation bands at days 155 and 162 for 1300 to 1500 m a.s.l. and 1500 to 1700 m a.s.l., respectively.However, observed lakes begin to appear 4 days earlier than simulated for the 1300 to 1500 m a.s.l.elevation band.Since the processes of lake drainage and refreezing are not incorporated within our model at this time, our results are restricted to comparisons of lake filling and are meaningful only until the date of drainage/commencement of refreezing.However, we can investigate model skill in simulating cumulative lake area, assuming no drainage, by comparing model output for daily lake area with an estimate of cumulative lake area derived from the satellite observations (Fig. 5).The model predicts the maximum cumulative area of lakes very well, with an overestimation of just 1.5 % of the maximum lake area between 1100 and 1700 m a.s.l.The maximum cumulative area occupied by lakes is reached around day 203 (21 July) in the observations -about one week later than the maximum daily runoff amount in the MAR data, which occurs at day 196 -and is reached around day 245 (2 September) in the simulation.
Although our model does not include drainage, it is reasonable to conclude that runoff extraneous to that ponded in lakes by the observed date of drainage or refreezing has passed into the ice sheet though englacial channels (e.g.crevasses or moulins).Above 1100 m a.s.l., up to the date of observed maximum simulated lake area (day 245), 11.89 km 3 of runoff are produced, according to the MAR model, of which 6 % are simulated to be stored in lakes.A maximum of 12 % of total runoff is estimated to be stored in observed lakes when considering a cumulative lake area, a conical lake and a diameter-depth ratio of 100 : 1 (Box and Ski, 2007).Using this same approximation, the corresponding value for modelled runoff would be 11 %.Fig. 6.Sensitivity analysis with respect to runoff amount for the 1100-1700 m elevation band of lake area.Presented here are lake area profiles simulated using the full runoff amount and scaling factors of 2, 1.5, 0.75, 0.5, and 0.25 applied to the full runoff amount.Also shown is the observed cumulative lake area profile (black dashed) and the maximum possible lake-covered area, calculated by filling all sinks (red dashed).
Because the estimated runoff data have not been independently evaluated, we tested the sensitivity of our model results to the quantity of runoff supplied.Simulations were performed, using a time step equal to 90 s, in which the runoff predicted by MAR was scaled by the following factors: 2.0, 1.5, 1.0, 0.75, 0.5 and 0.25.In locating lakes, little variation is observed, with a range of 62 to 66 % of observed lakes being co-located with observations for scalings of 0.25 to 1 for the MAR-simulated runoff.Similarly, little variation is seen in either the odds ratio or Heidke skill score, with values ranging from 17 to 19 and 0.27 to 0.33 across this range.Maximum simulated lake area is the metric upon which modifying the runoff amount has the most significant impact (Fig. 6).We use ArcGIS to fill the sinks in the DEM to their maximum capacity and find the corresponding lake-covered area to be 6.4 % of the ice sheet above 1100 m a.s.l.Even with double the MAR-simulated runoff supplied to our hydrology model, this value is not reached.Lake onset is initiated sooner with an increase in runoff amount occurring approximately on days 153, 154, 155, 157, 159 and 160 when the runoff is scaled by 2.0, 1.5, 1.0, 0.75, 0.5 and 0.25, respectively.
We ran the model at a resolution of 90 seconds for 2005, 2006 and 2007 (in addition to 2003) to investigate the performance of the model in capturing inter-annual variability and to provide a more robust validation.These years were chosen because they were the years for which observations were available using the same method as those for 2003 (Sundal  2).We compare the cumulative lake area as a percentage of area coverage as simulated (solid) and observed (dashed) in Fig. 7.The best agreement between modelled and observed maximum cumulative lake area occurs in 2003.The model overestimates this value slightly in 2005 and significantly more in the other years considered.Although the maximum lake area is over-predicted, the variability of cumulative lake area is reproduced well for 2005 and 2006; the simulated and observed cumulative area profiles show the same shape.Again we see reasonable correlation in simulated and observed onset day, with the exception of 2007, where correlation is 0.30, which may be considered insignificant.2003 and 2007 were particularly wet years, with a total runoff amount over the whole of this sector of the ice sheet of 12.34 km 3 and 12.99 km 3 , respectively (Table 2).This corresponds to a daily mean runoff depth of 16.52 and 20.77 mm.The 2005 and 2006 melt seasons have over 30 % less runoff than the other two years at 6.93 km 3 and 7.93 km 3 , respectively.

Discussion and conclusions
We have described and evaluated a new method model of modelling supraglacial lake evolution, which uses Darcy's law for flow through a porous medium and Manning's equation for open channel flow to route runoff simulated by the MAR RCM over an InSAR-derived DEM.This model has been used to simulate the onset and growth of SGLs in the Russell glacier catchment in Western Greenland.Our model performs well in predicting SGL location (co-locating 66 % lakes with respect to observations) and maximum cumulative lake area (with an overestimate of ∼ 1.5 %) within 1100 to Table 2. Statistical description of the forecast skill in predicting lake locations in our study area between 1000 and 1600 m elevation in 2003, 2005, 2006 and 2007.Also given are the maximum observed and modelled lake areas as a percentage of total area between 1100 and 1700 m elevation, MAR-simulated runoff and the number of satellite images available for comparison with the model in each of those years.Daily mean runoff relates to the period of days 130 to 250, which we observe to be the period during which supraglacial lakes are seen in this sector of the GrIS.The model predicts a progression of lake initiation inland as the melt season progresses, a pattern that has been reported previously in numerous satellite surveys of SGLs (e.g.Sundal et al., 2009;McMillan et al., 2007;Sneed and Hamilton, 2007).

Multi
We have also used the model of supraglacial lake evolution to estimate a maximum value for the area and volume of SGLs situated at elevations between 1100 and 1700 m a.s.l.within our study area.We agree with the findings of Lüthje et al. (2006) that the local topographic setting limits potential SGL area, and we find this limitation to be 6.4 % of the entire region in our model.We estimate that the maximum volume of water that can possibly be stored in SGLs in this sector of the Greenland ice sheet is 1.49 km 3 or 12 % of all runoff produced in 2003.This suggests that the vast majority of runoff passes through or over the ice sheet without being stored temporarily in lakes.We neither simulate nor observe lake area coverage of this magnitude, even with an assumed doubling of runoff.We note that although simulated lake area coverage stabilises, not all lakes are brimful by the end of the melt season with respect to the sink in which they are located, especially at the higher elevations where melting begins later and ends sooner.If the ablation season were to lengthen, particularly at these higher elevations, it is possible that this topographic limit could be reached.
In a previous study of the nearby Swiss Camp region of the GrIS, McMillan et al. (2007) estimate, using observations of lake area and mean depth, that 17 % of total meltwater volume produced had been stored in lakes by August.By this month, only 6 % of runoff volume produced is simulated to be stored in lakes using our model.We also note a similar discrepancy between our modelled volume estimate and that calculated using a conical area-volume approximation (e.g.Liang et al., 2012;Krawczynski et al., 2009).This, together with the fact that we do not simulate 34 % of the MODISobserved lakes in 2003, suggests that the DEM used in this study underestimates the amplitude of the short period variations in the ice sheet surface, which lead to the formation of lakes.This could arise from the DEM smoothing process or indicate the presence of refrozen lakes misinterpreted as part of the ice sheet surface.This also provides an explanation why a number of lakes we simulate are slightly offset from observed locations.As a result, our 12 % estimate for maximum storage volume should be taken as a lower bound, given a 6.4 % maximum lake-covered area.
We have assessed the ability of our model to reproduce inter-annual variations in supraglacial lake evolution in 2003, 2005, 2006 and 2007.The model shows closest agreement with observations of cumulative lake area described by Sundal et al. (2009) in 2003.Our model predicts the maximum cumulative lake area to be similar in years with comparable runoff; however, observations suggest, particularly between 2003 and 2007, that this is not the case.There are 28 days of observations for 2003 and only 13 per year for the period 2005 to 2007, in addition to which the observations taken in 2007 are clustered around the early and late melt seasons.We suggest that the sparseness of observations from 2005 to 2007 relative to 2003 and the clustering of observations in 2007 mean that observations in these years are not able to capture a comparable degree of variability in daily lake area.We conclude that denser temporal sampling is required to effectively investigate inter-annual variations.
We see good agreement in 2003, 2005, and 2006 between simulated and observed lake locations.These results confirm that the locations of SGLs coincide with intransient depressions in the ice surface topography and, in our experiment, it was not essential to have contemporaneous surface elevation and runoff observations to predict where lakes form.This is in agreement with the findings of Echelmeyer et al. (1991) and Selmes et al. (2011), who used a range of observational techniques to investigate the inter-annual distribution of SGLs on the GrIS.Both studies observed that SGLs appear in the same locations in different years, also indicated in our modelling study.
Our study shows that the location of supraglacial lakes shows little dependence on the amount of runoff supplied to the ice sheet surface when that amount is between half and double that predicted by MAR.For runoff fractions less than 0.5, there is clearly not enough runoff produced to allow all lakes to be filled.These sensitivity tests indicate that lakes would still form in the spatial patterns observed under a wide range of runoff scenarios, allowing for relative uncertainty in the forcing data.The MAR model can, for example, overestimate the quantity of runoff due to the use of a constant bare ice albedo (∼ 0.45) for the entire bare ice area of the GrIS (Fettweis et al., 2011).However, it is unlikely that this would be as much as 100 %.Modelled lake onset days are slightly earlier in the season when runoff is increased, which suggests that in high-runoff years we would expect to see SGLs appearing earlier than in low-runoff years.
It is important to note that our model simulates only lake growth and at present does not incorporate processes leading to rapid lake drainage, which is known to be an aspect of the seasonal cycle of some lakes (e.g.McMillan et al., 2007;Georgiou et al., 2009;Sundal et al., 2009;Selmes et al., 2011).In consequence, while our model provides information about the location, onset, and cumulative area of SGLs, it cannot fully simulate the evolution of lakes that drain.On the other hand, differences between modelled and observed lake volumes can provide useful information as to the quantity of water that has drained from lakes.In future work, processes such as drainage and refreezing will be investigated in the context of this model.Supplementary material related to this article is available online at: http://www.the-cryosphere.net/6/1077/2012/tc-6-1077-2012-supplement.pdf.

Fig. 1 .
Fig. 1.Map of Greenland showing elevation contours (Bamber et al., 2001) at 500 m intervals from 500 m.MAR-simulated total runoff for 2003 is shown in blue (white areas represent zero runoff) and the DEM used in this study is bounded in black.

Fig. 3 .
Fig. 3. (a) Composite plot of all MODIS-observed lakes (red) and all lakes simulated using the transient version of the model (blue) at a 90-s resolution for the 2003 melt season.Coincident lake area is shown in purple.(b) East-west lake surface (blue) and bed (black) profiles for modelled lakes indicated by number in the main plot.The transect taken is shown in black in (a).

Fig. 4 .
Fig. 4. Correlation between simulated and observed onset day for 159 supraglacial lakes.Vertical error bars indicate uncertainty in observed onset day due to non-uniform temporal sampling, size and weight of symbol indicate number of points, n, with indicated value.Line of best fit is shown in red and a 1 : 1 correlation is shown in blue.

Fig. 5 .
Fig. 5. Comparison between modelled (solid), observed (dashed) and cumulative observed (dotted) fractional lake area in 2003 for three altitude bands in the Russell glacier catchment of the Greenland Ice Sheet.

Fig. 7 .
Fig. 7. Multi-year comparison of cumulative lake area (1100 m to 1700 m a.s.l.) as a percentage of total area.Modelled values are given with solid lines and observed values are given dashed.

Table 1 .
Contingency table detailing the number of predicted and observed lake cells in the study area.