Seasonal Variability in In-situ Supraglacial Streamflow and Drivers in Southwest Greenland in 2016

Greenland ice sheet surface runoff is evacuated through supraglacial stream networks, which influence surface mass balance as well as ice dynamics. However, in-situ observations of meltwater discharge through these stream networks are rare. In this study, we present 46 discharge measurements and continuous water level measurements for 62 days spanning 13 June to 13 August 2016 for a 0.6 km supraglacial stream catchment in southwest Greenland. The result is an unprecedented long record of supraglacial discharge capturing both diurnal and seasonal variability. By comparing in situ hydraulic geometry 15 parameters with previous studies, we find that significant heterogeneity exists such that estimating stream discharge using these parameters over ungauged supraglacial catchments could lead to substantial errors. A comparison of surface energy fluxes to stream discharge reveals shortwave radiation as the primary driver of melting (78% of melt energy). However, during high melt episodes, the shortwave only contributes to 50% of melt energy. Instead, the relative contribution of longwave radiation, sensible and latent heat fluxes to overall melt increases by 16.5%, 4%, and 7% respectively. Our data also show a 20 seasonal variation in the timing of daily maximum discharge during clear sky days, shifting from 16:00 local time (i.e., 2.75 hours after solar noon) in late June to 14:00 in late July, then rapidly returns to 16:00 in early August coincident with an abrupt drop in air temperature. These changes in peak daily flow timing can be attributed to a changing effective catchment area, resulting in a smaller stream network supplying water to the outlet at the end of the season throughout the melt season. Further work is needed to uncover how widespread rapid shift in the timing of peak discharge is across Greenland supraglacial streams, 25 and thus their potential impact on meltwater delivery to the subglacial system and ice dynamics. https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Mass loss from the Greenland ice sheet increased six-fold from the 1980s to the 2010s (286±20 Gt yr -1 ) (Mouginot et al., 2019;Shepherd et al., 2019;Velicogna et al., 2020). This mass loss was dominated by enhanced surface melting and runoff 30 . The increase in runoff raised Greenland's contribution to global sea-level rise from less than 5% in 1993 to more than 25% in 2014 (Chen et al., 2017). Increased surface melting also influences ice sheet basal properties (Das et al., 2008;Colgan et al., 2011;Bell et al., 2017;Kingslake et al., 2017;Flowers, 2018), and ice dynamics ( Van de Wal et al., 2008;Shepherd et al., 2009;Schoof, 2010;Hoffmann et al., 2011;Hewitt, 2013;Andrews et al., 2015). Despite the importance of Greenland's surface runoff to surface mass balance and sea level rise, only a handful of studies using in situ 35 supraglacial stream discharge to characterize current conditions (Holmes, 1955;Chandler et al., 2013;McGrath et al., 2011;Smith et al., 2015;Gleason et al., 2016;Smith et al., 2017), and these studies are limited to short periods (a couple of days to a maximum of 15 days).
Every melting season, numerous supraglacial stream/river drainage networks develop on the surface of the Greenland ice sheet ablation zone 2017;Yang and Smith, 2013;2016;Pitcher and Smith, 2019). These networks 40 transport runoff sourced from melting ice, snow, and/or slush within the stream catchment (Holmes, 1955;Karlstrom et al., 2014), and often terminate in moulins, from where meltwater moves within and beneath the ice sheet before emerging in proglacial rivers, lakes, fjords, and the ocean (Chu, 2014;Rennermalm et al., 2013).
Supraglacial stream discharge varies seasonally in concert with surface melting, with low flow in the beginning and end of the melting season and higher flow in the middle of the melting season (Holmes, 1955). Supraglacial discharge also 45 shows a pronounced diurnal variation (McGrath et al., 2011;Smith et al., 2017;Yang et al., 2018). While daily maximum discharge varies with catchment size and day of the season, discharge decreases when melt energy drops off at night (Marston et al., 1983;Mernild et al., 2006;McGrath et al., 2011;Yang et al., 2018). Diurnal variability and timing of meltwater delivery to the subglacial drainage system have been shown to influence ice sheet velocities in several studies (Bartholomew et al., layer of ice that retains meltwater temporarily and promotes subsurface flow (Karlstrom et al., 2014;Cooper et al., 2018), and therefore may slow the transport of meltwater to supraglacial streams (Munro, 2011;Karlstrom et al., 2014;Cook et al., 2016). 60 In Greenland's ablation zone, the seasonal and interannual variability in meltwater production, thereby surface runoff, is primarily driven by the variability of the absorption of shortwave radiation (Van den Broeke et al., 2011). Secondary melt drivers are turbulent fluxes of sensible and latent heat, particularly in the lower ablation zone of southwest Greenland (Van den Broeke et al., 2011, Fausto et al. 2016. Additional drivers are anomalously moist and warm air masses advected over the ice sheet by atmospheric rivers (Mattingly et al., 2018), and clouds with contrasting feedback to surface melt (Bennartz et al., 65 2013). While an increase in cloud cover enhances downward longwave radiation and hence melt (Van Tricht et al., 2016;Gallagher et al., 2020;Izeboud et al., 2020), it also limits the shortwave radiation thus decreasing summer melt in ablation areas (Hofer et al., 2017;Izeboud et al., 2020). Numerous studies examine the linkages between surface energy balance, surface melting, and runoff using regional climate models Noel et al., 2018), and automatic weather station data (van As et al., 2012), but relatively few compare surface energy fluxes with in situ observations of supraglacial stream 70 discharge in Greenland (Smith et al., 2017). Smith et al. (2015) and Gleason et al., (2016), examine how Greenland supraglacial streams can be characterized with hydraulic geometry principles. In terrestrial hydrology, at-a-station hydraulic geometry theory holds that channel width, depth, and velocity co-vary nonlinearly with discharge for a fixed stream cross-section (Leopold and Maddock, 1953;Gleason et al., 2015). This theory provides a set of equations with parameters that can be generalized to estimate discharge in ungauged rivers 75 (Smith et al., 1996;2015;Ashmore and Sauks, 2006;Andreadis et al., 2020). Similarly, a generalization of supraglacial streams' hydraulic geometries would open possibilities for scaling and modeling discharge. For example, Smith et al. (2015) used field-calibrated hydraulic geometry to estimate instantaneous discharges in 523 moulins in southwestern Greenland, yielding values ranging from 0.36 to 17.72 m 3 s -1 with a mean value of 3.15 m 3 s -1 . In contrast, Gleason et al. (2016) andSmith et al. (2017) argued that unlike terrestrial systems, uniform hydraulic behavior cannot be expected from an ice 80 substrate. However, only a few studies have quantified hydraulic geometry of supraglacial streams.
In this study, we present a 62-day time series of supraglacial streamflow, in southwest Greenland, spanning most of the 2016 melting season in terms of total meltwater production. The supraglacial stream drainage network was mapped using unmanned aerial vehicle (UAV) imagery and manual GPS observations, and surface energy fluxes were calculated or measured using meteorological observations from a nearby automatic weather station. Using these data, we examine supraglacial 85 discharge uncertainties, covariance between discharge and surface energy balance drivers, seasonal changes in daily peak discharge timing, and hydraulic geometry parameters. Finally, we compare our hydraulic geometry parameters to previous work on supraglacial streams. https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License.

Study Area
The study area is a 0.6 km 2 internally drained supraglacial catchment in southwest Greenland, hereafter called '660 90 catchment' (after "Point 660" where a gravel road ends at the ice sheet margin). The catchment is located ~1-2 km upstream of the ice edge between two outlet glaciers, Isunnguata Sermia and Russell Glacier, and roughly 35 km east of the town of Kangerlussuaq (Fig. 1). The stream network terminates in a moulin (67.1562°N, 50.0064°W). About 850 m upstream of the moulin, a gauging station for monitoring water level and discharge was installed at 67.1573°N, 49.9951°W. Elevations in the catchment span from 610 m near the gauging station to 660 m (above the WGS84 ellipsoid) at the catchment's highest point 95 ( Fig. 1).
The catchment surface consists of a rugged bare-ice landscape with small supraglacial ponds and an incised stream network (Fig. 2a). The clean bare-ice surface with minimal impurities has an albedo of 0.57±0.04  and had a thin (~0.1-0.3 m) surface layer of weathering crust comprising porous ice and cryoconite holes (Fig. 2c). These cryoconite holes are partially filled with water and accumulate "cryoconite," consisting of dust, sediment, and biological matter 100 (Takeuchi et al., 2000;Cooper et al., 2019). Cryoconite deposits are widespread in streams and ponds throughout the catchment. The catchment is situated in a region where winter snow accumulation is relatively low, and that experiences extensive melting from June through August so that little to no snow cover remains on the bare ice early in the melting season (Rennermalm et al., 2013;Ryan et al., 2019). stage. Stream water stage was measured every 5 min using a setup of two Solinst loggers (pressure transducers), a Levelogger in a weighted steel enclosure resting on the stream bed tied to a pole (Fig. 2b) and a Barologger (Solinst, 2020) installed 25-30 m northeast of the gauging station. Stage is calculated after barometric pressure correction. 115

120
Discharge was calculated with the velocity-area method using inputs of cross-sectional area and stream water velocity (e.g. Herschy, 1993a). In total, 46 observations of velocity and cross-sectional area were made, including 27 measured every hour from 15:30 on 26 July 2016 to 17:30 on 27 July 2016 (local time) to capture the entire diurnal range.
The remaining 19 observations were collected over the entire study period and sampled on average every 3-7 days between 12:00 to 17:00 local time. Cross-sections of stream depth were measured at 0.2 m intervals across the 1.7-3.3 m wide 125 stream. Though measurements were collected in the same location throughout the season, continuous thermal erosion of the bed resulted in small changes in cross-sectional geometry (Fig. 3a). Stream velocity was measured at each 0.2 m interval at 60% of the depth, with either a General Oceanics current meter or Price Type-AA current meter.
The discharge rating curve was generated with a best fit power-law (e.g. Herschy, 1993b): (1) 130 where p and β are constants estimated by fitting the curve to observations of discharge (Q) and water level, also called stage height (H), and α is the water level sensor offset from the stream bottom. In this study, the box with the level logger was placed on the stream bed and therefore α=0. The rating curve (Q = 3.925 H 2.44 , R 2 = 0.94) was then used to generate continuous discharge values from stage 140 measurements recorded every 5 min throughout the season. These data were in turn averaged to yield hourly discharge data ( Fig. 4).
Therefore, we conclude that our rating curve is sufficiently robust to estimate discharge. Four uncertainty estimates were calculated (see Appendix A for more details): 1) uncertainty for the 46 individual discharge measurements, 2) uncertainty https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License. due to the rating curve, 3) uncertainty of daily mean discharge, and 4) uncertainty due to the stream bed incision into the ice 145 over the melting season. While the uncertainty due to measurement errors (Ume) was estimated as 10.8%, uncertainty due to the rating curve (URC) was estimated at 17%. The measurement uncertainties were encompassed in the envelope of uncertainty due to the rating curve (Q±URC) (Fig. 4). The averaging of hourly discharge to daily mean discharge generated an uncertainty of 25% (Xdm). Uncertainty due to the stream bed incision was estimated using the cross-sectional profiles of the stream bed ( Fig. 3). 150

