Hypsometric amplification and routing moderation of Greenland ice sheet meltwater release

Abstract. Concurrent ice sheet surface runoff and proglacial discharge monitoring are essential for understanding Greenland ice sheet meltwater release. We use an updated, well-constrained river discharge time series from the Watson River in southwest Greenland, with an accurate, observation-based ice sheet surface mass balance model of the  ∼  12 000 km2 ice sheet area feeding the river. For the 2006–2015 decade, we find a large range of a factor of 3 in interannual variability in discharge. The amount of discharge is amplified  ∼  56 % by the ice sheet's hypsometry, i.e., area increase with elevation. A good match between river discharge and ice sheet surface meltwater production is found after introducing elevation-dependent transit delays that moderate diurnal variability in meltwater release by a factor of 10–20. The routing lag time increases with ice sheet elevation and attains values in excess of 1 week for the upper reaches of the runoff area at  ∼  1800 m above sea level. These multi-day routing delays ensure that the highest proglacial discharge levels and thus overbank flooding events are more likely to occur after multi-day melt episodes. Finally, for the Watson River ice sheet catchment, we find no evidence of meltwater storage in or release from the en- and subglacial environments in quantities exceeding our methodological uncertainty, based on the good match between ice sheet runoff and proglacial discharge.


Introduction
The majority of recent Greenland ice sheet mass loss can be attributed to increases in surface melting (Enderlin et al., 2014). Several methods exist to determine ice sheet mass balance, including remotely sensed changes in gravity (e.g., Velicogna, 2009). Yet gravimetric estimates alone cannot distinguish between changes in surface mass balance (SMB) and dynamic mass loss through solid ice discharge. Conventional tools for determining Greenland-wide SMB are atmospheric and land surface models forced by atmospheric reanalysis data (Fettweis, 2007;Ettema et al., 2009;Langen et al., 2015). Regionally, SMB models can be driven by in situ weather station observations . These models resolve the components of the surface energy and mass budgets, and therefore allow evaluation of physical processes impacting SMB, such as snow accumulation, surface melting, refreezing in snow and firn, and sublimation.
Another method for studying regional SMB is by monitoring proglacial river discharge comprised of melt-and rainwater exiting the ice sheet without getting retained in supra-, en-, and subglacial environments. Along the land-terminating sectors of the Greenland ice sheet, hundreds of proglacial rivers transport sizeable meltwater fluxes to surrounding seas. Studies of proglacial river discharge were initiated during the 1970s and 1980s in response to the growing interest in Greenland hydropower, but in recent years attention has turned to understanding ice sheet surface mass balance and hydrology (Mernild and Hasholt, 2009;Bartholomew et al., 2011;Rennermalm et al., 2012;Banwell et al., 2013;Van As et al., 2014;Overeem et al., 2015;Smith et al., 2015). Yet still relatively little is known about Greenland ice sheet freshwater discharge in terms of delays imposed by limitations in hydrological efficiency or because of supra-, en-, and subglacial storage.
The complexity of freshwater discharge delay partly stems from the various supra-, en-, and subglacial transit processes that are involved, most notably subglacial drainage channel opening and closure in response to meltwater supply (Fountain and Walder, 1998). Studies quantifying the freshwater transit time between ice sheet surface generation and release at the margin are scarce. In terms of surface routing before entering a moulin, Zuo and Oerlemans (1996) hypothesized that drainage of surface water takes between 1 and 26 days depending on surface slope. For several moulin sites, Chandler et al. (2013) determined en-and subglacial flow velocities to be in the 0.2-1.6 km h −1 range, resulting in routing delays from hours to days, depending on the distance traveled and the efficiency of the subglacial drainage system (Meierbachtol et al., 2013). Subglacial lakes are known to exist and drain in Greenland (e.g., Palmer et al., 2015), but we have no precise estimate of the contribution of such retention to regional ice sheet discharge on intra-and interannual timescales. Efforts have been made using the input-output method, in which the difference between ice sheet surface runoff and proglacial discharge provide an estimate of retention. For instance, Rennermalm et al. (2013) calculate retention of up to half of the water running off the surface of the ice sheet.
In this study, we apply the input-output method to determine the routing delays as a function of meltwater origin elevation to quantify ice sheet hydraulic transmission efficiency. Our study area is the relatively arid sector of the Greenland ice sheet east of the Kangerlussuaq settlement, located in southwest Greenland (Fig. 1). We make use of time series of supraglacial meltwater production as calculated by an in situ, observation-driven SMB model and observed proglacial river discharge. Both time series stem from earlier research, but have been revised to increase absolute accuracy. Furthermore, by studying the large "Kangerlussuaq catchment" of the ice sheet, the relative impact of errors in catchment delineation is reduced. In interpreting the discharge time series, special attention is given to how ice sheet discharge responds to climate variability. Previous studies have identified the role of the area-elevation distribution, also known as hypsometry McGrath et al., 2013;Mikkelsen et al., 2016), but have not quantified its effect as we do here. Hypsometry is expected to be an amplifier of meltwater runoff from the Greenland ice sheet in a warming climate.

