Intra-seasonal variability in supraglacial stream sediment on the Greenland Ice Sheet

On the surface of the Greenland Ice Sheet, the presence of low-albedo features greatly contributes to ablation zone meltwater production. Some of the lowest albedo features on the Ice Sheet are water-filled supraglacial stream channels, especially those with abundant deposits of consolidated cryoconite sediment. Because these sediments enhance melting by disproportionately lowering albedo, studying their spatial extent can provide a better understanding of Greenland’s contribution to global sea level rise. However, little is known about the spatial distribution of supraglacial stream sediment, or how it changes in response to seasonal flow regimes. Here, we surveyed a supraglacial stream network in Southwest Greenland, collecting imagery from seven uncrewed aerial vehicle (UAV) flights over the course of 24 days in 2019. Using Structure-from-Motion-generated orthomosaic imagery and digital elevation models (DEMs), we manually digitized the banks of the supraglacial stream channels, classified the areal coverage of sediment deposits, and modeled how the terrain influences the amount of incoming solar radiation at the Ice Sheet surface. We used imagery classified by surface types and in-situ spectrometer measurements to determine how changes in sediment cover altered albedo. We found that, within our study area, only 15% of cryoconite sediment was consolidated in cryoconite holes; the remaining 85% was located within supraglacial streams mostly concentrated on daily inundated riverbanks (hereafter termed floodplains). Sediment cover and stream width are highly correlated, suggesting that sediment influx into supraglacial drainage systems widens stream channels or darkens previously widened channels. This reduces albedo in floodplains that already receive greater solar radiation due to their flatness. Additionally, the areal extent of stream sediments increased in August following seasonal peak flow, suggesting that as stream power decreases, more sediment accumulates in supraglacial channels. This negative feedback loop for melting may delay Greenland’s runoff to the latter end of the melt season. This study shows in unprecedented detail where and when sediment is deposited and how these deposits potentially impact the Ice Sheet surface energy balance. These findings may allow for better prediction of how supraglacial floodplains, and the microbiomes they contain, might change in response to increased melting.


Introduction
Greenland Ice Sheet mass loss is a major contributor to global sea level rise. Between 2002-2017 Gt yr −1 to sea level rise (Zou et al., 2020), corresponding to 46% of total Arctic land ice contributions (Box et al., 2018). Over half of the observed ice loss results from negative surface mass balances (Shepherd et al., 2020). This mass loss process is exacerbated by the prevalence of low albedo features on the Ice Sheet surface including supraglacial streams and sediment Ryan et al., 2018). Supraglacial stream drainage networks form in the ablation zone of glaciers around the world (Holmes, 1955;Ferguson, 1973;Marston, 1983;Smith et al., 2015;Watson et al., 2016;Kingslake et al., 2017;Ryan et al., 2018;Lu et al., 2020). These supraglacial stream networks form a series of supraglacial lakes, meandering rivers, and floodplains (flat, sediment-covered regions along the stream banks that are often diurnally inundated with meltwater) as meltwater flows toward moulins (Rennermalm et al., 2013a;Karlstrom et al., 2013;Chu, 2014;Leidman et al., 2021a).
On the Greenland Ice Sheet, supraglacial streams cover between 2 and 7% of the ablation zone surface Lu et al., 2021). Despite their relatively low spatial coverage, supraglacial streams explain 12% of Greenland's albedo variability in the ablation zone . One reason for this is that supraglacial streams consolidate low albedo (0.09) sediment along their beds Leidman et al., 2021b). Sediment is transported onto the surface of the Greenland Ice Sheet by windborne transport from local proglacial moraines and valleys (Bøggild et al., 2010;Nelson et al., 2014;Humbert et al., 2020) as well as, but to a lesser extent, accumulation of fine dust particles from long-range sources such as Asia (Bory et al., 2002;Bory et al., 2003), surface melting of ice containing outcroppings of Pleistocene era long-range dust (Bøggild et al., 2010;Cooper et al., 2018), landslides (Svennevig, 2019), and thrust faulting (Moore et al., 2010). Once on the Ice Sheet, sediment stays distributed on the surface, melts into cryoconite holes, or washes into supraglacial streams (Leidman et al., 2021a).
Surface impurities from algae and sediment can be widespread in regions such as the so-called Dark Zone in Southwest Greenland resulting in algae-rich regions having albedo values~0.4 lower than clean ice (i.e., ice with a minimal amount of surface impurities) regions Yallop et al., 2012;Stibal et al., 2017;Williamson et al., 2018;Wang et al., 2020). As more of the Dark Zone becomes snow-free for longer periods due to climate change, algal blooms growing on the surface also increase and accelerate melting (Wang et al., 2018;Ryan et al., 2019;Wang et al., 2020). Distributed sediments on the ice surface often consolidate and, due to their dark coloration, absorb solar radiation and melt into the ice, forming cryoconite holes roughly 0.15-0.30 m deep (Bøggild et al., 2010). Sediment within cryoconite holes harbors cyanobacteria, algae, and micro-invertebrates (Uetake et al., 2016;Uetake et al., 2019). As these organisms grow, they wrap around sediment grains and form extracellular polymeric substances causing grains to clump together into granules (Hodson et al., 2010;Nagar et al., 2021). As a result, cryoconite can have organic contents of >6% dry weight at higher Ice Sheet elevations where cryoconite holes have more time to develop (Stibal et al., 2010).
Once cryoconite holes form, surface albedo may be 0.6 higher than areas with distributed sediment as the bottom of the cryoconite hole is largely shaded from direct sunlight due to the hole geometry and ice lids coving the holes after nighttime freezing (Bøggild et al., 2010). However, sediment granules can flush out of cryoconite holes during heavy melt events (often associated with cloudy conditions or rainfall; Takeuchi et al., 2018;Tedstone et al., 2020;Muthyala et al., 2022) and deposit in supraglacial streams (Leidman et al., 2021b). The similar granular structure and composition of sediment found in supraglacial streams also indicates that it predominantly originates from these cryoconite holes (Bøggild et al., 2010).
Supraglacial water bodies (lakes, streams) have a significantly lower albedo than the surrounding clean ice (0.21 for water bodies compared to 0.55 for clean ice) and, as such, disproportionately contribute to negative surface mass balances (Pope et al., 2016;Ryan et al., 2018). Previous research using synoptic numerical models show that for shallow streams (0.5 m depth), 38% of the incoming spectrally integrated solar radiation is absorbed by the bed compared to 51% by the water column (Bray et al., 2017). As a result, for terrestrial streams, dark sediment along the streambed can facilitate enough energy absorption to increase water temperatures by over 2°C (Bray et al., 2017). Isenko et al. (2005) observed that this process also occurs in supraglacial streams resulting in stream temperatures up to 0.4°C when sediment was present. Energy absorbed by supraglacial stream sediment may be partially re-radiated back to the atmosphere or exported from the glacier's terminus before contributing to increased melting. However, the long transport distance between the absorption and export locations and the similar hydraulic properties of these systems to terrestrial rivers (Marston, 1983;Smith et al., 2015;Gleason et al., 2016;Muthyala et al., 2022) suggests that much of the solar energy absorbed by sediment goes towards increasing melt rates. The highly turbulent flow characteristics of supraglacial streams (Gleason et al., 2016) would also facilitate a more efficient transfer of energy to the ice for melting. The fact that supraglacial streams display extensive incision compared to the background ablation rates (Ferguson, 1973;Karlstrom et al., 2013;Leidman et al., 2021a) also points to an efficient transfer of energy to the bed. While this incision can lead to shadowing of the stream channel (Leidman et al., 2021b), this effect is only significant in areas of high relief (such as near immediately upstream of moulins). Despite this incision induced shadowing, hydrological modeling of supraglacial stream meandering indicates that the majority of the melt energy goes towards lateral rather than vertical melting thereby continuing the exposure of the stream bed to direct solar radiation (Karlstrom et al., 2013). In this study, we therefore assume that all energy absorbed due to changes in albedo translate to changes in meltwater production.
This study investigates how sediment within supraglacial stream networks change throughout the melt season and postulate how it might impact albedo in response to changing meltwater inputs. We use repeat photography from uncrewed aerial vehicle (UAV) flights to classify supraglacial streams and sediment in a supraglacial catchment in Southwest Greenland. In total, seven classified imagery mosaics were used to determine how sediment cover changes throughout the stream length and melt season. Additionally, we use digital elevation models (DEMs) derived from the UAV imagery and in-situ hyperspectral measurements of various surfaces to determine the amount of incoming solar radiation exposure for each surface class; bare ice, water, stream sediment, and cryoconite holes. We compare results with stream discharge to examine whether sediment cover is dependent on Frontiers in Earth Science frontiersin.org 02