Figure 4: Rating curve (black line) determined from the best fit of power law (Eq. 1) to observations of stage and discharge (blue dots). Error bars show measurement error uncertainty (Ume). Rating curve uncertainty (URC) is shown in the grey shaded area (see Appendix A for uncertainty calculations).
Reconstructing hydrographs for supraglacial streams with high seasonal and diurnal variations with a rating curve is 155 typically unreliable (Smith et al., 2017;Pitcher et al., 2019). In terrestrial rivers, shifts in the rating curve are a reflection of either a datum adjustment or changes in channel cross-section. Unlike terrestrial rivers, the bed under supraglacial streams is constantly melting and incising into the ice resulting in an ever-changing cross-sectional profile. To examine if our rating curve is robust despite channel cross-sectional profiles changes, we compared depth profiles collected while velocity was measured (Fig. 3). By separating profiles collected over high flows through the season from profiles collected over the 26-hr 160 period of hourly measurements and assuming small incision over the 26-hour-period, depth measurement errors can be isolated from incision errors. While profiles collected over the season show a 3.7 cm standard deviation (Fig 3a), hourly https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License.
profiles collected over the 26-hr period show a 1.9 cm standard deviation (Fig 3b). The uncertainty in stream discharge (here we use the 95% confidence interval) due to non-uniform stream bed incision and depth measurement errors, Uin is 10.9 % (of the average stream depth), and the depth measurement alone is 5.9 %. Despite these errors, the channel geometry incises 165 uniformly through the season with all the cross-sections lying inside the uncertainty levels ( Fig. 3a).