River discharge
The Watson River drains an ice sheet area of ∼ 12 000 km 2 (Lindbäck et al., 2015) through confluent branches into the "Kangerlussuaq fjord" (Fig. 1; Nielsen et al., 2010). River stage has been monitored near the bridge in Kangerlussuaq, southwest Greenland, since 2006 (Fig. 2;Hasholt et al., 2013). The location of pressure transducers 140-150 m upstream of the bridge in between the two main river channels was chosen to be as close to the river bottom as possible to capture low water levels and sufficiently far upstream to avoid influence of the drop in water level across the stable control section of the river. Water pressure measurements are adjusted for barometric pressure before converting into hourly averages of water stage. River discharge measurements taken over a range of water stages are required to construct a rating curve for converting stage into discharge (e.g., Rennermalm et al., 2012). For this purpose, we measure discharge in three different ways : 1. At very low water levels in spring we wade across the river measuring water velocity with an OTT C2 propeller current meter. The depth-and width-integrated water flux has an uncertainty < 10 % combining sensor error, and bank-to-bank and depth integration of the point measurement, when the full width of the river is surveyed.
2. From the bridge, we release a tethered float to estimate the depth-and width-averaged flow velocity based on travel time and distance. From 2007 to 2015 ∼ 200 float measurements were obtained from which we select the 13 measurements for which we can constrain the crosssectional area from depth probing. Only a few measurements pass this criterion because during intermediate to high flow, depth probing by tethered weight or iron rod is impossible, whereas we know there to be depth variations in channel 1 (Fig. 2) in excess of 2 m due to erosion and deposition of bed load (sediment and gravel; Hasholt et al., 2013). One high-discharge 2011 measurement is available in which a boat was used as a float in a wider, echo-sounded section of the river. Hasholt et al. (2013) estimate the float method to be accurate within 15 %, representing the combined uncertainties in cross-sectional area, surface velocity determination, and calculation of depth-average velocity, although below we argue this to be an underestimate.
3. In 2012, 74 river crossings were done to perform acoustic Doppler current profiler (ADCP) measurements by means of a 600 kHz WorkHorse Monitor ADCP from Teledyne RDI, mounted downward-looking over the side of a boat. Concurrently geographic position of the boat is determined using a hand-held GPS receiver. For  each river crossing the discharge is calculated by integrating the true flow velocities from bank to bank, taking into account the direction and the length of the survey path. Where the flow velocities could not be measured, primarily at the surface, bottom, and banks, they are assumed to be equal to the mean flow velocity at the location at hand. The combined uncertainty of each depth-and width-integrated river discharge value, including sensor error in sediment-loaded water, is estimated to be at most 5 % + 100 m 3 s −1 .
The three methods provide 91 near-instantaneous river discharge values to construct a stage-discharge relation. In Fig. 3 we plot discharge versus stage measured at the pressure transducer site. Least-squares fitting of a power function reveals the stage-discharge relation: with river discharge D in m 3 s −1 and stage H in m units. The exponent value has theoretical foundation in falling within the 2-3 range that is common in hydraulic and fluvial morphology (Hershey, 1999). We find a fit correlation of r = 0.994 and a root mean square difference of 72 m 3 s −1 . An uncertainty of 8 % encompasses all ADCP measurement uncertainties, but we use a conservative uncertainty value of 15 % for converting stage into discharge. The best fit suggests that the bottom of the stable bedrock control cross section is 1.97 m below the lowest pressure transducer. Here no sediment is observed to accumulate towards the end of the melt season. In autumn, ice does accumulate in the control cross . Measured Watson River discharge versus stage. Symbols denote discharge measurements by three methods. The black line represents the best power-function fit, and the grey area a 15 % uncertainty range. Dotted lines illustrate earlier versions (v1 and v2) of the stage-discharge relation.
- Figure 4. Hourly values of river discharge versus 10-day smoothed air temperature in Kangerlussuaq for July-September data 2008-2015 (dark grey). The solid line denotes a power-function fit to the data, with the light grey area the 70 % uncertainty range. Icedammed lake drainages (jökullaups) stand out as spikes.
section that allows for a stable stage-discharge relation but melts shortly after the river starts to flow in spring. ADCP measurements were largely unavailable to previous studies of Watson River discharge (Mernild and Hasholt, 2009;Hasholt et al., 2013). The availability of ADCP data considerably changes the stage-discharge relation (Fig. 3), effectively doubling discharge from Hasholt et al. (2013) and quadrupling the values by Mernild and Hasholt (2009). The two earlier relations differ from each other because of a revision in the cross-sectional area and larger availability of floatderived measurements at high stage - Hasholt et al. (2013) used ∼ 140 float measurements to construct a rating curve. Yet because of the persisting problem of measuring channel depth during intermediate and high discharge, we have to reject the (majority of the) float-derived discharge values for our study because the cross-sectional area could not be determined simultaneously. We speculate that our new rating curve exceeds the previous one due to remaining uncertainties in cross-sectional area in Hasholt et al. (2013) and/or uncertainties in deriving the channel average velocity from surface float measurements in the rapid, supercritical flow through the irregular and seasonally heavily sedimented channel 1 (Fig. 2). Neither error source affects ADCP measurements taken elsewhere in the river.
An independent method to validate river discharge values is presented by the occasional freshwater outbursts (or jökullaups) of an ice-dammed lake at the ice sheet margin within the Kangerlussuaq catchment (Russell et al., 2011), causing distinct peaks in river discharge (Fig. 4;Mernild and Hasholt, 2009;Mikkelsen et al., 2013). Lake volume changes during these outbursts can be compared with the time-integrated spike in river discharge at the Kangerlussuaq bridge, after adjustment for background flow. Russell et al. (2011) determined ice-dammed lake drainage volumes to be 39.1 ± 0.8 × 10 6 in 2007 and 12.9 ± 0.3 × 10 6 m 3 in 2008. Discharge values calculated from the previous stagedischarge relations underestimated these jökullaup volumes. Our updated stage-discharge relation provides jökullaup volume estimates close to those by Russell et al. (2011), totalling 43.1 ± 8.6 × 10 6 in 2007 and 9.4 ± 1.9 × 10 6 m 3 in 2008. The 2008 event is underestimated even using the updated stage-discharge relation, suggesting that a further increase in the river discharge calculation at low river stage may be appropriate.