FIGURE 1
Map of the study site with UAV imagery from 22 July 2019 clipped to the large area overlap footprint. Orange lines show the location of manually digitized shorelines and the pink line indicates the centerline of the main channel. The location of the stream gauge used to measure water depth throughout the observation period, the furthest downstream location of the August 7th observed main channel, and the moulin for the main channel are shown as a yellow triangle, a green star, and a red circle respectively. The inset map shows the location of the field site in southwestern Greenland. Background imagery is a Maxar Technologies image obtained from Google Earth collected on July 24, 2012 overlain on the ArcticDEM (Polar Geospatial Center et al., 2018). Because the time offset between the background image (2012) and the data collection (2019), the movement of the ice has shifted the channel location and the Aug 7 Endpoint and moulin are outside the 2012 channel.
TABLE 1 Statistics of the spatial distribution of sediment and water within each flight scene including the date (a) and average time the images were collected (West Greenland local time) (b), UAV platform (c), area of the UAV scene footprint (d), structure-from-motion output cell size (e), the percent of the scene covered by river water rivers (includes both river beds covered by sediment and clear ice) (f), the percent of the scene covered by bare ice (g), the percent of the scene covered by sediment (includes sediment on bare ice and in the river) (h), and the percent of the scene where rivers are covered by sediment. Albedo was calculated as the weighted average of the extent of the scene covered in river water (f minus i) and bare ice (g) using and albedo value of 0.26 and 0.55 respectively from Ryan et al. (2018) as well as the extent of the scene covered in sediemnt using an albedo of 0.07 derived from Figure 2 stream power and, if so, discuss how sediment might impact the input of meltwater and nutrients to the subglacial environment.

UAV imagery collection
To investigate the intra-seasonal variation of supraglacial stream sediment distribution, we conducted repeat UAV flights with a DJI Phantom 3 or a DJI Mavic 2 Pro over a single supraglacial catchment in Southwest Greenland (67.15⁰N, 50.00⁰W) (Figure 1). In total, seven flights were conducted (July 22, August 2, 4, 7, 11, 14, and 15, 2019) (Table 1). All imagery was collected within 3 h of solar noon. The DJI Phantom 3 was mounted with an 18 megapixel FC3610 camera whereas the Mavic 2 Pro was mounted with a 17 megapixel L1D-20c camera. All images were collected at a maximum height of 60 m from the surface and had at least 80% image overlap. Markers were placed around the study site to use as ground control points including orange plastic trays (0.35 x 0.25 m) (Supplementary Figure S1), orange painted rocks (~0.1 m diameter), and red tarps with black X's to demark the center point (1.5 x 1.5 m).

