Evidence of meltwater retention within the Greenland ice sheet

Greenland ice sheet mass losses have increased in recent decades with more than half of these attributed to surface meltwater runoff. However, the magnitudes of englacial storage, firn retention, internal refreezing and other hydrologic processes that delay or reduce true water export to the global ocean remain less understood, partly due to a scarcity of in situ measurements. Here, ice sheet surface meltwater runoff and proglacial river discharge between 2008 and 2010 near Kangerlussuaq, southwestern Greenland were used to establish suband englacial meltwater storage for a small ice sheet watershed (36–64 km 2). This watershed lacks significant potential meltwater storage in firn, surface lakes on the ice sheet and in the proglacial area, and receives limited proglacial precipitation. Thus, ice sheet surface runoff not accounted for by river discharge can reasonably be attributed to retention in suband englacial storage. Evidence for meltwater storage within the ice sheet includes (1) characteristic dampened daily river discharge amplitudes relative to ice sheet runoff; (2) three cold-season river discharge anomalies at times with limited ice sheet surface melt, demonstrating that meltwater may be retained up to 1–6 months; (3) annual ice sheet watershed runoff is not balanced by river discharge, and while near water budget closure is possible as much as 54 % of melting season ice sheet runoff may not escape to downstream rivers; (4) even the large meltwater retention estimate (54 %) is equivalent to less than 1 % of the ice sheet volume, which suggests that storage in enand subglacial cavities and till is plausible. While this study is the first to provide evidence for meltwater retention and delayed release within the Greenland ice sheet, more information is needed to establish how widespread this is along the Greenland ice sheet perimeter.