Gap filling: temperature-based discharge
Due to the risk of frost damage when encased in solid ice, the pressure transducers recording stage cannot remain installed through winter. Consequently, early-and late-season periods exist during which water stage is not recorded. This includes instances during which water stage falls below the level of the pressure transducers. A significant data gap also exists in the 2006 data record when water stage exceeded the pressure transducer's measurement range.
To estimate river discharge during data gaps, we use hourly air temperature data collected in Kangerlussuaq (Cappelen, 2016) as an ice sheet melt proxy. Whereas local melting can be approximated by a linear temperature-index model, for the whole catchment draining in the Watson River we find that a power law approximates the relation between river discharge and air temperature (Fig. 4): where D T is temperature-derived discharge (in m 3 s −1 ) and T 0 is Kangerlussuaq air temperature (in • C). To roughly account for transit time between the ice sheet surface and the river monitoring site, we apply a 10-day smoothing and 5day delay to the temperature data. We distinguish between the first half of the year (up to and including June) and the second half because of differences in winter-accumulated snow on tundra and ice sheet, impacting meltwater retention and retardation. Figure 4 illustrates that the temperature response factor F T equals 0.31 during the peak and late melt season (July and after, when the ice sheet ablation area has little to no snow cover), yielding a correlation of r = 0.74. F T is smaller (0.17) during the first half of the year (r = 0.65). The correlations are good considering that temperature cannot serve as a proxy for rain, rapid lake drainages, and other factors governing river discharge variability. We set a large, conservatively chosen ±70 % uncertainty range (Fig. 4) to encompass most discharge measurements except during icedammed lake drainages. We return to these equations later when we use them for quantifying ice sheet hypsometric amplification of meltwater release.
In gap filling, D T makes up 69 % of the 2006 river discharge total but adds up to just 1.2 % of the discharge spanning 2007-2015. Hasholt et al. (2013) estimated the river discharge to be 4-11 % of the annual values during these data gaps in 2007-2010. We find a smaller contribution during data gaps largely due to the revised stage-discharge relation.

Surface mass balance modeling
Supraglacial runoff is calculated by an improved version of the SMB model by  and has been used successfully for different glacier settings around the globe (e.g., Van As et al., 2005;Van den Broeke et al., 2008). The model has the advantage of being forced with local observations, as opposed to, e.g., regional climate models that are constrained at remote boundaries. The model inter-/extrapolates meteorological and radiative measurements from three automatic weather stations placed at different ice sheet elevations into 100 m elevation bins and calculates the surface energy and mass components in each bin, ranging from the margin to 2000 m above sea level (a.s.l.), above the surface runoff limit. For every time step, the model iteratively solves the surface energy balance for the surface temperature. If surface temperature is limited by the melting point, the surplus energy is used for melting of snow or ice.
In calculating turbulent heat fluxes, we set aerodynamic surface roughness for momentum to commonly accepted and observed values of 0.02 and 1 mm for snow and ice, respectively (Van As et al., 2005;Smeets and Van den Broeke, 2008). We adopt a snow density value of 500 kg m −3 after Van den Broeke et al. (2008). Snow and firn densify and gain heat in the model through the refreezing of meltwater that percolates from the surface, provided cold content is available (Charalampidis et al., 2015) and that no "impenetrable" ice layers exceeding a 1 m thickness are encountered. The model is initialized in April 2009 with linearly thickening firn with elevation (0.14 m m −1 ) on top of solid ice above the long-term equilibrium line altitude at ∼ 1550 m a.s.l. There is a lack of precipitation measurements over the Kangerlussuaq catchment of the ice sheet. Therefore our model includes a precipitation parameterization in which a constant precipitation rate is assumed for snowfall (air temperature below freezing) and rainfall (above freezing) when downward long-wave radiation exceeds the blackbody emissions calculated from air temperature. This rough estimate is partly justified by the small impact precipitation has on the outcome of this study because of the arid climate that governs the lower elevations of the Kangerlussuaq catchment (Van den Johansson et al., 2015). The precipitation rate is tuned to provide optimal results in terms of winter accumulation. Importantly, the model does not use surface albedo measured at the weather stations, as it is spatially heterogeneous while highly influential in model calculations. Instead we use Moderate Resolution Imaging Spectroradiometer (MODIS) Terra MOD10A1 albedo within the Kangerlussuaq catchment, averaged over the 100 m elevation bins utilized by the model, after removing spikes. One of the two major improvements made to our model is that we calibrated MODIS albedo to 5 years of daily albedo measured at the weather stations of the Greenland Climate Network (GC-Net) and the Programme for Monitoring of the Greenland Ice Sheet (PROMICE). This is needed because we find from linear regression that MODIS on average gives albedo values lower than those derived from in situ observations for solar zenith angles below ∼ 74 • (mean bias of 0.043): where α is albedo and θ is the solar zenith angle. The second large improvement to the model is an increase in temporal resolution. While for Van As et al. (2012) the daily MODIS resolution was an argument for running the SMB model in daily time steps, we recognize the need to resolve the daily cycle in ice sheet runoff. Therefore the current model version runs at an hourly time interval with a fixed daily albedo. The two above-mentioned important changes to the Van As et al. (2012) model mostly increase the melt and runoff in our calculations. The SMB model is not in any way tuned to match river discharge.
To test the model's performance, we compared its calculations of ablation and accumulation with independent in situ measurements from the three weather stations. Over the course of seven melt seasons, the accumulated SMB at any time is modeled within 1.5 m ice equivalent of the measured values in spite of 20-40 m elevation differences (Fig. 5). At the end of the model run, the accumulated model error at the lower and middle weather station site is negligible; at the upper site the model underestimates SMB by 1 m of snow, i.e., half a meter in water equivalent. The model annually overestimates winter snow accumulation by 0-0.5 m snow/ice equivalent at low elevation; Van den Broeke et al. (2008) suggest that not all accumulation gets recorded by low-elevation weather stations because snow mostly collects in the depressions between the ∼ 5-10 m diameter ice hummocks. At high elevation, winter accumulation and summer ablation (melt) are overestimated by up to 1 m ice equivalent in some years, which partially cancel each other out in terms of SMB (Fig. 5). Note that observations as well as model results place the elevation at which the climatological mass budget is in balance (SMB + refreezing = 0) around 1800 m a.s.l. for 2009-2015, which is ∼ 250 m higher than reported for the 1990-2011 period, and similar to elevations in the top-ranking years 1995, 1999, 2003, 2007(Van de Wal et al., 2012. Lindbäck et al. (2015) delineated the "Kangerlussuaq catchment" by means of hydraulic potential analysis, which is a steady-state proxy for routing of subglacial water (Shreve, 1972). The catchment stretches from the ice sheet margin to the ice divide with a total area of ∼ 12 000 km 2 , or 0.7 % of the total ice sheet area. This delineation method is superior to methods entirely dependent on surface slope (e.g., Van As et al., 2012) as moulins and crevasses are abundantly present far from the ice sheet margin in the Kangerlussuaq region, yielding subglacial meltwater routing over large distances (Yang et al., 2016). There is evidence of meltwater reaching the glacier bed at over 130 km from the ice sheet margin where the surface is 1840 m a.s.l. ). Yet ice sheet catchment delineation remains a relatively uncertain undertaking, e.g., due to uncertainties in bedrock topographical maps (Carroll et al., 2016) that were minimized by Lindbäck et al. (2015) by constructing a detailed bedrock map from ice-penetrating radar measurements. To our advantage, the Kangerlussuaq catchment is relatively wide at 30+ km near the ice sheet margin (Fig. 1), making the total runoff from the area less sensitive to boundary shifts. Also, the large majority of meltwater is generated at low elevation, where catchment delineation is least uncertain due to its proximity to the watersheds at the ice sheet margin.