UAV imagery processing
The imagery collected from each of the seven UAV flights was georeferenced and processed using Agisoft Metashape Pro v1.6.5 Structure-from-Motion software to create high quality dense point clouds that were then converted to DEM and orthomosaic images. Outputs were exported at their maximum raster resolution (0.026 ± 0.016 m for the DEMs and 0.008 ± 0.004 m for the orthomosaics; Table 1). All DEMs and orthomosaics were georeferenced to the UTM 22 N projection, a conformal, shapepreserving projection with minimal distortion in this area. Due to logistical challenges of collecting data on the Ice Sheet, imagery coverage varied from 0.007 to 0.151 km 2 (Supplementary Figures S2-S8; Table 1).

Surface reflectance data
To relate changes in the extent of sediment and stream channels to changes in surface energy balance of the field site, we created an empirical relationship between water depth, sediment cover, and broadband reflectance collected with an Analytical Spectral Devices (ASD) Fieldspec4 Pro spectrometer. The ASD contiguously samples the spectral range 350-2500 nm for each measurement. In total, 41 ASD measurements were collected, all of which were under consistent illumination conditions within 2 h of solar noon. Sampling sites were selected quasi-randomly throughout the area of overlapping flight footprints (Supplementary Figure S9) in order to measure a range of different surfaces with varying sediment coverage and water depths. The spectral reflectance at each location was calculated by taking the ratio between the average of ten upwelling solar radiance measurements of the ground surface and the average of ten upwelling solar radiance measurements of a Spectralon panel reflector. Spectral irradiance was TABLE 2 Statistics for the coverage of the streams within the UAV scenes derived from manually digitized shorelines (b) as well as the relative coverage of clean ice (c) and sediment within (d) and outside of (e)observed in each flight scene. Stream albedo was calculated as the weighted average of columns c and d using albedo values of 0.24 and 0.07 respectively.

FIGURE 2
Scatterplot of reflectance values for each of ASD measurement and the sediment cover derived from classified handheld imagery of the ground surface within the ASD field of view. Colors indicate the water depth of the measurement site. The least-squares trendline was used to determine the albedo of 0.07 for 100% sediment cover where y is reflectance, and x is the percent sediment cover. R 2 is the coefficient of determination of the trendline.
Frontiers in Earth Science frontiersin.org 04 used to calculate the spectrally weighted broadband reflectance value. Co-located ground surface images were taken in tandem with ASD measurements. Surface images were collected in triplicate with a RICOH WG-30W 5-megapixel digital camera.
Surface images were cropped using a centered circular geometry to match the area covered by the 7.5°viewing angle of the spectrometer at the height of the camera during image capture. Cropped images were then converted to black and white images in Adobe Photoshop v22.5.1. All images were manually inspected to ensure that they matched the expected spectrometer coverage and areas of contiguous sediment patches that displayed speckling or glare artifacts were manually edited. The proportion of black versus white pixels was used to determine the average sediment cover for the triplicate images. Sediment cover could be obtained this way due to the high visual difference between the clean ice (converted to white) and dark sediment (converted to black). All images were taken within 2 h of solar noon under similar cloud conditions. Clean ice and dark sediment were spectrally quite distinct and we determined that any errors from slight tilting of the camera were likely significantly higher than any errors associated with illumination. Reflectance was plotted against sediment cover and fitted with a least-squares best fit exponential trend line ( Figure 2). The 100% sediment cover intercept of this relationship was used to determine the reflectance of sediment observed in UAV imagery, which was found to be 0.07 ( Figure 2).

Water-level and runoff data
Stream water pressure was measured with a Solinst Levellogger pressure transducer placed in a RACO-2 metal project box filled with rocks with conduit holes knocked out to allow water to easily flow into the box. The box was anchored to a PVC pole, which was drilled into the thalweg of the main stream channel (denoted as the stream gauge in Figure 1). The Levellogger continuously measured water pressure every minute between July 13th and August 10th, 2019 (Supplementary Figure S10). Water level was calculated by correcting water pressures for atmospheric pressure. Atmospheric pressure was measured using a Solinst Barologger fixed in place 2 m above the ice on the shore of the largest channel within the study area (hereafter termed the main channel). Manual measurements of water depth were recorded with a meter stick periodically throughout the field season to determine if any measurement drift was present and to correct for systematic instrumentation error. The metal conduit box was thick enough to prevent accelerated melting of the underlying bed faster than the surrounding ablation rate, hence the logger stayed stationary relative to the bed throughout the measurement period. Additionally, stream discharge was measured with a Sontek Flowtracker2 Acoustic Doppler Velocimeter immediately upstream of the Levellogger. In total, 33 measurements of discharge were recorded throughout the study period and used to produce a rating curve to convert water level measurements to discharge throughout the observation period (Supplementary Figure S11).
Additionally, we examined runoff estimates gathered from the Modéle Atmosphérique Régional (MAR) regional climate model version 3.11.5 (Fettweis et al., 2017) in order to extend the runoff record to the length of the UAV observation period. Average modelled daily runoff over the course of the observation period was calculated from the four 15 km grid cells closest to our study area (Supplementary Table S1). The center points of the grid cells were 10.8 ± 2.2 km (mean ± standard deviation) from the 7-Aug endpoint, had an average elevation of 627 m and average ice coverage of 74%. Runoff values were highly correlated with stream discharge measurements (p<<0.01, Supplementary Figure  S11). This data was used to compare changes in reflectance of the stream system to changes in runoff (Figure 7).