Calculation of hydraulic geometry parameters
Hydraulic geometry parameters were calculated to examine the relative importance of width, velocity, and depth in controlling discharge, and to compare with other studies reporting hydraulic geometry data for supraglacial streams. Hydraulic geometry is a set of equations that describes the relationship between discharge and stream channel geometry, namely, width, 170 depth, and velocity (Leopold and Maddock, 1953): where, w, d, and v are stream width, depth, and velocity of the cross-section, respectively. The exponents b, f, and m represent 175 the slopes of the power laws. The magnitude of the exponents represents the rates of change of each variable with respect to the independent variable, discharge. The coefficients a, c, and k represent the y-axis intercepts. The law of conservation of mass implies that the product of coefficients a, c, and k should equal unity. Likewise, the sum of the exponents b, f, and m should also equal unity (Leopold and Maddock, 1953).

Automatic weather station observations in the 660 catchment 180
To identify the timing of daily maximum melt during clear sky days, a shortwave pyranometer (HOBO S-LIB-M003) was mounted at 2 m height on an automated weather station (AWS) ~20-25 m from the gauging station. Other meteorological data were measured at this station and are a part of the dataset accompanying this paper. However, better data continuity and accuracy is provided by the nearby AWS, KAN_L, located ~7-8 km southeast of our study area. Therefore, meteorological data from KAN_L were used in our analysis. The KAN_L AWS is maintained by the Geological Survey of 185 Denmark and Greenland (GEUS) (van As et al., 2011).