Kangerlussuaq catchment delineation
We subdivide the catchment area into 100 m elevation bins to match the output of the SMB model, allowing multiplication with calculated ice sheet runoff (in m) to derive volumetric units. We do not take the ice-free catchment of Watson River into account (Weidick and Olesen, 1978). Tundra makes up a small fraction (5 %) of the total catchment area (Mikkelsen et al., 2016), contributing on average less than 0.1 km 3 yr −1 (in the order of 1 %) to discharge from precipitation, not counting losses from evaporation.

Meltwater runoff delays
The time series of proglacial river discharge and ice sheet surface runoff enable the calculation of routing delays per elevation bin, introduced by meltwater transiting the supra-, en-, sub-, and proglacial environments. We achieve this by finding the highest correlation (r) between (observed, not temperature-derived) river discharge and catchmenttotal runoff at time t (R(t) = N n=1 R n (t)) after introducing different time delays for every elevation bin: R delay (t) = N n=1 R n (t − d n ), shifting the R n time series in elevation bin n by d n (full) hours before summing. In calculating correlations, we loop through all possible d n values. We limit the search by setting d n ≤ d n+1 , i.e., it cannot take shorter for meltwater to transit from higher elevation, and d n+1 −d n ≤ 36, i.e., the added time delay is at most 36 h compared to that of the neighboring lower elevation bin. For 18 elevation bins (covering 50-1850 m a.s.l., up to the observed maximum elevation of the surface runoff line) and 37 delay hours to test per bin, this requires 37 18 calculations of correlation, indicating the need for simplification to reduce computing time. Therefore we test routing delays only every second or third elevation bin, interpolating delays linearly for intermediate bins. We apply 1-, 2-, or 4-hourly time steps depending on the elevation of the bin.
Optimal discharge delays are defined as the mean of those for which correlation r falls within 0.01 of maximum correlation. If more than 100 solutions fulfil this criterion, we calculate the average of the cases with the 100 highest correlations. The standard deviation serves as a measure of spread in the results. We determine the optimal routing delays for every June, July, August, and September of every year with overlapping river discharge and ice sheet runoff time series (2009)(2010)(2011)(2012)(2013)(2014)(2015). We also determine the optimal delays for the entire 7-year time series. This correlation-based method is insensitive to errors in the magnitude of ice sheet runoff and river discharge (and thus in the stage-discharge relation), as correlation is a measure of covariability.

Hypsometric amplification
In deriving an equation to quantify hypsometric amplification of ice sheet discharge, we assume that river discharge can approximate the product of meltwater M (in vertical units) and surface area A, integrated over the melt zone of the ice sheet ranging from the margin at Z 0 up to the upper melt limit at Z. This is a fair assumption since meltwater retention in firn is of second-order magnitude for catchmenttotal runoff in this catchment even during extreme melting (Machguth et al., 2016). Further assuming that (1) M is a linear function of (positive) temperature T and multiplier F M , (2) T is a function of elevation z using linear lapse rate λ and T (Z 0 ) = T 0 , and (3) A is a power function of surface elevation with multiplier F A and exponent p, we find Thus, river discharge can be approximated from air temperature T 0 at Z 0 provided the ice sheet hypsometric amplifier p is known along with the melt (F M ) and area (F A ) factors. For the Kangerlussuaq region, air temperature is measured at the approximate elevation of the lowest glacier margin (∼ 50 m a.s.l.; Cappelen, 2016). We derive p and F A from the area-elevation distribution of the Kangerlussuaq catchment as delineated by Lindbäck et al. (2015). We also determine p for the entire Greenland ice sheet and its northern (> 77.13  3 Results

Surface runoff and river discharge variability
Ice sheet surface melt in the Kangerlussuaq ice sheet catchment, and thereby Watson River discharge, is confined to the May to September period, with minor episodic exceptions ( Figs. 6 and 7). River discharge in our observational period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) peaks between 11 July and 2 August each year (Table 1), with a median of 25 July. Peak discharge ranges by a factor 2.7 from a low peak of 1.2 × 10 3 in 2015 to a high peak of 3.2 × 10 3 m 3 s −1 in 2012. Annual totals range by a factor of 3 from 3.8 in 2015 to 11.2 km 3 in 2010 (Table 1, Fig. 6). These values match those of ice sheet runoff within uncertainty; both time series show an equally large interannual variability (Fig. 6b). Also intra-annual variability in river discharge is large, as values can increase by over a factor of 2 over the course of a few days (Fig. 7). Discharge peaks roughly double their normal values for the time of year are explained by intense ice sheet melt episodes, such as those in late summer 2011 (Doyle et al., 2015) and mid-July 2012 Mikkelsen et al., 2016). Diurnal variability in river discharge is on the order of 200 m 3 s −1 , or 10-20 % of the total signal. Meanwhile, Kangerlussuaq catchment ice sheet melt regularly exceeds 3000 m 3 s −1 during mid-day while often halting at night, thus displaying a 10-20 times larger diurnal variability than river discharge. Many other aspects of the runoff and discharge time series can be identified in Fig. 7, such as the influence of rain or rapid drainage by supraglacial and ice-dammed lakes. We will return to these topics in the discussion section. For now, we focus on two science questions that arise from the results in terms of spatial and temporal variability: 1. Can we explain the large variability in ice sheet meltwater release by quantifying hypsometric amplification?
2. Can we explain the moderation of ice sheet meltwater release by quantifying routing delays in the supra-, en-, sub-, and proglacial environments?