Feature mapping and network analysis
The shorelines for all streams at least 0.1 m wide from all seven flights were manually digitized using ArcGIS Pro v2.7.2. While flights were taken at different times of day (Table 1) and different water levels, analysis of stream water level showed that water heights were all within 30% of each other. Importantly, all of the images were collected when water level was high enough to fully inundate the observed floodplains. Channel walls were relatively steep outside of the floodplain boundary and therefore the variability in water level for the different image acquisition times resulted in insignificant variability in stream width and only a minimal impact on our analysis of sediment extent. Channel centerlines were found by calculating Euclidean distance from the shorelines and running a flow accumulation algorithm on the distance raster to the 7-Aug end point. The centerlines were inspected manually and proved to be more accurate than those generated by the ArcGIS Pro's centerline tool.
The generated centerlines were used to determine how sediment cover on the stream bed changed along the length of the main channel. We quantified sediment cover for each zone. Zones subdivided the stream area every 1-m along the stream centerline from the 7-Aug endpoint (the furthest observed downstream location) to the August 11th start-point (the furthest observed upstream location). To create these zones, we generated points every 0.2 m along the stream centerlines. Then, the Network Analyst Closest Facility tool in ArcGIS was used to find the distance along the centerline of each point to the endpoint for the corresponding imagery from that flight date. The path distances determined for each point were adjusted so that they all represented the distance to the 7-Aug endpoint. We then used these points to create 1m zones that spanned the width of the stream channel and extended 1-m along the centerline of each UAV scene. Finally, we applied an inverse distance weighting (IDW) algorithm on the network distance point values. The resulting interpolated output raster, clipped to the stream extent of the hand digitized stream shorelines, was then used to find 1-m contours. These contours were used as the boundaries for the 1-m zones used for quantifying sediment cover along the stream length.
Additionally, the 1-m zones were used to analyze how elevation, stream width, water-surface slope, and exposure to sunlight changed as you moved up the main stream channel. Average elevation of each 1-m zone was calculated using zonal statistics on the UAV derived DEM. Elevation data was also extracted to each centerline point and a linear regression was fitted to each point contained within each 1-m zone to determine the water surface slope (Supplementary Figure S12). The average width of the channel for each zone was determined by creating a Euclidean distance raster from the channel shorelines and doubling the resulting extracted values at the center-point locations. Average width Frontiers in Earth Science frontiersin.org 05 for each 1-m zone was therefore calculated as twice the distance to the shorelines averaged over all the center-points within the 1-m zone.

Surface feature classification
The glacier surface observed in each UAV image was classified into three types: bare-ice, sediment, and water bodies, where sediment could cover either ice or stream. The classification was achieved by first applying a Maximum Likelihood supervised classification algorithm in ArcGIS Pro to each orthomosaic image to identify stream water, bare-ice, and sediment. Twentyfive training polygons were used for each of the three classes of each image totaling 600 polygons averaging 1.0 m 2 in area each. This is the same method for classifying sediment areas as published by Leidman et al. (2021a) and, since each classification was made using training polygons from the same image, this method effectively accounts for variable solar illumination of the different mosaic images. The location of each classified pixel was then analyzed to determine whether it was found within or outside of the manually digitized stream extent. Visual inspection showed that while the classification of sediment areas was highly accurate, this process led to an overestimation of water outside of streams as well as pixels within the streams classified as bare-ice. This was likely due to the spectral similarities between bare-ice and water in RGB imagery as well as variable illumination throughout the scene due to shadows and scattered clouds. To account for this error, we used the manually digitized stream extent to convert all pixels identified as bare-ice within a stream channel to water pixels and all pixels classified as water outside of the stream as bare-ice. This correction needed to be applied to roughly 30% of stream pixels and 7% of non-stream pixels but resulted in a classification that better matched observations. Zonal statistics were then used to determine the average sediment concentration for each 1-m zone along the channel centerlines.

Influence of terrain on solar irradiance
In our study area, the walls of deeply incised supraglacial stream can cast shadows the Ice Sheet surface and cause large variability in incoming solar radiation (Leidman et al., 2021b). To examine how shadowing from local topography affects the radiation exposure of sediment within stream systems compared to sediment within cryoconite holes, we used ArcGIS Pro's Area Solar Radiation tool (Kausika and van Sark, 2021) and the Structure-from-Motion DEM to model irradiance for each pixel within the study area over the course of a 24-h clear-sky day. Specifically, solar irradiance was calculated every 15 min for 24-h using the DEM retrieved on July 22nd as well as a transmissivity value of 0.545 which is representative for the atmospheric conditions in this area (van den Broeke et al., 2008). The distribution of solar radiation values for each surface cover class inside and out of the stream extent (i.e. sediment, ice, and water; section 3.2) were examined. This allowed us to determine if sediment inside the stream received more incoming solar radiation than sediment consolidated on bare-ice. Incoming solar radiation was also correlated with the average stream width of each 1 m zone, however the relationship is weak (Supplementary Figure S13).

Meteorological data
Daily averaged wind speeds were gathered from the PROMICE Automatic Weather Station network (Fausto et al., 2021) for KAN_L. KAN_L is located about 7 km south of our study site at 67.0955 N, 49.9513 W with an elevation of 670 m. Wind speed is measured at approximately 3.1 m above the ice with a R. M. Young 05103 wind monitor every 10 min. Trends in wind speed data were used to determine if local weather conditions potentially had an impact on sediment deposition rates within the observation period.