Introduction
Greenland ice sheet mass losses derive from increased ice discharge and meltwater runoff since the early 1990s . Increases in meltwater runoff are estimated to be twice as large as cumulative ice discharge anomalies between 2000(Van den Broeke et al., 2009a. This meltwater loss intensification is consistent with observations of rising mean annual near-surface air temperature (+1.8 • C since the Little Ice Age; Box et al., 2009) and expanding melt area on the ice sheet surface (Mote, 2007;Tedesco, 2007), including a record near-total surface melting on 12 July 2012 (Hall et al., 2013;Nghiem et al., 2012). Furthermore, declining ice sheet surface albedo amplifies temperature sensitivity, driving half of the runoff increase between 2000 and 2012 .
Continued mass loss from the Greenland ice sheet has the potential to raise global sea levels by 9 ± 4 cm by 2050  and between 17 and 54 cm by 2100 (Pfeffer et al., 2008). These estimates project that 13-58 % of this increased mass loss occurs as a result of declining surface mass balance (Pfeffer, 2008;Rignot et al., 2011). The wide range in this fraction can partly be attributed to uncertainties associated with quantifying the link between surface mass balance and ice dynamics (Bindschadler et al., 2013), and surface meltwater retention (e.g., Rennermalm et al., 2013).
It has long been known that meltwater transport from the ice sheet surface to its margin occurs through a complex, poorly understood system of supra-, en-and subglacial pathways consisting of crevasses, moulins, fractures, conduits and supraglacial stream channels (Fountain and Walder, 1998), referred to as the cryo-hydrologic system (CHS) (Phillips et al., 2010(Phillips et al., , 2013. Meltwater passage through the CHS and temporary storage in englacial cavities, and proglacial and supraglacial lakes, can retard gravity-driven flow of meltwater from its creation on the ice sheet surface to its appearance at the ice margin (Cuffey and Paterson, 2010;Fountain and Walder, 1998). Meltwater can be retained if it refreezes or is accumulated in firn layers (Boggild, 2007;Boggild et al., 2005;Brown et al., 2011;Fausto et al., 2009;Greuell and Konzelmann, 1994;Pfeffer et al., 2008;Reeh, 1991), which may buffer meltwater losses for several decades , leaving a lesser portion that truly contributes to sea level rise.
Here, ice sheet englacial meltwater retention and duration are estimated for a small watershed (36-64 km 2 ) in the region near Kangerlussuaq, southwestern Greenland. To this end, modeled ice sheet surface meltwater runoff volumes are compared with proglacial river discharge measured 2 km downstream from the ice margin to quantify runoff deficit and surplus. Ice sheet runoff deficit and surplus are directly indicative of internal ice sheet meltwater retention and release, because (1) there are no significant proglacial or supraglacial lake storages, (2) terrestrial fluxes from groundwater and/or snowmelt are limited due to a very small proglacial area (only 6-12 % of the total watershed area), and (3) modeled runoff represents the ice sheet conditions very well as it is estimated from a surface energy balance model forced by nearby automatic weather station (AWS) data Van den Broeke et al., 2011b).

Study region and in situ data
The study catchment is situated along the Greenland ice sheet's southwestern margin ∼ 30 km northeast of Kangerlussuaq, between the Isunnguata Sermia and Russell Glacier (Fig. 1). The study watershed is small (36-64 km 2 ), and mostly ice covered with a small proglacial area (4 km 2 ) consisting of exposed bedrock, tundra, loess deposits, lakes, and river alluvium. The ice sheet surface elevations within this catchment range from 500 to 860 m a.s.l, well below the mean equilibrium line altitude, which over the period 1991-2011 averages 1553 m a.s.l. (Van de Wal et al., 2012), and has no major supraglacial lakes.
Since 2008, hydrometeorological observations have been collected at three sites along the Akuliarusiarsuup Kuua River (Rennermalm et al., 2012) including (1) river discharge, stream water level and stream temperatures (Site AK4, Fig. 1); (2) lake stage fluctuations (Site AK5, Fig. 1); and (3) near-surface air temperatures (Site AK1 and AK2, Fig. 1). Site AK4 is upstream and therefore unaffected by episodic jökulhlaups from a large proglacial lake in this region (e.g., Russell et al., 2011). Meteorological observations on the ice sheet were acquired at three AWSs labeled S5, KAN_L, and S6 in Fig. 1 at 490, 670, and 1020 m a.s.l., respectively (KAN_L from Van As et al., 2011and S5 and S6 from Van den Broeke et al., 2011a), and observations on the tundra were from Kangerlussuaq airport (data access from National Climatic Data Center, 2011). These stations are situated along the K-transect, which is an array of eight surface mass balance observation points operational since 1990, that extends 141 km into the ice sheet ( Van de Wal et al., 2005).

Methods
River discharge was determined by relating 36 discharge measurements with half-hourly stream water level observations (Rennermalm et al., 2012) (rating curve F 1 , Fig. 2b). Data were collected with high-precision Price Type-AA current meters, and Solinst Level-and Barologgers. River discharge uncertainty is 22 % of half-hourly observations, updated from Rennermalm et al. (2012). At Site AK4, river discharge was estimated between 9 June 2008 and 20 April 2011. At Site AK5, which is downstream a neighboring watershed, water levels were measured between 2 June 2007 and 19 August 2008. To include the 2008 melt season onset at the AK4 site, its time series was retroactively estimated with a regression model using AK5 water levels and parameterized with data from an overlapping period between 9 June and 19 August 2008. This relationship was described with a linear relationship (L AK4 = 4.6L AK5 − 15.7 where L AK4 is AK4 water level, and L AK5 is AK5 water level, R 2 = 0.64). While winter discharge uncertainty could not be quantified with conventional methods (Pelletier, 1990;WMO, 2010a, b) due to lack of wintertime sampling, development of a simple uncertainty model provides a lower boundary during freezing temperatures. This model is based on Manning's equation (e.g., Dingman, 2002) simulating a conservative scenario with an ice filled stream channel below the pressure transducer (i.e., no flow beneath the sensor) (Fig. 2a). Manning's equation is simplified similar to Rantz (1983) by assuming that Manning's roughness coefficient, n, and stream surface slope, S, are constant, and by relaxing the 2/3 exponent for hydraulic radius, R h , (i.e., cross sectional area of flow divided with the wetted portion of the flow perimeter): where Q is discharge, a is a constant (equal to S 1/2 /n in Manning's equation), A is stream cross sectional area, and b is the relaxed 2/3 exponent in Manning's equation. Power law parameters a and b were estimated with MATLAB nonlinear least squares fitting routine using all 36 open water observations of discharge and hydraulic radius (F 1 rating curve, R 2 = 0.95). Q during freezing conditions was determined using Eq. (1) with R h and A modified to only include the area above the pressure transducer installation (see sketch in Fig. 2a), and multiplying with an arbitrary factor of 0.8 to account for a situation where S and n would reduce overall winter flow. These modified values of Q were used to construct a power-law discharge rating curve (F 2 rating curve) using MATLAB non-linear least squares fitting (Fig. 2b).
Two time series of discharge were constructed by (1) applying F 1 rating curve, hereafter Q; and (2)  was filtered with a 10 day moving average) when F 2 was applied, hereafter Q LOW . Three possible ice sheet watershed realizations upstream of site AK4 were derived to quantify uncertainty in drainage area delineation. The first watershed was determined with methods derived from terrestrial land surface hydrology approximating hydraulic potential with elevation (e.g., Dingman, 2002) by using ice sheet surface topography to dictate boundaries and water flow directions. The second and third watershed delineations capture the influence of ice overburden pressure on the hydraulic potential by using both surface and basal topography (Cuffey and Paterson, 2010;Lewis and Smith, 2009;Shreve, 1972). Both methods are based on simplifications of ice sheet meltwater flow. The former method assumes that any en-and subglacial flow mirrors supraglacial flow patterns. The latter method is more sophisticated by assuming that water pressure equals ice overburden pressure, but effectively considers the ice sheet to be floating on a film of water (Cuffey and Paterson, 2010). In reality, hydraulic potential near the ice sheet margin is more likely to be spatially heterogeneous, where meltwater flows in a channelized system with lower water pressures (Cowton et al., 2012) as well as in preferential flow paths in pipes and fractures (Cuffey and Paterson, 2010). Thus, none of these three watershed delineations may be representative for the actual ice sheet catchment upstream of AK4, but rather provide a range where it may be found.
Watershed 1 (W1, ice area 60 km 2 ) was delineated from surface topography alone, similar to Mernild et al. (2010a), using the ASTER GDEM2 surface elevation dataset (ASTER GDEM Validation team, 2009). ASTER GDEM2 provides high-resolution gridded elevations (30 m horizontal resolution) with an average vertical precision of 17 m (95 % confidence interval), and agrees well with ICESat (Ice, Cloud, and land Elevation Satellite) surface elevation observations near the ice margin (Carabajal, 2011). Watersheds 2 and 3 (W2 and W3, ice area 31 km 2 and 50 km 2 , respectively) were delineated based on the hydraulic potential surface calculated using ASTER GDEM2, and interpolated point basal topographic data from Multichannel Coherent Radar Depth Sounder (MCoRDS) L2 Ice Thickness (IRMCR2) data (Allen, 2010). Two delineations were made due to uncertainties in basal topography data in this study area, which are sparse, and include low quality and possibly unreliable data. Erroneous basal topography data west of 49.9 • W were identified as points with extremely shallow ice (< 20 m) where thicker ice is expected (judgment from field visits) and ignored in the delineation of W3. However, removal of these points resulted in patchy spatial coverage west of 49.9 • W, which may disproportionately influence any remaining, and unidentified erroneous data (Fig. 3). Thus, W2 was created only using MCoRDS data east of 49.9 • W, of which the high data density ensures better delineation quality (Fig. 3). The W3 ice sheet area is only 57 % of W1 area and should be considered a lowest boundary for the true watershed area. variability is correlated with cloudy conditions and precipitation) and validated against solid precipitation measurements (Van As, 2011; Van As et al., 2012). Any modeled energy surplus results in melt. When the surface is ice covered, or if snow pore space and temperatures in a snow layer do not facilitate refreezing, meltwater and liquid precipitation percolates to a deeper layer and may runoff once it reaches a layer with no pore space. By providing snow depth as a model input, snow refreezing capacity is realistically constrained. Superimposed ice layers (ice forming on top of snow and firn) or internal ice layer causing perched water tables, are marginal processes in this area due to shallow snow depth cover and therefore were not examined in detail. Modeled daily runoff volumes (m 3 d −1 ) were interpolated in 50 m wide elevation bands, and watershed runoff (R W1, R W2, R W3 for W1, W2, and W3, respectively) was calculated as a function of elevation band runoff and its respective area fraction. The 2008 melting season onset predates the Van As et al. (2012) dataset, and was estimated with a linear regression model based on a similar surface energy balance model applied only at K-transect AWS S5 and S6 ( Van den Broeke et al., 2008, 2009b, 2011a. While AWS S5 is on the Russell glacier tongue and known to be influenced by air heated over the surrounding tundra (Van den Broeke et al., 2011a), inclusion of GAP/PROMICE KAN_L and AWS S6 ensures that the watershed energy balance is representative of colder conditions further inland. Given that S5, KAN_L, and S6 are located just outside of W1, W2, and W3, we implicitly assume that elevation is a dominant control on ice sheet runoff.
Characterizations of meteorological conditions on the ice sheet and in the proglacial tundra environment were examined with near-surface air temperatures obtained from the Kangerlussuaq AWS 042310 (data accessed from National Climatic Data Center, 2011), AWS S5 (due to its full time series during the study period), and AK1 and AK2 Barologgers (unshielded and only providing reliable temperature data for winter and early spring). Snow depth was retrieved from AWS S5 and AWS S6.

Results
The lower bounds for winter discharge (Q LOW ) and runoff from the smallest watershed realization (W3) can produce a match between cumulative ice sheet runoff and river discharge over the three-year period (Fig. 4). However, Q LOW and R W3 are conservative estimates and should be considered absolute lowest boundaries, thus water budget closure is unlikely. The other, perhaps more realistic catchments (W1 and W2) produce a surplus of ice sheet runoff relative to river discharge at the end of the three-year period (Fig. 4). Uncertainties of winter discharge measurements exacerbate deficit in river discharge relative to R W1 and R W2 . Thus, two scenarios are plausible: (1) runoff exceeds river discharge over a three-year period suggesting ice sheet meltwater retention (W1 and W2); (2) ice sheet runoff is less than river discharge at the end of the three-year period suggesting an englacial water source. In subsequent paragraphs, analysis of W1 and W3 will be presented because results for W2 and W1 are similar.
For each melt-year within the three-year period, watersheds switch from net release of ice sheet meltwater at the end of melt-year 2008 to net retention at the end of meltyears 2009 and 2010 (Fig. 5). An exception to this occurs in W1 during 2008 when winter discharge uncertainties (Q LOW ) are considered. Winter discharge uncertainty reduces meltwater release, but only in W1, in melt-year 2008, does Q LOW reverse balance between ice sheet meltwater retention (i.e., ice sheet runoff surplus relative to river discharge) and release (i.e., ice sheet runoff deficit). Meltwater retention in W3 is less than in W1, and the W3 water budget is almost in balance in 2009 and 2010 with ice sheet runoff exceeding river discharge resulting in retention of 1-14 % and 15 % of melt-year ice sheet runoff volume, respectively. Given that actual ice sheet runoff probably exceeds R W3 , relatively large net ice sheet meltwater retention in 2009 and 2010 is more likely. In contrast to W3, W1 ice sheet runoff exceeds river discharge resulting in retention of 44-51 % of melt-year ice sheet runoff volume in 2009 and 54 % in 2010. This meltwater generated on the ice sheet surface and not escaping to the river is presumably retained in en-or subglacial storages, given lack of supra-and proglacial lake storage.
Daily summer river discharge and ice sheet runoff during periods with surface melting co-vary strongly with a linear relationship in all years (the correlation coefficient is 0. 46, 0.85, 0.85 in 2008, 2009, and 2010 for Q and R W1 ). Strongest covariability between daily river discharge and ice sheet runoff is found during the warm season (Figs. 6,  7). proportional to watershed size. R W1 tends to be larger and more variable, while R W3 is almost equal to Q for most of the melting season. Watershed runoff and river discharge agreement improve with consideration of winter discharge errors in the early melting seasons of 2008 and 2009. However, considerable discrepancies between river discharge and ice sheet runoff are noted: spring 2008; fall 2008; and spring 2009. Regardless of how watershed runoff is determined, river discharge onset generally precedes ice sheet runoff. In 2008 and 2009, river flow commences within 0-2 days of isolated ice sheet runoff pulses, but is followed by 18-23 days without significant runoff (Table 1). In 2010, river discharge lags ice sheet runoff onset by 4 days. In 2008, river discharge onset is 31 days earlier than ice sheet thaw onset at S5, but occurs within 3 days of proglacial thaw at Kangerlussuaq (Table 1). In 2009 and 2010, opposite patterns are detected when river discharge onset coincides with thaw at S5 (within 1-3 days), but is considerably delayed compared to thaw at Kangerlussuaq (32-35 days).
Three prolonged periods of river flow occur at times with no significant ice sheet surface meltwater production, suggesting a cold-season release of englacial meltwater or proglacial snowmelt (river surplus is shown as negative values in Figs. 6b, c and 7b, c). One period occurs in winter of 2008 between 31 October and 28 November. The two others occur in pre-melting season months (March and April) between 23 March and 16 April 2008 and between 27 April and 15 May 2009. While cold-season river flow event magnitudes are reduced by factoring in wintertime discharge estimation errors, they are still present (Figs. 6c, 7c).
The 2008 winter runoff event is identified as a water pressure anomaly between 30 October and 28 November (Fig. 8). It is accompanied by a rapid temperature increase in the stream channel from −5.5 to −0.5 • C between 26 and 28  October, followed by a period of constant near-zero (average is −0.5 • C) stream channel temperature until 28 November (Fig. 8) indicative of a water phase change from frozen to liquid. This river discharge anomaly is preceded by a small ice sheet runoff pulse between 27 and 30 October 2008,  T AK4 T AK2 T S5 Fig. 8. Ice sheet, river, and proglacial conditions during the October-November 2008 runoff event, including discharge with (Q LOW ) and without winter error compensation (Q), average hourly ice runoff from W1 (R W1 ) and W3 (R W3 ) and sensor water pressure (ρ AK4 ) (a); and temperatures measured in the stream channel (T AK4 ), in the proglacial area by site AK2 (T AK2 ), and on the S5 ice sheet AWS site (T S5 ) (b). A short ice sheet runoff pulse (R W1 and R W3 ) that accompanied above 0 • C ice sheet (T S5 ) and proglacial (T AK2 ) surface air temperatures precedes a longer river discharge event (Q and Q LOW ). Comparison of surface air and stream channel temperatures reveals when the stream channel pressure transducer sensor is exposed to air (15-28 October), and when it is not (28 October-1 December). On 1 November, stream channel pressure increases above background noise and stream channel temperatures rise rapidly to near-zero suggesting that the AK4 stream channel sensor is submerged by liquid water after this date.
coinciding with above-freezing near-surface air temperatures at both AK2 and S5. River discharge peaks on 5 November, seven days after runoff peaks on the ice sheet. Two river discharge events observed in March/April 2008 and April/May in 2009 are also registered as pressure anomalies (Figs. 9, 10). In March/April 2008, discharge is inferred from a sensor at AK5 installed at a lake bottom and cannot be used to determine stream temperatures at AK4 during this anomaly. In contrast, April/May 2009 discharge is determined from a sensor at AK4. In 2009, stream temperature increases by 5 • C on 27 April, in concert with rising water pressure. Although AK4 stream temperatures remain T AK1 T S5 Fig. 9. Same as Fig. 8, but for the 2008 pre-melting season runoff event (due to lack of data T AK5 are shown instead of T AK4 ). A short ice sheet runoff pulse between 22 and 23 March (R W1 and R W3 ) precedes a marked increase in water pressure at AK5 (ρ AK5 ; AK4 site was not operational at this time and Q is extrapolated from site AK5). The ice sheet runoff between 22 and 23 March coincide with a short period with above 0 • C ice sheet (T S5 ) and proglacial temperatures (T AK1 ) at Sites S5 and AK1, respectively. Constant lake temperatures (T AK5 ) during this period indicate unfrozen conditions at the AK5 lake bottom.
below freezing, this may be indicative of a shallow layer of thawing ice by the pressure transducer. Similar to the 2008 October/November discharge event, these two pre-melting season events are preceded by small runoff pulses (22-23 March 2008 and25 April 2009), and short periods (hours) of above 0 • C air temperatures on the ice sheet and in the proglacial environment. Ice sheet melting preceding and during these three coldseason events are considerably smaller than river discharge regardless of watershed delineation, discharge calculation method and uncertainty bounds (Table 2). Even considering the most conservative scenario, using proglacial river discharge with winter discharge errors (Q LOW ) and ice sheet runoff from W1, cold-season events are 297, 7, and 300 % larger than runoff in winter 2008, fall 2008, and winter 2009, respectively. While watershed configuration has little impact on total meltwater losses during these events (runoff from W1 was only 6-13 % larger than W3 during these events), www.the-cryosphere.net/7/1433/2013/ The Cryosphere, 7, 1433-1445, 2013 River discharge Q (10 6 m 3 ) 29 ± 0.2 34 ± 0.2 23 ± 0.2 River discharge Q LOW (10 6 m 3 ) 9.1 ± 0.6 10 ± 0.6 10 ± 0.6 Ice sheet runoff R 1 (10 6 m 3 ) 0.03 ± 0.005 1.4 ± 0.2 0.03 ± 0.005 Ice sheet runoff R 3 (10 6 m 3 ) 0.03 ± 0.004 1.3 ± 0.2 0.03 ± 0.005 Maximum water level (m) 0.96 0.76 0.76 Hypothetical snow depth (m) 9.6 7.6 7.6 (10 % of water density) Hypothetical snow depth (m) 1.9 1.5 1.5 (50 % of water density) T AK4 T AK2 T S5 Fig. 10. Same as Fig. 8, but for the 2009 pre-melting season runoff event. A short ice sheet runoff pulse (R W1 and R W3 ) on 25 April that accompanied above 0 • C ice sheet (T S5 ) and proglacial (T AK2 ) surface air temperatures precedes river runoff onset by 2 days (Q and Q LOW ). Stream channel pressure on 27 April increases simultaneous with a 5 • C rise of stream channel temperatures, which are suggestive of phase change and liquid water flow after this date.
the choice of discharge calculation method has a large impact (discharge without corrections were 123-229 % larger than with corrections). Although precipitation, snowmelt, and/or measurement uncertainties and errors cannot be ruled out entirely, coldseason release of englacially stored meltwater is a likely explanation for stream channel pressure increases recorded during cold-season events. Recorded pressure anomalies correspond to maximum water depths of 0.96, 0.76, and 0.76 m respectively (Table 2), which is well within site AK3's typical flow depths (0-1.4 m). Snow accumulation is a less likely explanation, requiring depths of 7.6-9.6 m and 1.5-1.9 m, using respectively 10 and 50 % of water density (Table 2). Such depths are markedly greater than typical winter snow accumulation at the S5 and S6 AWS sites of 0.16 and 0.33 m, respectively.

Discussion
This study finds that a large fraction of ice sheet surface runoff produced each melt-year may be retained somewhere within the watershed; 52 ± 1.6 % in 2009 and 54 ± 1.4 % in 2010 (Fig. 5a). In 2008 and 2009, 4-27 % (12-37 % with Q, 10-20 % with Q LOW ) of ice sheet meltwater is released with a 1-6 month delay in the cold season, but cannot close meltwater input and output gaps, except if the ice sheet catchment areas were very small (W3). Uncertainties associated with runoff and discharge estimation prohibit discerning whether net meltwater retention or release took place in 2008 (Fig. 5a). Large net annual meltwater release in 2008 is only possible from an upstream ice sheet watershed that is very small (i.e. W3, Fig. 5b). Few basal topographic data and low quality data near the ice margin suggest that the true watershed is larger than W3. Even with availability of high quality basal topography data, calculations of hydraulic potential ignore preferential flow paths in en-and subglacial conduits thereby only providing an approximation of the true watershed. Assuming that the true catchment is larger than W3 and smaller than W1, and that effective ice sheet watershed area remains the same from one year to another, meltwater is likely to be retained in the ice sheet during the three-year study period.
Lack of substantial surface meltwater lakes on the ice sheet and in the proglacial study area indicate that surface meltwater not exiting in the AK River must have been stored in subglacial till or ice sheet macropores (Fountain and Walder, 1998). Glacier macropores, i.e., the volume fraction of crevasses and fractures, is a potential storage mechanism for some retained meltwater. Typical macroporosity values in glaciers are ∼ 1 % (Anderson, 2004;Fountain and Walder, 1998;Harper and Humphrey, 1995), and must be augmented by storage in subglacial cavities and till, since the volume fraction of ice sheet runoff surplus in 2009 (0.3 %) and 2010 (0.5 %) would saturate macropores in 2-3 yr. Some retained meltwater probably refreezes during winter, but cold-season discharge events show a substantial amount of unfrozen water throughout winter. Some ice sheet meltwater may surpass the river channel altogether and be transported as groundwater flow and in the river hyporheic zone, but the presence of permafrost probably limits the magnitude of this loss. Evidence for delayed ice sheet surface meltwater release is found on daily, seasonal and multiyear timescales. First, ice sheet surface runoff and downstream river discharge have high co-variability in summer months. Amplitude of river discharge is dampened relative to ice sheet runoff, which is typical of meltwater transformation in supra-, en-, and subglacial environments during its passage to rivers at the ice terminus (Fountain and Walder, 1998). Second, three coldseason meltwater releases during times with insignificant ice sheet surface melting suggest that the CHS system can retain meltwater for months after the end of the melting season in September. Finally, ice sheet surface runoff and river discharge are not in balance, resulting in net meltwater storage two years in a row, suggesting meltwater retention on timescales beyond one year. This is verified by modeling of ice sheet hydrology in a western Greenland watershed, where unfrozen ice sheet meltwater was retained internally over multiple years (Colgan et al., 2011). Delays between ice sheet runoff and downstream river discharge on monthly and multiyear timescales are not unique to the Greenland ice sheet, and have been identified for other Arctic glaciers (Hagen et al., 2003;Hodson, 2005;Jansson et al., 2003;Stenborg, 1965;Wadham, 2000).
Three cold-season meltwater releases were only accompanied by modest ice sheet surface runoff, but two pieces of evidence show that river channel flow indeed did occur between 31 October and 28 November 2008, 23 March and 16 April, and 27 April and 15 May. First, stream channel pressure anomalies during cold-season events cannot be explained by the weight of dense snowpacks; instead, derived water levels and runoff are well within AK4's natural range. While a very conservative method to adjust winter discharge errors reduce cold-season river discharge, it does not eliminate flow at this time of year. Second, stream channel temperatures suggest liquid water, and/or phase change. During October 2008, stream channel temperatures co-vary with nearsurface air temperatures until 31 October because there is no water evident at AK4 (Fig. 8). However, after 1 November, stream channel temperatures rise and air temperatures remain below freezing, coinciding with stream channel pressure increase indicating flowing unfrozen water. During the 2009 pre-melt season, a sudden 5 • C stream channel temperature increase coincided with raised stream channel pressure to indicate melting and flowing water (Fig. 10). In the premelting season of 2008, meltwater release precedes ice sheet surface melting and above-freezing near-surface air temperatures (Fig. 9). Thus, despite inherent uncertainties of wintertime river low-flow observations (Pelletier, 1990) and assuming extremely reduced flow, existence of cold-season meltwater discharge is evident. Measurement uncertainties and errors are unlikely to completely account for this given that the sensor operated as expected before and after the three cold-season discharge episodes, and coherence with brief preceding ice sheet runoff events providing a potential triggering mechanism. Presence and release of unfrozen meltwater in the Greenland ice sheet at subfreezing temperatures suggests that parts of the CHS are intact in the cold season. Indeed, this possibility is confirmed by modeling studies (Colgan et al., 2011;Phillips et al., 2010Phillips et al., , 2013, and observations with ground penetrating radar (Catania and Neumann, 2010).
Above-freezing air temperatures may provide a triggering mechanism for release of stored meltwater. During these times, a short thaw period (1-4 days) started at the S5 AWS site 30-75 h before they were detected in the stream. Although bursts of ice sheet runoff accompany all events, their durations are too short (1-3 days versus 18-30 days) and volumes are too small (never exceeding 14 % of river discharge volume even when the most conservative estimation techniques were used) to explain all river runoff observed during these periods. A lack of significant ice sheet runoff during these months indicates that cold-season river discharge was produced during warm summer months. Additional drivers for winter release may stem from buildup of subglacial pressure as englacial meltwater drains to subglacial cavities and drainage exits close (Irvine-Fynn et al., 2011).
It is known that meltwater may be retained on the ice sheet surface in lakes (Selmes et al., 2011;Sundal et al., 2009) and percolate into firn layers . However, evidence presented in this study suggests that meltwater transported through englacial conduits (Catania and Neumann, 2010;Catania et al., 2008) to storage cavities in englacial and subglacial environments may remain there until short thaw events presumably trigger its release. Such subglacial storages, and delayed release, have been observed for high Arctic glaciers in Svalbard (Hodson, 2005;Wadham et al., 2001), but not previously identified for the Greenland ice sheet. In contrast, good agreement between ice sheet runoff and river discharge was found for the large Watson River watershed (∼ 12 500 km 2 ) in southwestern Greenland in 2009 and 2010 . While ice sheet runoff in the Watson River watershed was larger than observed discharge (8 % in 2009, and 19 % in 2010), the AK4 study site examined here exhibited much larger ice sheet runoff (up to 118 % more than river discharge in [2009][2010]. Because the AK4 site is a subwatershed within the Watson River, this discrepancy could suggest that the AK4 site retention cannot be upscaled to the entire Watson River watershed. However, this cannot be established without Watson River winter discharge observations, which are currently unavailable. Regardless, an increase in ice sheet runoff surplus relative to river discharge from 2009 to 2010 in the Watson River indicate a non-systematic underlying driver, for example, en-and subglacial meltwater storage. Cold-season release may in fact be a consequence of meltwater storage within the Greenland ice sheet according to the following hypothesis: meltwater retention builds up subglacial pressures such that CHS, in the cold-season near the ice margin, remains largely intact and can rapidly reactivate in response to short-lived melt events. This hypothesis is not unique to this study. Six years of energy and mass balance studies from a Svalbard glacier reveal that years with lower than expected ice sheet meltwater export (due to internal storage) preceded years when less than usual energy was required for activating the subglacial system, perhaps due to increased subglacial pressure from larger than usual internal storage (Hodson, 2005). It is unclear how meltwater retention influences CHS evolution. Modeling and observational studies suggest that CHS seasonally evolves from an un-channelized state dominated by distributed cavities to an efficient channelized system with large conduits controlled by meltwater fluxes to the subglacial hydrologic drainage system (Schoof, 2010;Sundal et al., 2011). This influences ice sheet velocities, so that in years with strong surface melting, and thus more efficient subglacial drainage, ice sheet velocities decelerate early in the melt season (Sundal et al., 2011). To understand the importance of Greenland ice sheet meltwater retention on these processes, more observational and modeling studies must establish how common this is in other parts of Greenland, and how it is related to subglacial pressure, cold-season releases, and CHS development.

Conclusions
Greater cumulative ice sheet surface runoff relative to downstream river discharge suggests meltwater retention in enand subglacial storages up to as much as 54 % of annual ice sheet runoff. However, estimation errors, particularly associated with watershed delineations, may reduce this amount and explain discrepancies between runoff and downstream river discharge resulting in very small retention. While the upper estimate of meltwater retention is a large fraction of annual ice sheet runoff (e.g., 54 %), it is less than 1 % of total watershed ice volume. This is similar to observed values of glacier macroporosity, thus macropores, and sub-and englacial storage provide reasonable storage mechanisms for excess ice sheet runoff. Furthermore, alternative storages are unlikely due to lack of ice sheet firn, surface lakes, and limited proglacial area.
Three observed river runoff events outside of the regular summer melting period, and preceded by very small ice sheet runoff fluxes, are examples of 1-6 month delayed release of en-and subglacially stored water. These cold-season events cannot be explained by proglacial lake drainage (there are no significant proglacial lakes) or snow melt. Instead, stream channel temperature recordings provide strong evidence for water flow during these events. Regardless, these delayed meltwater releases do not close the watershed water budget, indicating multiyear retention of ice sheet meltwater, and evidence that parts of the Greenland ice sheet CHS may remain active year round.