Hypsometric amplification
To investigate the hypsometric amplification of meltwater release through increases in ice sheet area with elevation, we return to Eq. (4), which allows us to quantify the temperature response of ice sheet discharge provided we determine the value for hypsometric amplifier p. For a linear surface slope p = 0, while p < 0 for a convex hypsometry typical for mountain glaciers. The higher the p value, the more sensitive the ice mass is to atmospheric temperature increase. Hypsometric amplification (p > 0) is known to occur for the Greenland ice sheet McGrath et al., 2013;Mikkelsen et al., 2016). From the geometry of the Kangerlussuaq catchment (Lindbäck et al., 2015) we deduce that A = F A · (z − Z 0 ) p ≈ 160 · (z − 50) 1.4 is a good (r = 0.993) approximation for the area below 1350 m a.s.l., which generated 85 % of all surface runoff for 2009-2015 according to our model calculations. Correlation reduces to r = 0.80 when including the nearly twice as narrow region between 1350 m and the maximum observed runoff elevation around 1850 m a.s.l. (Fig. 1) Total discharge 5.4 ± 2.9 7.5 ± 1.2 5.5 ± 0.9 4.9 ± 0.9 11.2 ± 1.7 7.8 ± 1.2 10.7 ± 1.6 4.3 ± 0.7 6.8 ± 1.0 3.8 ± 0.6 (km 3 ) Peak discharge 1.  tative of the lowest elevation of the regional ice sheet margin indicates that p falls in the range 1.2-1.6 for the Kangerlussuaq catchment. An independent method using the elevations of six on-ice weather stations and their distance to the margin gives p = 1.3 ± 0.2, ignoring changes in catchment width with elevation. Plotting Fig. 4 with logarithmic axes confirms the relation between river discharge and air temperature with exponent p = 1.4. For the entire Greenland ice sheet we find that p = 1.5 (Fig. 8) for elevations up to 2550 m a.s.l., well above the runoff area. This p value indicates that the Kangerlussuaq catchment can be considered a representative section of the ice sheet in terms of its temperature sensitivity of meltwater release. At higher elevations the ice sheet converges towards its topographic peaks, altering the area-elevation distribution, yet not of relevance until runoff starts occurring much higher on the inland ice. We divide the ice sheet in four portions of equal surface area and find that p = 1.5 in the north (> 77.13 • N), p = 1.3 in the south (< 68.74 • N, including the Kangerlussuaq catchment), and p = 1.8 for the western slope of the ice sheet (Fig. 8). The largest hypsometric amplification is found for the eastern slope of the Greenland ice sheet (p = 2.4). This relatively high factor is likely caused by the eastern ice sheet being bordered by high mountains, forcing the ice sheet to converge into valley glaciers with a relatively small area at lower elevations, as opposed to the generally less irregular margin elsewhere in Greenland. Nevertheless, the hypsometric amplifier of 2.4 does suggest that increases in temperatures in east Greenland yield the largest increases in ice sheet meltwater discharge into the oceans.
To determine the impact of the hypsometric amplifier, we calculate using Eq. (4) how much temperature variability is required to produce a factor of 3 discharge variability (see previous section) for p = 1.4. Applying this temperature variability to a reference scenario with a linear surface slope, we find that meltwater release from the Kangerlussuaq catchment is 56 % larger for p = 1.4 than for p = 0. Extrapolating this methodology to other regions of the ice sheet, we find a hypsometric amplification of 62 % for the entire ice sheet (p = 1.5) and a 115 % amplification for the eastern slope of the Greenland ice sheet (p = 2.4). Figure 4 features the estimated temperature response of Watson River discharge using the p value of 1.4 and a discharge factor F D for bare ice (July-September) of 0.31. Applying a mid-range free-atmospheric adiabatic lapse rate λ of −7 × 10 −3 • C m −1 , we can solve Eq. (4) to find melt factor F M = 8.2 mm water equivalent • C −1 day −1 , which can be considered the catchment-average positive degree day factor. This factor is only of relevance to the gap filling of the Watson River discharge time series; it does not influence this study's main conclusions. Note that melt is estimated from atmospheric temperature outside the shallow, stable, surfacecontrolled boundary layer over the ice sheet, which is argued Figure 7. Hourly Watson River discharge (black), and hourly (light blue) and daily (dark blue) catchment-total ice sheet surface meltwater runoff. The red line gives meltwater runoff with elevation-dependent routing delay. Dark grey lines are river discharge estimated from temperature (see Methods section). Light grey areas represent discharge uncertainty. Plus symbols give the timing of precipitation events in Kangerlussuaq exceeding 1 mm in 6 h. to be preferable because it better represents atmospheric forcing of melt (Lang, 1968;Ohmura, 2001).
The Watson River discharge time series confirms the value of the hypsometric amplifier in the Kangerlussuaq catchment and thus that variability in meltwater release from the ice sheet is disproportional to the variability in atmospheric forcing of melt. Yet the uncertainty in determining river discharge from air temperature (i.e., not the topic of this study) remains large at a conservative 70 % (Fig. 4). Part of the scatter in Fig. 4 is due to the fact that meltwater generated at the ice sheet surface takes time to transit the supra-, en-, sub-, and proglacial routing environments, which is not properly accounted for in the figure. Therefore we turn towards our second science question: how to quantify such routing delays.