Stream geometry
In the study area, the main channel was 4.6 ± 2.6 m wide (mean ± standard deviation), with a maximum width of 17 m across (including diurnally inundated floodplains). Tributaries were, on average, one-third the size of the main channel with widths weighted by mapped stream length of 1.7 ± 2.2 m (mean ± standard deviation) wide with a median width of 0.8 m. There were no endorheic lakes within the study area besides a few small melt pools (less than 1-m wide). Instead, the supraglacial hydrology within the study area was dominated by streams. Supraglacial streams covered 8.4 ± 0.4% (mean ± standard deviation) of the large scene area (hereafter defined as the overlapping region of the four orthomosaics with the widest coverage, Figure 1 and Table 2). This is greater than the 2% reported by Ryan et al. (2018) and 1-2% reported by Bøggild et al. (2010) for northeast Greenland as well as the seasonal maximum supraglacial river area fraction of 4.5% reported by Lu et al. (2021) for the Devon Ice Cap.

Sediment coverage
Within the large area scene (overlapping stream area for the four UAV scenes with the largest coverage), sediment covered 26.5 ± 3.3% (mean ± standard deviation) of the stream bed (Table 1, range 23-31%). Sediment in streams on average comprised 85% of the areal coverage of all sediment visible to the UAV within the study area. In total, sediment covered 2.6% of the large scene area on average and, although not measured, the granule size of sediment within the streams appeared similar to that found in cryoconite holes. Sediment coverage along the main channel reached up to 91% cover within the 1-m cross-sectional contour zones ( Figure 3A). While sediment cover within the stream increased throughout the observation period, no consistent trend was observed as a function of distance to the endpoint. Instead, changes in sediment cover within the main channel ranged from -3.0 to 2.6% per day (standard deviation = 0.72% d −1 ) ( Figure 3B). Within the small scene area (overlapping stream area for all seven scenes), sediment concentrations decreased from July 22nd to August 4th and then steadily increased thereafter at a rate of 0.75% per day on average (10% in total) (Figure 4, Supplementary Figure S14). The change was smaller for the large scene area where the August sediment cover only increased by 0.09% per day. Sediment outside of the stream system (mainly cryoconite holes) within the large scene area showed a much Frontiers in Earth Science frontiersin.org 06 lower trend increasing in areal coverage by only 0.02% per day in August.
Sediment cover also seemed to have an impact on the roughness of the stream bed. Areas of total sediment cover were uniformly flat whereas areas of patchwork sediment cover displayed a rougher ice surface with hummocks and holes multiple centimeters deep where sediment congregated. While not directly measured, this geometry seemed constant throughout the observation period and likely did not have an effect on the seasonal change in sediment cover. Areas near the thalweg that were covered by water throughout the night had smoother bed surfaces than areas with diurnal inundation therefore this process likely had a minimal impact on the roughness of the streambed.

Relationship between width, sediment cover, and solar radiation
Analysis of the spatial distribution of sediment floodplains revealed that wider sections of the stream contained more sediment than narrower sections (p-value << 0.01, Figure 5, Supplementary Figure S15). As a result, stream sections with above average width were 21% covered in sediment whereas stream sections with below average width were only 2% covered in sediment on average. The relationship between stream width and sediment cover persisted throughout the observation period (Supplementary Figure   S15). A statistically significant relationship was also present between channel slope and sediment cover (p < 0.01), however the lower reliability of the elevation data from the structure-from-motion DEMs used to derive slope gives us limited confidence in this relationship (Supplementary Figure S12).
The wider widths and flatter terrain of sediment covered areas of the stream resulted in fewer shadows precluding incoming solar radiation from reaching the surface. As a result, sediment within the river system on July 22nd, on average, received 7.8% more incoming solar radiation compared to sediment outside of the stream ( Figure 6, Supplementary Figure S16).

Change in albedo over the melt season
Albedo of the study area was calculated using values for bare-ice and shallow water albedo from the literature (0.55 and 0.26 respectively, Ryan et al., 2018) as well as measured reflectance values for cryoconite sediment (0.07, Figure 2) weighted using the spatial coverage for each class from the classified imagery. The maximum albedo of the large scene area before peak melt (July 22nd) was 0.522. The trend in albedo after peak melting was determined by the linear regression of the albedo values of the three August flights. The trend in albedo before peak melting was determined as the linear regression between the July 22nd albedo and the August 1st albedo extrapolated from the aforementioned post-peak trend. After peak melt, the albedo decreased by 4.0 x 10 -4 per day, largely due to darkening of the stream channel (stream albedo decreased by 1.6 x 10 -3 per day) (Figure 7). If no sediment was present in the stream system, the July 22nd albedo would have been 0.53.

Supraglacial stream coverage
The greater coverage of supraglacial streams in our study is likely partly due to the higher imagery resolution (0.01 m) compared to previous studies (0.15 m used by Ryan et al. (2018) and 10 m used by  Figure 1). (B) The rate of change in sediment cover (slope of the linear regression of sediment cover, % d −1 ) for each 1-m zone of the main channel with at least three observations. Sections of the stream where the change in sediment cover was statistically significant (p<0.05), are indicated with red dots. On average, sediment cover of the main channel increased by 0.07% per day.