Surface energy balance model
To examine surface energy drivers of supraglacial discharge, energy balance components were obtained from a surface energy balance model (described in . This model uses forcing data from in-situ meteorological and radiative observations from KAN_L to calculate the surface energy balance components net shortwave radiation, net longwave radiation, sensible heat flux, latent heat flux, sub-surface conductive heat flux, and heat flux from rain. While incoming shortwave and longwave radiation are gathered from the AWS, turbulent heat fluxes are calculated from nearsurface gradients of meteorological variables, air temperature, humidity, and wind speed using Monin-Obukhov similarity theory.

Catchment delineation 195
The catchment boundary and supraglacial stream network were manually digitized using two sources. Firstly, we used WorldView-1 (WV1) panchromatic imagery (spatial resolution of 1 m) acquired on 16 August 2016 to manually digitize the stream network. We also collected 20,000 handheld GPS points of the catchment boundary in the field, by walking along the catchment divide separating the stream network from adjacent catchments. We did not observe a change in catchment size during the study period. We estimate the catchment area to be accurate within 5% given that it was manually identified in the 200 field. However, the precise delineation of the catchment is not relevant to the outcome of the study.

Results
Our catchment has a dendritic drainage pattern (Strahler order = 4), and is internally drained, meaning that all surface meltwater is routed through streams, tributaries, and ponds to a terminal moulin (Fig. 1). Repeated visits to the study site between 2013-2019, suggest that this network, with streams ~ 0.1-15 m widths and ~ 0.1-2 m depths, is an annually recurrent 205 feature of the local hydrological system. A new tributary joining the network was observed in 2017 but does not affect the present study.

Hourly and daily variations in supraglacial discharge
Stream discharge between 13 June and 13 August varies strongly both diurnally and seasonally (Fig. 5a). Hourly discharge fell as low as 0 m 3 s -1 at night and daily peaks exceeded 0.3 m 3 s -1 on most days. Three distinct melt episodes with 210 larger discharges were recorded on 13 June (0.81 m 3 s -1 ), 19 June (0.94 m 3 s -1 ), and 16 July (0.93 m 3 s -1 ). During these times peak flows were almost double the long-term average of daily maximum discharge (0.5 m 3 s -1 ) and occurred around 15:00 local time.
Daily maximum discharge varies from 0.05-0.94 m 3 s -1 through the season, with the highest values around the three melt episodes (Fig. 5b). Daily minimum discharge has much less seasonal variability than daily maximum discharge but 215 exhibits two occurrences with anomalously larger flow on 23 June and 19 July (Fig. 5b). During the melt episodes, these positive anomalies in daily minimum discharge follow a steep decrease in daily maximum discharge, meaning, the anomalously large low flows at night follow a dip in day-time streamflow. Compared to the first and second episodes, daily minimum discharge behaves differently after the third episode. After the second episode, and before the third episode begins, night-time low-flow discharge occasionally falls as low as 0 m 3 s -1 but rises to 0.04 m 3 s -1 after the third episode. Finally, the 220 diurnal amplitude (daily maximum minus daily minimum) tracks daily maximum discharge except for the second and third https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License. melt episodes with large daily minimum discharge (Fig. 5c). After the third melt episode, there is a steady decline in diurnal amplitude from 0.64-0.33 m 3 s -1 (21 July-13 August). The daily mean discharge varies from 0-0.51 m 3 s -1 over the 62 days (Fig. 6a) and co-varies seasonally with daily maximum discharge except for the timing of the peak daily mean discharge (Fig. 5b). The daily mean discharge peaks on 19 July, three days after the second-largest melt episode in hourly discharge. In contrast, the second largest peak in daily mean 230 flow occurs on 20 June, which is the same day as the largest episode in hourly discharge. In both cases, hourly maximum discharge is accompanied by several days of very high daily minimum flows (Fig. 5b), which explains the discrepancy between the timing of the daily mean and daily maximum episodes.

Surface energy balance
Throughout the season, net shortwave radiation exceeds all other surface energy fluxes and thus is the primary driver 240 of stream discharge (Fig. 6b). However, the second and third melt episodes coincide with peak longwave radiation (65 Wm -2 increase) and turbulent heat fluxes (40-80 Wm -2 increase) along with a drop in shortwave radiation (110-120 Wm -2 decrease) ( Fig. 6). Thus, during high-melt episodes, longwave radiation and turbulent heat fluxes become more pronounced drivers of streamflow. Among all energy fluxes, sensible heat flux correlates most with daily mean discharge (R=0.88, p-value = 0.00001). During the third melt episode, the hourly peak discharge coincides with a peak in shortwave radiation on 16 July 245 (Fig. 7a). However, the peak daily mean discharge occurs three days later on 19 July 2016 due to high net longwave radiation and turbulent heat fluxes from 16-20 July (Fig. 6). This high net longwave radiation was caused by overcast conditions (Fig.   7b) and resulted in a high low flow at night (Fig. 7a). The consistently large low-flow persisted from 17 to 20 July (Fig. 7a) and resulted in the seasonal peak discharge on 19 July (0.51 m 3 s -1 ) (Fig. 6a). This can also be seen in the hourly variation of https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License. surface energy balance components and stream discharge (Fig. 7a). Between 17-20 July, night-time streamflow is much higher 250 than before and after (Fig. 7a), and coincides with increased net longwave radiation. While a dip in shortwave radiation on 18 July decreases the high flow during the day, the low flow during the night increases due to a spike in net longwave radiation (Fig. 7a). To further examine the difference in energy balance components' contribution to high versus low stream discharge, we aggregated components for the second and third melt episodes, and compared to all-season data (Fig. 8). The shortwave radiation proportion of melt energy fell from 78 % to 50 % (normalized) during the melt episodes. Simultaneously, the 260 contribution of longwave radiation and turbulent heat fluxes increased during those days. The longwave radiation's proportion of melt energy increased from a seasonal average of -0.32 to -0.085 during the peak flow days (increased by 16.5%), while sensible and latent heat fluxes' proportion of melt energy increased by 0.06 (4%) and 0.11 (7%) during the melt episodes, respectively. https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License.

Timing of peak daily discharge
To examine if the transport of meltwater from its production on the ice sheet surface to the discharge observation site 270 varies over the season, we calculated time to daily maximum discharge, following 'time-to-peak' methodology in traditional terrestrial hydrology (Chow, 1964). As the season progresses, the time to daily maximum discharge will reflect changes in network storage and transport efficiency. In contrast, during clear-sky days, when solar radiation drives melt, the timing of daily maximum surface meltwater creation is not expected to change over the season and be proportional to solar noon. While the incoming solar radiation peaks at solar noon around 13:15 local time, the timing of the daily maximum discharge varies 275 from 16:00 in late June to 14:00 in late July (Fig. 9a). The time lag between the solar and discharge peaks changes from 3 hours to 1 hour from 30 June to 31 July and has a statistically significant negative trend (R 2 = 0.793, p-value = 0.00011). After 31 July, the time lag abruptly shifts back to initial season conditions and stabilizes at 3 hours in early August. This shift in the lag coincides with the sudden decrease in daily mean temperatures from 4.3 °C on 31 July to 2.5 °C on 3 August (Fig. 9b).
Though the daily mean temperature started to decrease after the third melt episode on 19 July, the sudden drop in early August 280 coincides with the shift in time of daily maximum discharge due to a drop in night-time temperatures close to 1 °C (Fig. 9b).
These air temperature measurements were collected at 2 m above the surface and therefore the skin temperatures are expected to be below freezing.

Hydraulic Geometry
Hydraulic geometry parameters are determined by generating a power law between stream discharge, and width (R 2 290 = 0.87, p-value = 0.00001), depth (R 2 = 0.94, p-value = 0.00001) and velocity (R 2 = 0.88, p-value = 0.00001) (Eq. 2-4 respectively). For the 660 catchment, the exponents b, f, and m are 0.194, 0.387 and 0.373 respectively, and coefficients a, c, and k are 3.44, 0.541, and 0.632 respectively. While the sum of the exponents of these power laws, representing sensitivity of discharge to the individual variable, equals 1, the product of the coefficients must also equal 1 (Leopold and Maddock, 1953). In our study the sum of exponents equals 0.95 and the product of coefficients equals 1.17.
The oldest study we are aware of, Holmes (1955), reports supraglacial stream discharge of 0.14-5 m 3 s -1 (about five times larger than at the 660 catchment) at a catchment of 25-50 km 2 (much larger than the 660 catchment) with a daily maximum discharge occurring between 16:00-20:00 local time in southwest Greenland. In a study by Smith et al. (2017) in a much 310 larger (63 km 2 ) supraglacial catchment, daily maximum discharge occurred between 18:00-20:00 local time, and discharge varied between 5-26 m 3 s -1 (Table 1). In all studies providing time-series of stream discharge (Marston, 1983;Mernild et al., 2006;McGrath et al., 2011;Chandler et al., 2013;Smith et al., 2017), the lag between solar noon and daily maximum discharge is larger for the larger magnitude of stream discharge, due to runoff generation from a larger catchment area, thus increasing the distance of surface routing and a later peak of daily discharge at the catchment outlet (Smith et al., 2017). 315 In addition to the diurnal variability, supraglacial stream discharge exhibits strong seasonal variability at the 660 catchment. We observed three melt episodes of high discharge, one peaking on 13 June (350% of mean discharge), a second on 19 June (400% of mean discharge), and a third on 16 July (400% of mean discharge) (Fig. 5a). The timing of these episodes corresponds with periods of anomalous high river discharge observed ~35 km downstream in the Watson River, Kangerlussuaq (van As et al., 2018). The 660 catchment is a very small subcatchment of the Watson River catchment (12000 320 km 2 ).
Over the 62-day study period, while net shortwave radiation provides the majority of melt energy and is the primary driver of streamflow, net longwave radiation and turbulent heat fluxes (sensible and latent) become more dominant melt drivers during the three melt episodes (Fig. 8). These findings agree with Greenland-wide studies identifying a link between longwave radiation and enhanced surface melting (Van Tricht et al., 2016;Gallagher et al., 2020) while disagreeing with studies 325 suggesting that overcast conditions and lower solar radiation reduce surface melt in the ablation zone (Hofer et al., 2017;Izeboud et al., 2020). Furthermore, we find that daily average sensible heat flux has the highest correlation with stream discharge (R 2 =0.83, p-value=0.00001) of all energy balance components (e.g. net shortwave radiation, R 2 =0.23, p-https://doi.org/10.5194/tc-2020-314 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License. value=0.035). Correlating energy components with stream discharge is justified here due to the quick routing in this catchment (1-3 hours). This contribution of sensible heat flux is consistent with Fausto et al. (2016), who show that peak melting occurs 330 at times with anomalously large turbulent energy fluxes.   and Gleason et al. (2016) have common data sets. Range of width and discharge are the same for both the studies. However, Gleason et al. (2016) primarily discussed the hydraulics of these streams. Therefore, this data set is mentioned as Gleason et al. (2016) in the discussion.
We find that the time of maximum diurnal discharge varies from 16:00 in late June to 14:00 in late July (Fig. 9a), 335 which translates to a lag between maximum diurnal melt and maximum diurnal discharge from 1-3 hours. This lag is consistent with the previous studies, where it varies from 3-9.5 hours (Holmes, 1955;Mernild et al., 2006;McGrath et al., 2011;Chandler et al., 2013;Smith et al., 2017). The difference in lag is most likely due to the difference in the size of catchment. While large catchments in Holmes (1955), Mernild et al. (2006), Chandler et al. (2013) and Smith et al. (2017) show lags between 3-9.5 hours, smaller catchments in Marston (1983) and McGrath et al. (2011) show lags between 2-4 340 hours (Table 1). For the 660 catchment, the lag decreases through the season from 3 hours in late June to 1 hour in late July, consistent with the change in lag through the melting season shown by Mernild et al. (2006) from 5-7 hours in May to 3-4 hours in August. Furthermore, at the 660 catchment, in early August, the lag abruptly increases to the initial season conditions and stabilizes at 3 hours, coinciding with a sudden drop in air temperature from 4.3 °C to 2.5 °C from 31 July to 3 August (Fig. 9b). With the night-time air temperatures close to 1 °C, the skin temperatures could be below freezing affecting 345 the efficiency of the stream network, thereby increasing the temporary storage of melt generated in the weathering crust, and delay the transport of melt to the moulin. This phenomenon of a shift in the timing of daily maximum discharge agrees with the model study by Yang et al. (2018), who shows that a poorly developed stream network has delayed timing of peak and up to half the diurnal variation in stream discharge relative to a well-developed and an efficient network.
Our work confirms the findings by Gleason et al. (2016) that hydraulic geometry parameters cannot be generalized 350 for supraglacial rivers in Greenland, despite having a common ice substrate. Comparing our data with parameters from previous studies (Knighton, 1981;Marston, 1983;Gleason et al., 2016) in a ternary diagram reveals three clusters (Fig. 10). These three clusters can be grouped based on their b-values (width power-law exponent). The first cluster has high b-values (b>=0.35) and includes the downstream station of Gleason et al. (2016), which are smaller streams with discharge varying between 0.006 to 0.402 m 3 s -1 . The second cluster has low b-values (b<=0.05) and includes at-a-station data from Gleason et al. (2016), which 355 are larger streams with discharge varying between 4.58 to 23.12 m 3 s -1 . Finally, the third cluster has moderate b-values (0.05<b<0.35), and includes this study, and Knighton (1981) and Marston (1983). Though the discharge from Knighton (1981) is two orders smaller than ours, Marston (1983) has similar discharge values (Table 1). The streams with discharge of same order of magnitude i.e., varying between 0-1 m 3 s -1 , from Knighton (1981), Marston (1983), and the current study show moderate sensitivity to stream width (moderate values of b) and streams with higher magnitude of discharge show very small 360 sensitivity to stream width (smaller values of b). However, small streams from Gleason et al. (2016)   Our analysis of a season-long record of streamflow and its drivers has several implications for large-scale Greenland ice sheet hydrology. First, we confirm Gleason et al. (2016)'s work by finding that hydraulic geometry parameters vary across ice sheet catchments. Thus, using generalized parameter space to estimate supraglacial streamflow from ungauged catchments may not be valid. Second, we find that longwave radiation and turbulent fluxes have an increased contribution (adding up to 370 50% of melt energy), in governing stream discharge during the melt episodes. Given that several regional climate models underestimate turbulent fluxes, they also will underestimate melt episodes (van den Broeke et al., 2011;and Fausto et al., 2016). Since the most accurate way to estimate surface runoff from the entire Greenland ice sheet is through regional climate/surface mass balance models Mernild et al., 2018;Noel et al., 2018), underestimating turbulent heat fluxes also underestimates runoff. Third, with the lack of lengthy records of supraglacial stream discharge over the 375 Greenland ice sheet, the 62-day long time series could expand climate/surface mass balance models validation capability to a great extent. Currently, the only study to validate model simulated runoff used field observations of streamflow covering just 3 days (Smith et al. 2017). Lastly, the accurate timing and magnitude of peak discharge can help understand its influence on subglacial water pressure and ice velocities, although the impact of stream velocities are probably minor from this catchment since it is located very close to the ice sheet margin.

Conclusions
We present one of the longest records of Greenland supraglacial stream discharge, spanning 62-days of the 2016 melting season for a 0.6 km 2 catchment in southwest Greenland. These long-term observations could be used in validating regional climate models, currently the best tools to estimate surface runoff from the entire Greenland ice sheet. Concurring with the previous studies, we show that hydraulic geometry parameters cannot be generalized for supraglacial rivers in 385 Greenland. The observed stream discharges vary both seasonally and diurnally. Our record includes three distinct episodes of large discharge, one from 13-14 June, the second one from 19-25 June, and the third one from 15-21 July. The daily maximum discharge and amplitude show similar seasonal and diurnal variation except during the second and third melt episodes when nighttime melting reduced the daily amplitude. During the third melt episode, the large daily discharge minimums, caused by high longwave radiation, drive the seasonal peak discharge on 19 July with continuous high flow in 390 the day and night for 3-4 days (16-19 July). The stream discharge is primarily driven by net shortwave radiation through the melting season, except during the high melt episodes when net longwave radiation and turbulent heat fluxes show an increased contribution (16% and 4-7% respectively) to melt energy. The lag between the time of daily maximum discharge and time of daily maximum melt (assumed to be solar noon here) varies through the season from 3 hours in late June to 1 hour in late July and goes back to 3 hours in early August due to a sudden drop in air temperature. This change in lag could 395 be due to the expansion and contraction of the stream network, caused by probable freezing skin temperatures at night, affecting network efficiency and therefore, change in routing time of the melt. Though this theory was explained by a model simulation from Yang et al. (2018), further work is required to reveal if the rapid shift in the timing of peak discharge, we notice at the 660 catchment, also takes place across Greenland supraglacial streams and analyze the influence of stream network development on timing and magnitude of peak discharge. 400

Appendix A: Streamflow uncertainty
The uncertainty of the 46 individual discharge observations due to seven different types of measurement errors in velocity and stage height was calculated following Herschy (2002) and WMO (2010), and assumes that segment discharges and standard uncertainties are approximately equal in each of the segments in the cross-section: 405 ( ) = √ 2 + 2 + 1 { 2 + 2 + 2 + 2 + 2 } (A1) where Ume is uncertainty due to measurement errors, Q is discharge, K is the so called coverage factor (K=2 gives the 95% confidence interval), um is uncertainty in determination of mean velocity for number of verticals (M), us is uncertainty due to calibration errors, ub is uncertainty in width, ud is uncertainty in depth, up is uncertainty in determination of mean velocity for number of points in the vertical, uc is uncertainty in determination of mean velocity for current meter rating, and ue is 410 uncertainty in determination of mean velocity for time of exposure. A summary of all nomenclature used in this study is found in Table A1.  (Table A2), we find Ume = 10.8%. 415 The uncertainty in hourly stream discharge due to the rating curve was determined with a statistical method by calculating the standard error of estimate, Se, where the quadratic rating curve was linearized with a logarithmic transformation following Herschy (1994): where t is Student's t correction (for 95% confidence t=2), n is number of discharge measurements, Qi is discharge measured, 420 and Qc is discharge estimated using a rating curve.

Uncertainties
Values used in this study um 4.5 % for a minimum of 10 verticals us 1 % ub 0.15 % (for width range 0 -100 m) ud 0.65 % (for depth range 0.4 -6 m) up 7.5 % (number of verticals = 1) uc 1 % (for average velocity of around 0.25 ms -1 ue 4% (for time of exposure between 30 -60 s) M 10 (number of verticals varied between 10-16 K 2 for 95% C.I Using Equation A2, the uncertainty of hourly stream discharge due to the rating curve, URC is calculated to 17%. 425 Lastly, the uncertainty in daily mean discharge due to the averaging of hourly data is estimated in two steps using the methodology of Dymond and Christain (1982). First, Smr, the standard error of the daily mean, is calculated using a logarithmic transformation of discharge, Q from Equation 1, in order to make it a linear relation: where, N is the number of data points in the sample (here, 24 samples in a day). Second, Xdm is uncertainty in the daily mean 430 discharge is calculated as: where, Swl is the standard error of log of water level measurement (calculated using Equation A2). Using Equation A3 and A4, the uncertainty in the daily mean discharge, Xdm is calculated to 25 %. All uncertainties are expressed as the 95% confidence level. 435

Author contribution statement
RM performed data analysis and wrote the manuscript with support from AR. RM and AR conceived and planned the study.
SL and MC did the majority of field data collection, with support from RM, AR, and SC. DvA performed surface energy balance modeling. All authors discussed the results and contributed to the final manuscript.

Data availability
KAN_L weather station data from the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) and the Greenland Analogue Project (GAP) were provided by the Geological Survey of Denmark and Greenland (GEUS) at http://www.promice.dk. All other data will be available at the PANGEA repository (in the meanwhile this data is made available to the reviewers as supplementary material). 450