Routing delays
Watson River discharge trails ice sheet surface meltwater runoff by several days (Figs. 6a and 7) caused by routing delays in transiting the supra-, en-, sub-, and proglacial environments. Figure 9 illustrates in grey the average routing delay per elevation bin for which the ice sheet surface meltwater best matches river discharge in terms of correlation (r = 0.90-0.91). Routing delays increase with elevation, more so at higher elevation where the distance to the Watson River bridge measurement site becomes increasingly large due to the widening of the elevation bins. The dependency of the delay (t d in h) on surface elevation (z in m a.s.l.) can be approximated (r = 0.997) by the polynomial t d = 65.6 · 10 −6 · z 2 − 18.1 · 10 −3 · z + 7.3.
This delay also includes the transit through the 26-33 km long proglacial river system, which we estimate to typically take 4-5 h based on measurements taken during a 2010 ice-dammed lake drainage  and is thus small compared to the total delay, which exceeds 1 week at 1800 m a.s.l. We calculate the average effective, horizontal travel velocity through the supra-, en-and subglacial environments to be 0.7-0.8 km h −1 for water originating from the majority (> 90 %) of the runoff area (650-1850 m a.s.l.). This spatially rather constant water velocity indicates that for the multi-year average time delay as presented by Eq. (5), hydraulic efficiency is similar throughout most of the runoff area of the ice sheet. Our velocity values are below the average but within the 0.2-1.6 km h −1 range reported by Chandler et al. (2013) and Cowton et al. (2013), who used a methodology that bypasses supraglacial routing. Mikkelsen et al. (2016) found a general discharge delay of between 1 and 5 days, which we confirm for the ice sheet area between 650 and 1450 m a.s.l., making up a dominant portion of 58 % of the runoff area. The routing delay Eq. (5) is used in Fig. 7 to adjust the hourly catchment-total values in light blue to delayed runoff in red. We apply an additional 10 h smoothing to adjust the meltwater runoff record from being a composite of discrete elevation bins to something that also represents the spread in routing delays within the 100 m elevation bins. Delayed runoff matches river discharge in terms of amplitude and variability, especially during the low discharge years 2009, 2013, and 2015. However, during certain periods agreement is lower. For instance, in 2012 the delay in ice sheet runoff is too small before the peak of the melt season and too large during and after, likely related to a rapid development of the englacial drainage system (Bartholomew et al., 2011;Palmer et al., 2011). To discover intra-seasonal changes in routing delays and thus hydraulic efficiency of the ice sheet, we also apply the correlation procedure to each June-September month in the 2009-2015 period. Figure 9 illustrates that the elevation-dependent delays display a large range, with delays during some years being 2-3 times larger than during other years. June routing delays are, almost exclusively, larger than the multiyear average (Fig. 9a) because of slow percolation through winter-accumulated snow on the ice sheet surface (Bøggild et al., 2005) and an underdeveloped englacial drainage system . In July, commonly the peak river discharge month (Table 1), routing delays and the spread therein are smaller as surface snow is largely melted and presumably the englacial drainage system develops rapidly in response to increases in water supply. The reduced delays are most relevant at the lower and mid-elevation bands from which most meltwater originates. Our results suggest that development of the englacial drainage system can occur over the course of mere days; for instance in the first half of July 2012, the drainage system shifted from below-average efficiency (larger delays in Fig. 9a) to above-average (smaller delays in Fig. 9b). After the peak of the melt season, in August, Fig. 9 suggests that the englacial drainage system remains capable of efficiently routing the dropping water vol- Figure 9. Optimal meltwater routing delays for the supra-, en-, sub-, and proglacial environments determined for June (a), July (b), August (c) and September (d). As a reference, grey illustrates the values for the entire period with overlapping time series of river discharge and ice sheet surface meltwater runoff (June 2009-September 2015), as used for the red line in Fig. 7. The dashed line illustrates optimal routing delays for the 9-15 July 2012 period, encompassing the 10-14 July extreme discharge event.
umes given the fact that delays are typically similar to those in July (Fig. 9c). In September, routing delays increase, presumably as drainage channels close and hydraulic efficiency reduces, most notably at lower and mid-elevation where hydraulic efficiency is rapidly lost as water supply diminishes (Fig. 9d).
We note that not all optimal solutions at monthly time intervals represent glacial drainage delays accurately because both the ice sheet runoff and river discharge time series include features of other processes. For instance, the August and September 2010 delays are unrealistic due to a modeled melt underestimate during extreme melt periods  and an ice-dammed lake drainage that cannot be captured by the model. Likewise, August and September 2014 delays are inaccurate due to overestimated ice sheet runoff during a rain event (see Discussion section). We therefore regard the monthly panels in Fig. 9 as ensemble solutions and do not overemphasize results from specific months. However, some outliers are realistic, such as the July 2012 delay at elevations above 1200 m a.s.l. when hydraulic efficiency was higher than in any other year, as a consequence of unprecedented melting even at the highest plotted elevation bands (Mikkelsen et al., 2016).

Record-setting discharge July 2012
The highest discharge measured at Watson River occurred during the period 10-14 July 2012. Presumably this was the highest discharge in nearly 60 years given the 1955 Watson bridge road dam washout on 11 July 2012. This high discharge event coincided with a large melt episode impacting the entire Greenland ice sheet (Nghiem et al., 2012). In determining the causes of the extreme discharge, in addition to high ice sheet melting, we look into the timing of event. The KAN weather stations (Fig. 1) reveal that regionally the high-melt episode started around 11:00 UTC on 8 July and lasted 3.5 days. The Watson River reached its peak stage 3.3 days after the start of high melt. Calculating optimal discharge delays for the period 9-15 July, encompassing the extreme discharge event, we find that in 3.3 days the ice sheet meltwater could have traveled from as high as 1500 m a.s.l. (Fig. 9b). In other words, all meltwater passing at the bridge during the road dam washout originated from the multi-annual ablation area (defined in Van de Wal et al., 2012). Only towards the end of the 10-14 July extreme river discharge event did ice sheet contributions from the upper runoff area (∼ 1800 m a.s.l.) arrive at the bridge site. Therefore we consider the role of largely impermeable ice layers in the firn of the high-elevation accumulation area (Machguth et al., 2016) to be minor in causing peak river discharge and road dam washout (Mikkelsen et al., 2016).
In our interpretation, it appears relevant for the peakdischarge event that uninterrupted high melt persisted for a period of more than 3 days , allowing the (hypsometrically amplified) meltwater totals generated in the upper ablation and lower accumulation area on 8 and 9 July to exit the ice sheet along with that generated closer to the margin 2-3 days later (Fig. 9b). Figure 7g illustrates that transit times for the upper elevations were shorter than in other years. We speculate that this is partly because the 2012 peak melt was preceded by another high-melt episode around day of year 175 (dark blue line in Fig. 7g), increasing the efficiency of the englacial drainage system, allowing for faster transit during the following high-melt episode.
Altogether, we find the keys to record-setting discharge in the Watson River to be (1) intense ice sheet melting that (2) continues for several days, (3) amplified by ice sheet hypsometry, and (4) is preceded by another high-melt episode. Because of the multi-day transit time for meltwater to reach the community of Kangerlussuaq, an early warning system of future bridge flooding could in principle be deployed, relying on in situ ice sheet melt measurements.