FIGURE 4
Sediment cover of the main channel between 582 and 757 m from the 7-Aug endpoint (c.f. Figure 1). This stream reach corresponds to the section of the main channel where all seven flights overlap. The blue line indicates the average daily water level recorded within the main channel. The black dots indicate the percentage of the stream channel covered by sediment for each flight.

Frontiers in Earth Science
frontiersin.org 07 Lu et al. (2021). Additionally, our study area is closer to the ice margin (2 km instead of~15-45 km for Ryan et al. (2018) and 0-120 km for Lu et al. (2021)) and is subject to higher runoff production rates (Mernild and Liston, 2012;van As et al., 2017) than those studied previously, potentially contributing to increased channel density. Our field site has a substantially lower ice velocity compared to Nioghalvfjerdsfjorden glacier in Northeast Greenland (Joughin et al., 2010). However, the spatial coverage of supraglacial streams was similar to values reported by Lu et al. (2021) in region of the Ice Sheet dominated by high velocity basal sliding. This indicates that regions near the ice margin can display similar supraglacial drainage properties as high-velocity ice streams potentially due to high melt and strain rates causing denser channel networks. This is likely due to the fact that high-velocity ice streams can contain extensive areas of compressive surface-parallel mean stress resulting in areas of ponding that are hydraulically isolated by creep closure and refreezing (Chudley et al., 2021).

Sediment dynamics
The supraglacial sediment cover observed in this study (26.5%) is similar to values found by Leidman et al. (2021a) who found sediment covered 24% of the streambed in the same region. To better understand the drivers of the observed sediment distribution and how topography affected the frequency of floodplains, a Fast Fourier Transform (FFT) analysis was conducted in (using the Microsoft Excel FFT function). The FFT analysis was done on a 1024 m section of the main channel (223-1247 m in Figure 3) using the average sediment cover for all seven scenes for each 1-m zone. The magnitude of the Fourier transform peaked at a recurrence interval of 128 and 512 m (Supplementary Figure S17), roughly 27 and 111 times the average width of the channel (4.6 m) respectively. While such a short section of stream is likely not long enough to fully determine the periodicity of sediment cover peaks along the main channel, it shows that sediment fluctuates at longer wavelengths than the 81 m that would be expected for terrestrial river systems based on the average stream width of our system (Williams, 1986). This indicates that meandering formation is not only determined by thermal erosion processes (similar to bank erosion in many terrestrial bedrock channels). Instead, large-scale topographic influences from subglacial undulations and crevassing have a significant impact on meander geometries. The impact of large-scale topography is supported by the fact that stream channels were predominantly oriented perpendicular to the glacier flow direction and therefore parallel to the direction of stress fractures (Supplementary Figure S18). Stress fractures likely cause areas of preferential erosion in this region (Scott and Wohl, 2019). As a result, the fracturing caused by subglacial topography is likely a major driver of channel geometry. Sediment, however, still can shape supraglacial stream channel geometry. This is evidenced by the fact that wider sections of streams had significantly more sediment than narrow sections ( Figure 5). This suggests that sediment is widening stream systems or darkening already wide channels. In both scenarios, sediment leads to a greater absorption of incoming solar radiation.
Within our study area, the spatial coverage of consolidated sediment was more than five times greater within stream channels compared to sediment within cryoconite holes (consolidated sediment classified outside of the digitized stream extent). Thus, it appears that in this region of the Ice Sheet, sediment is efficiently flushed from the ice and deposited within floodplains where it can more effectively absorb solar radiation. However, the proportion of sediment found in streams is likely overestimated due to the slight non-nadir viewing angle of the UAV when imaging cryoconites holes. Many studies assume that cryoconite is largely located within cryoconite holes and is therefore predominantly shaded or ice covered throughout the melt season (MacDonell and Fitzsimons, 2008). This assumption, however, likely underplays the importance of cryoconite sediments in changing ablation zone albedo.
The August increase in sediment cover (0.75% per day, Figure 4, Supplementary Figure S19) suggests that diminished stream power after peak melt facilitated sediment deposition within streams. Intraseasonal changes in supraglacial stream sediment cover is therefore likely controlled by stream power rather than changes in sediment loading.

FIGURE 6
Box plot of incoming solar radiation for sediment in the stream (right) and outside of the stream (left). The difference in the average values of the two distributions was 14 Wm −2 . This is equivalent to an increase in incoming solar radiation of 7.8% for sediment in the stream.

FIGURE 5
Relationship between stream width and average sediment cover for each meter of the main channel that falls within the small scene where all seven flights overlaps.

Frontiers in Earth Science
frontiersin.org 08 The changes in sediment cover inversely covary with seasonal changes in water level within the study area (Figure 4). Supraglacial stream discharge in this region has a strong seasonal cycle (Muthyala et al., 2022), and peaked on August 1st in 2019 ( Figure 4). The relationship between sediment cover and stream discharge therefore agree with sediment transport principles that state that greater flows have greater stream power and sediment carrying capacity. Thus, greater stream flow has a greater ability to scour sediment from the bed whereas lower flows allow for deposition to occur (Dingman, 2009). We did not observe any delay between changes in water level and sediment cover. This nearly instantaneous response indicates that supraglacial streams can likely supply sediment to the bed as needed as the critical shear stress decreases. Therefore, supraglacial streams in this region are likely not significantly sediment starved systems.
Aeolian sediment loading to the Ice Sheet was not likely responsible for the observed change in sediment cover during the study period. While daily average wind speeds varied between 2-7 ms −1 at the nearby KAN_L PROMICE station, trends during the study period were insignificant and slightly decreasing (-0.03 ms −1 d −1 , p = 0.28) (Supplementary Figure  S20) with no correlation with observed sediment cover. As such, sediment loading to the Ice Sheet was likely constant or slightly decreasing after peak melting. Additionally, the observed expansion of floodplains is probably not caused by the release of interstitial sediment trapped within the ice due to the decreasing runoff rates in August (Figure 7). Some of the lowest concentrations of sediment is found nearest the moulin (<500 m from the 7-Aug endpoint in Figure 3A) where the channel is more incised and steeper. Following hydrological principles, the stream reaches closest to the moulin therefore likely have the greatest stream power and the greatest potential for sediment scouring (Dingman, 2015).