Stage-discharge relation and road dam washout
The washout of the road dam at the Watson River bridge during the July 2012 peak discharge event potentially presents a challenge for deriving a single stage-discharge relation that applies to water levels recorded before the event and after, when the newly formed third channel became a bridge segment (Fig. 2). For water levels below ∼ 7 m over our lowest pressure transducer (Fig. 3), channel 3 remains dry and no change in the rating curve is expected, but for higher river stages and thus discharge values exceeding ∼ 1200 m 3 s −1 a change may have occurred. We have no indication that water stage at our pressure transducer site 150 m upstream of the bridge is affected by the road dam washout, but it cannot be ruled out, mostly at extreme discharge levels that could have yielded upstream pooling in years before the washout. If a change in the rating curve due to the road dam washout occurred, we would expect overestimated river discharge at high river stage pre-2012. Judging from the 2010 and 2011 comparison of discharge with ice sheet surface runoff in Fig. 7, overestimated discharge at high stage is at least plausible. An argument for keeping a single rating curve for the entire 2006-2015 period is that although the curve is established using ADCP data retrieved after the road dam washout, the calculated discharge during pre-washout icedammed lake drainages agrees well with values reported by Russell et al. (2011). It is likely that the 15 % uncertainty in our discharge calculation encompasses most or all of the changes inflicted by road dam washout and bridge segment construction. In any case, changing the rating curve and thus the amplitude of the discharge signal does not impact our primary conclusions on hypsometric amplification and routing moderation.

Rapid lake drainages
Although the total supraglacial lake volume that is released in drainage events is on the order of a few percent of the annual discharge for the Kangerlussuaq region of the ice sheet, rapid lake water release can contribute tens of percents to instantaneous discharge, primarily in the early melt season Mikkelsen et al., 2016). For instance, a surge in river discharge resulting from clustered supraglacial lake drainage was suggested to have occurred during days of year 180-184 in 2010 (Doyle et al., 2013;Hasholt et al., 2013). In this period, we find river discharge to exceed modeled ice sheet runoff as the latter does not take into account rapid lake drainage (Fig. 7e). From our records, we find no further evidence of a clustered drainage event of similar magnitude in the period 2009-2015, although smaller river discharge spikes around day of year 180 in other years (e.g., 2008, 2009, 2011) may very well also be related to supraglacial lake drainage (Bartholomew et al., 2011). Mikkelsen et al. (2016) determined the contribution of supraglacial lake drainage to annual Watson River discharge to be negligible. Since supraglacial lakes are relatively abundant in the wide melt area of the Kangerlussuaq catchment, the impact of lake drainages on studies using our methodology would logically be smaller elsewhere in Greenland.
Rapid drainages of the ice-dammed lake at the margin of the Russell Glacier, the northernmost glacier within the Kangerlussuaq catchment, produce discharge peaks lasting under 2 days (Russell et al., 2011). At least five such jökullaups can be identified towards the end of the melt season in 2007, 2008, 2010, 2012, and 2013 (Fig. 7). They contribute less than 1 % to the annual-total discharge, yet they can yield high peak values. For instance, the 2007 jökullaup resulted in a higher peak discharge than all melt-induced values in 2006-2015, except during high-melt years 2010 and 2012.

Rain events
An exception to the otherwise good agreement between ice sheet surface runoff and river discharge in Fig. 7 occurs for the mid-July 2012 extreme discharge event when river discharge is seen to exceed runoff. Underestimating modeled ice sheet melting has been identified to be a potential issue during this unprecedented episode  and could therefore explain mid-July mismatch. On other occasions, ice sheet runoff exhibits distinct peaks not entirely captured by the river discharge measurements in Fig. 7, such as in July 2010 (days 205-208), September 2013 (day 248) and 2015 (day 247), and, most notably, in late August 2014 (day 231-235). These peaks coincide with model-generated rainfall across the entire elevation range. Rain events can yield intense ice sheet surface melting largely due to increases in long-wave radiative and turbulent heat fluxes (Doyle et al., 2015;Van Tricht et al., 2016). Meteorological measurements in Kangerlussuaq ( Fig. 7; Cappelen, 2016) and near the ice sheet margin  confirm precipitation in these periods, with only the 2015 event not exceeding a precipitation rate of 1 mm in 6 h in Kangerlussuaq (Fig. 7j), and with positive temperatures yielding the possibility of liquid precipitation over the ice sheet. Given that little snow is left on the ice sheet surface to retain water towards the end of the melt season, the modeled rain and enhanced meltwater can produce distinct runoff peaks (dark blue lines in Fig. 7). Delayed runoff values (red lines) exceeding river discharge following these rain events indicate that rainfall is likely overestimated during these often short-lived events, most notably in late August 2014. Figure 7 illustrates that during most precipitation events that generate spikes in ice sheet runoff, the model performs well in capturing these events without overestimating river discharge. A clear example of a heavy rain event that was modeled with accuracy occurred in late August 2011 as described in detail by Doyle et al. (2015) (Fig. 7f). Climatologically, however, Kangerlussuaq is arid in terms of precipitation due to blocking topography to the southwest (Van den Johansson et al., 2015), providing this study with the possibility to study routing delays in an environment where complications by rain are minimal.