Impact of albedo changes on surface melting
We estimate how our seasonal change in albedo relates to changes in surface melting by using findings from previous studies, specifically Konzelmann and Braithwaite (1995) and Riihelä et al. (2019). First, we examined Konzelmann and Braithwaite (1995) as their field sight more closely resembled our study site near the ice edge whereas Riihelä et al. (2019) was focused on the accumulation zone. Konzelmann and Brathwaite (1995) used measurements from nine ablation stakes in northeastern Greenland to determine a relationship between albedo and ablation energy. To extrapolate this relationship to our field site, we assumed that our study area is representative of the entire ablation zone (likely an over-prediction of supraglacial floodplain coverage). Based on the observations indicating that albedo decreased linearly throughout August, we found the linear regression of albedo during that period (Figure 7). We also assumed that the ablation zone covers 12% of the Ice Sheet as determined by Ryan et al. (2019). It is safe to assume that the large majority of the energy absorbed by the sediment eventually contributes to melting of ice since less than 8% of incoming radiation that reaches the streambed is reflected away before being absorbed by a shallow water column (Bray et al., 2017) and because runoff can be retained within the Ice Sheet for 1-6 months before export (Rennermalm et al., 2013b). Thus, there is plenty of time for heat to transfer from the water to the ice. Based on these assumptions, we use the Konzelmann and Braithwaite (1995) relationship to find that our observed decrease in albedo of 0.012 would cause an increase in melt energy of 2.6 W m -2 during the observation period. This translates to an increase in runoff of 4.1 Gt of water equivalence in August compared to earlier in the melt season (1.4% of annual Ice Sheet mass balance, Mouginot et al., 2019). The 4.1 Gt is likely an upper bound to the potential melt volume produced by sediment covered supraglacial floodplains. If we instead use the 2% stream cover reported by Ryan et al. (2018) is used for calculating albedo, sediment deposition within stream channels would result in 1.6 Gt of additional melting in August (0.5% of annual mass balance). Second, we used Riihelä et al. (2019) who derived surface albedo measurements from the CLARA-A2 dataset to relate large scale changes in albedo to runoff. Both methods, Konzelmann and Braithwaite (1995) and Riihelä et al. (2019), yielded similar results, albeit runoff volumes were slightly higher using the Riihelä et al. (2019) relationship (an increase of 4.6 Gt of water equivalence in August or 1.8 Gt assuming 2% stream coverage). These results suggest that decreasing stream flow after peak melting drives increasing deposition of sediment within supraglacial streams resulting in a more negative surface mass balance in the latter periods of the melt season. As supraglacial stream sediment transport capacity changes throughout the season, we postulate that sediment deposition is coupled with melting of the ablation zone creating a negative feedback loop. In this negative feedback loop, increased melt rates due to greater surface water coverage and frictional heating is counteracted by decreased coverage of stream sediment. Through that same process, as observed in this study, after melting peaks and surface water coverage begins to decrease, the resultant increase in albedo is dampened by deposition of sediment.

FIGURE 7
Albedo of the main stream channel for the overlapping area of the four largest coverage flights (i.e., large scene area). Albedo was calculated based on the percent coverage of sediment and water using, an albedo value of 0.26 for water , and an albedo value of 0.07 for sediment. A segmented line (where the second part of the line is a regression calculated with the equation shown above) has been added to support the hypothesis made in Section 5.2 that stream albedo increases before peak discharge as sediment is eroded and then decreases after peak discharge as sediment is deposited. Runoff values (blue stipple line) are daily averaged runoff values from the MAR v3.11.5 regional climate model.

Frontiers in Earth Science
frontiersin.org 09

Impact on ice dynamics
The inverse relationship between stream flow and sediment cover within supraglacial streams could potentially impact ice sheet dynamics. Meltwater routed through supraglacial streams is transported to the englacial and subglacial environments via moulins (Das et al., 2008;Andrews et al., 2014;Smith et al., 2015;Smith et al., 2017;Smith et al., 2021). After passing through the moulin, meltwater is routed through englacial channels to the bed (Catania and Neumann, 2010;Cowton et al., 2013). Water delivered to the bed then forms subglacial channel networks (Bartholomew et al., 2010;Chandler et al., 2013). Variability in meltwater supply to these subglacial channel networks can then have significant impacts on ice sheet dynamics (Fitzpatrick et al., 2013;Shannon et al., 2013;Stevens et al., 2016;Stevens et al., 2021;Tsai et al., 2021). The relationship between supraglacial stream inputs and ice dynamics is dependent on a number of factors such as moulin density (Banwell et al., 2016) and the distance from the terminus (Stevens et al., 2021), however, seasonal timing of discharge fluctuations also plays a major role in determining velocity responses . The change in sediment deposition described in this study would cause meltwater inputs to the bed to decrease in July and increase later into the melt season when subglacial channels are more developed and can more efficiently drain excess water (Bartholomew et al., 2010;Hoffman et al., 2011). As a result, the dynamic response to surface runoff might be dampened by supraglacial sediment scouring and deposition throughout the melt season. Due to the relatively small change in albedo and resultant meltwater inputs however, this effect is likely fairly minor.

Supraglacial stream sediment under a changing climate
With climate change, increasing spatial coverage of supraglacial streams on the Greenland Ice Sheet is expected to reduce ablation zone albedo (Ignéczi et al., 2016). Unprecedented atmospheric conditions caused exceptional melting during 2019, when this study was conducted . The intraseasonal changes in sediment cover observed in this study therefore might be more representative of future dynamics than current averages. However, increased surface melting may have been more related to an expansion of the melt area rather than increased melting at the ice edge. Ryan et al. (2019), for example, observed that the snowline on the Greenland Ice Sheet rose by 17 m yr −1 between 2001 and 2012 lowering albedo. Noël et al. (2019) showed that since the early 1990s, the ablation zone expanded by 25% in southern Greenland and 46% in northern Greenland, covarying with changes in cloudiness. As the climate warms, albedo changes by 2100 under an RCP 4.5 scenario is estimated to cause a 113% increase in meltwater storage on the surface of the Greenland Ice Sheet (Ignéczi et al., 2016). These increases in meltwater availability on the ice surface will likely lead to more extreme changes in supraglacial stream sediment coverage throughout the melt season. Additionally, as the melt season extends later into the year, supraglacial floodplains will be exposed for longer periods of time instead of being covered by snow. This would amplify the effect of sediment deposition on melting. Inter-seasonal sediment cover changes will also likely be influenced by changes in sediment supply. Lower surface mass balance values result in greater rates of interstitial dust particle melt-out. Additionally, as the Ice Sheet retreats and exposes more glacial till at the ice margin and wind speeds continue to increase (Zhang et al., 2021), locally derived sediment loads will likely increase as well. It is still unclear if the increase in streamflow will be large enough to accommodate the expected increase in sediment load.

Impacts on Greenlandic carbon cycle
Supraglacial streams are also a major component of the Greenland Ice Sheet supraglacial ecosystem . While carbon fluxes from supraglacial streams due to loss of sediment cover was not investigated in this study, the relatively constant stream hydrological connectivity and decrease in sediment cover suggests that carbon within supraglacial stream sediment is being delivered to the bed via moulins. As such, we can speculate as to how this might impact the carbon cycle of the Ice Sheet. As sediment granules in supraglacial streams flow into moulins, the organic matter delivered to the bed can potentially lead to increased subglacial anaerobic productivity that causes an elevated export of methane from the Ice Sheet (Lamarche-Gagnon et al., 2019;Christiansen and Jorgensen, 2018). Previous studies have shown that several subglacial catchments in Greenland are significantly under-saturated in organic carbon (Pain et al., 2021). As such, as supraglacial streamflow increases and delivers more cryoconite to the base of the Ice Sheet, this may impact methanogenesis and further contribution to atmospheric greenhouse gas concentrations. Lee et al. (2020) found that increased dust and Ca 2+ concentrations in Greenland ice cores were highly correlated with increased methane concentrations potentially due to methanogenesis by extremophiles. Runoff generated at the surface can also lead to substantial subglacial erosion, which releases suspended sediment to the ocean and leads to a non-linear response of summertime marine productivity (Overeem et al., 2017;Hopwood et al., 2018). The extent of this process and the influence of supraglacial streams is still highly understudied and further research is needed to quantify how supraglacial carbon fluxes impact emissions.

Conclusion
Supraglacial streams and the sediment they contain play a crucial role in determining the albedo of the Greenland Ice Sheet. As far as we are aware, this is the first study to show that stream sediment concentrations change over the course of a melt season and how these changes are linked to streamflow conditions. As supraglacial streams discharge decreases after peak melt, floodplains respond by increasing sediment cover and lowering ablation zone albedo. This is, in effect, a negative feedback where the expanded sediment cover at the end of the season lowers albedo and potentially extends the melt season. This in turn can impact ice dynamics and mass loss. Compared to sediment in cryoconite holes, sediments in supraglacial streams are far more prevalent and receive substantially more solar radiation since most of this sediment is deposited in exposed floodplains. As climate changes, sediment deposition will likely respond to increased flow rates and seasonal variability potentially impacting supraglacial and Frontiers in Earth Science frontiersin.org subglacial ecosystems. Further research to model how supraglacial streams and floodplains respond to changing climates and how they impact the absorption of solar radiation on the Greenland Ice Sheet would therefore improve predictions of Greenland's contribution to greenhouse gas emissions and sea level rise.

Data availability statement
The data used for this article are available through the Arctic Data Center, including all DEMs and orthomosaic rasters generated from UAV imagery (Leidman et al., 2023), water stage data and discharge measurements for the supraglacial stream (Leidman et al., 2023), and spectrometry measurements taken with the ASD device (Leidman et al., 2023). 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) and can be found at https://dataverse.geus.dk/dataverse/ AWS. Modéle Atmsphérique Régional (MAR) data were provided by Xavier Fettweis, University of Liege, Belgium.

Author contributions
SL conducted the majority of the field data collection, data analysis, and writing with assistance from ÅR. RM and SS aided in field data collection. SL, ÅR, AG, and SS were instrumental in the conceptualization and editing of the paper.

Funding
Funding for this research was provided by the National Science Foundation Graduate Research Fellowship Program and the NASA Cryospheric Science Program Grants #80NSSC19K0942.