Piracy between catchments
Disparity between river discharge and ice sheet runoff may also be related to transient behavior between adjacent catchments, driven by seasonal changes in basal water pressure. Lindbäck et al. (2015) found that at ice sheet elevations above ∼ 1200 m a.s.l. changes in the subglacial hydrology can lead to large shifts in the Kangerlussuaq ice sheet catchment boundaries during the melt season, also known as piracy between catchments. Since piracy impacts the meltwater running off from the upper ablation area and above, the effect of a catchment boundary shift on catchment-wide runoff is expected to be largest during high-melt periods when a substantial amount of meltwater is generated at high elevations. However, Lindbäck et al. (2015) find that the Kangerlussuaq catchment above 1200 m a.s.l. shifts north in its entirety when subglacial water pressure builds, causing only a small change in surface area and therefore little (∼ 10 %) increase in catchment-wide runoff. Yet it remains plausible that the mismatch between river discharge and ice sheet runoff during the high-melt seasons of 2010-2012, and particularly during the mid-July 2012 extreme melt episode, is in part explained by a temporarily underestimated catchment area at high elevations. Figure 7 illustrates that modeled ice sheet runoff exceeds the river discharge values by 100-200 m 3 s −1 during springtime of all years except 2010 when accumulation during the preceding winter was well below average (Tedesco et al., 2011). We attribute this to a model underestimate of meltwater retention in winter-accumulated snow. Also, meltwater storage (and increases thereof) in supraglacial lakes  is not calculated by the model. Both processes provide plausible explanations given the better agreement between river discharge and ice sheet runoff in summer and autumn, when winter snow has melted and most supraglacial lakes have drained. Since the mismatch is smallest for spring of low-accumulation year 2010 (Tedesco et al., 2011), it is most likely that the SMB model underestimates retention in winter-accumulated snow, possibly due to underestimated snow accumulation that does not get recorded by weather stations when it collects in crevasses and in between ice hummocks in the lower ablation area ( Van den Broeke et al., 2008).

Retention in snow and glacial lakes
In all, river discharge and ice sheet runoff agree at an unprecedented level in this study because we (1) use an improved, validated, observation-based time series of modeled ice sheet runoff, (2) study a large catchment area implying a reduced relative uncertainty in delineation (Lindbäck et al., 2015), (3) constrain the river discharge calculations with superior ADCP measurements, and (4) introduce meltwater routing delays. In light of the high level of agreement in terms of variability and quantity, and the fact that we are able to provide plausible explanations for periods of mismatch, we find no evidence of meltwater storage in the en-and subglacial environments in amounts that surpass the detection limit as set by our methodological uncertainties. Using similar methods, such retention has been suggested to be significant for the Greenland ice sheet in previous studies (Rennermalm et al., 2013;Overeem et al., 2015;Smith et al., 2015;Mikkelsen et al., 2016) with values of up to half the meltwater availability reported. Whereas changes in supra-and subglacial storage are known to occur in Greenland as seen from rapid lake drainages, they are reported to make up only a few percent of the annual discharge in the Kangerlussuaq region www.the-cryosphere.net/11/1371/2017/ The Cryosphere, 11, 1371-1386, 2017 Palmer et al., 2015), i.e., well within our uncertainty for ice sheet runoff and river discharge.

Conclusions
The Watson River in west Greenland drains a ∼ 12 000 km 2 sector of the ice sheet, where the altitude at which the climatological mass budget is in balance (SMB + refreezing = 0) has increased from ∼ 1550 m a.s.l. for 1990-2011 to ∼ 1800 m a.s.l. for 2009-2015. We calculate ice sheet runoff and river discharge for a 7-and 10-year period, respectively, using an improved, validated, observation-based ice sheet SMB model and an updated river stage-discharge relation constrained by newly available ADCP measurements. Interannual variability in ice sheet meltwater release is found to be large; for instance river discharge ranges from 3.8 in 2015 to 11.2 km 3 in 2010, a factor of 3 difference. With hypsometry known to be an amplifier of ice sheet runoff, we deduce that discharge D can be approximated using regional air temperature T 0 using D ∼ T p+2 0 . Here p is the hypsometric amplifier, determined to be 1.4 ± 0.2 for the Kangerlussuaq catchment of the ice sheet. For p = 1.4 we calculate ice sheet meltwater release to be amplified by ∼ 56 % due to hypsometry. We determine p = 1.5 for the entire Greenland ice sheet, with regionally higher values and thus higher climate sensitivity, such as a value of p = 2.4 for the eastern slope of the ice sheet.
Diurnal variability in river discharge (∼ 100 m 3 s −1 ) is found to be more than an order of magnitude smaller than the variability in ice sheet surface meltwater runoff for the Kangerlussuaq catchment. The difference in diurnal variability is a result of the time lag involved in routing meltwater from its origin at the ice sheet surface through the supra-, en-, sub-, and proglacial environments to reach the river monitoring site. Introducing time lags to ice sheet runoff as a function of elevation results in a good agreement with river discharge. Calculated optimal delays reveal considerable changes in ice sheet hydraulic efficiency throughout the melt season, with time lags smallest shortly after high-melt episodes that overwhelm and develop the englacial drainage conduits. On average, the routing delays can be approximated by t d = 65.6·10 −6 ·z 2 −18.1·10 −3 ·z +7.3, which for instance yields that meltwater generated at 1500 m a.s.l. takes 5-6 days to be released from this relatively arid sector of the ice sheet. An implication of this result is that melt episodes are more likely to cause overbank river flooding, such as at Kangerlussuaq in mid-July 2012 when they last several days. Finally, due to the close agreement between river discharge and ice sheet surface meltwater runoff after the inclusion of routing delays, we find no evidence of meltwater retention in the en-and subglacial environments beyond the detection limit set by our methodological uncertainties. Data availability. All data can be obtained by contacting the first author.
Author contributions. D. van As and B. Hasholt conceived the study; B. Hasholt, A. B. Mikkelsen, and D. van As monitored river discharge; M. H. Nielsen provided ADCP data; K. Lindbäck delineated the catchment; L. C. Liljedahl enabled the in situ monitoring of the ice sheet; D. van As calculated river discharge and ice sheet runoff and wrote the text with contributions by all coauthors.
Competing interests. The authors declare that they have no conflict of interest.