A new surface meltwater routing model for use on the Greenland Ice Sheet surface

Large volumes of surface meltwater are routed through supraglacial internally drained catchments (IDCs) on the Greenland Ice Sheet surface each summer. Because surface routing impacts the timing and discharge of meltwater entering the ice sheet through moulins, accurately modeling moulin hydrographs is crucial for correctly coupling surface energy and mass balance models with subglacial hydrology and ice dynamics. Yet surface routing of meltwater on ice sheets remains a poorly understood physical process. We use high-resolution (0.5 m) satellite imagery and a derivative high-resolution (3.0 m) digital elevation model to partition the runoff-contributing area of the Rio Behar catchment, a moderately sized (∼ 63 km2) midelevation (1207–1381 m) IDC in the southwestern Greenland ablation zone, into open meltwater channels (supraglacial streams and rivers) and interfluves (small upland areas draining to surface channels, also called “hillslopes” in terrestrial geomorphology). A simultaneous in situ moulin discharge hydrograph was previously acquired for this catchment in July 2015. By combining the in situ discharge measurements with remote sensing and classic hydrological theory, we determine mean meltwater routing velocities through open channels and interfluves within the catchment. Two traditional terrestrial hydrology surface routing models, the unit hydrograph and rescaled width function, are applied and also compared with a surface routing and lake-filling model. We conclude that (1) surface meltwater is routed by slow interfluve flow (∼ 10−3–10−4 m s−1) and fast open-channel flow (∼ 10−1 m s−1); (2) the slow interfluve velocities are physically consistent with shallow, unsaturated subsurface porous media flow (∼ 10−4–10−5 m s−1) more than overland sheet flow (∼ 10−2 m s−1); (3) the open-channel velocities yield mean Manning’s roughness coefficient (n) values of ∼ 0.03–0.05 averaged across the Rio Behar supraglacial stream–river network; (4) interfluve and open-channel flow travel distances have mean length scales of ∼ 100–101 m and ∼ 103 m, respectively; and (5) seasonal evolution of supraglacial drainage density will alter these length scales and the proportion of interfluves vs. open channels and thus the magnitude and timing of meltwater discharge received at the outlet moulin. This phenomenon may explain seasonal subglacial water pressure variations measured in a borehole ∼ 20 km away. In general, we conclude that in addition to fast open-channel transport through supraglacial streams and rivers, slow interfluve processes must also be considered in ice sheet surface meltwater routing models. Interfluves are characterized by slow overland and/or shallow subsurface flow, and it appears that shallow unsaturated porous-media Published by Copernicus Publications on behalf of the European Geosciences Union. 3792 K. Yang et al.: A new surface meltwater routing model flow occurs even in the bare-ice ablation zone. Together, both interfluves and open channels combine to modulate the timing and discharge of meltwater reaching IDC outlet moulins, prior to further modification by englacial and subglacial processes.

Abstract.Large volumes of surface meltwater are routed through supraglacial internally drained catchments (IDCs) on the Greenland Ice Sheet surface each summer.Because surface routing impacts the timing and discharge of meltwater entering the ice sheet through moulins, accurately modeling moulin hydrographs is crucial for correctly coupling surface energy and mass balance models with subglacial hydrology and ice dynamics.Yet surface routing of meltwater on ice sheets remains a poorly understood physical process.We use high-resolution (0.5 m) satellite imagery and a derivative high-resolution (3.0 m) digital elevation model to partition the runoff-contributing area of the Rio Behar catchment, a moderately sized ( ∼ 63 km 2 ) midelevation (1207-1381 m) IDC in the southwestern Greenland ablation zone, into open meltwater channels (supraglacial streams and rivers) and interfluves (small upland areas draining to surface channels, also called "hillslopes" in terrestrial geomorphology).A simultaneous in situ moulin discharge hydrograph was previously acquired for this catchment in July 2015.By combining the in situ discharge measurements with remote sensing and classic hydrological theory, we determine mean meltwater routing velocities through open channels and interfluves within the catchment.Two traditional terrestrial hydrology surface routing models, the unit hydrograph and rescaled width function, are applied and also compared with a surface routing and lake-filling model.We conclude that (1) surface meltwater is routed by slow interfluve flow (∼ 10 −3 -10 −4 m s −1 ) and fast open-channel flow (∼ 10 −1 m s −1 ); (2) the slow interfluve velocities are physically consistent with shallow, unsaturated subsurface porous media flow (∼ 10 −4 -10 −5 m s −1 ) more than overland sheet flow (∼ 10 −2 m s −1 ); ( 3) the open-channel velocities yield mean Manning's roughness coefficient (n) values of ∼ 0.03-0.05averaged across the Rio Behar supraglacial stream-river network; ( 4) interfluve and open-channel flow travel distances have mean length scales of ∼ 10 0 -10 1 m and ∼ 10 3 m, respectively; and (5) seasonal evolution of supraglacial drainage density will alter these length scales and the proportion of interfluves vs. open channels and thus the magnitude and timing of meltwater discharge received at the outlet moulin.This phenomenon may explain seasonal subglacial water pressure variations measured in a borehole ∼ 20 km away.In general, we conclude that in addition to fast open-channel transport through supraglacial streams and rivers, slow interfluve processes must also be considered in ice sheet surface meltwater routing models.Interfluves are characterized by slow overland and/or shallow subsurface flow, and it appears that shallow unsaturated porous-media 1 Introduction Supraglacial internally drained catchments (IDCs) are hydrologic units on the Greenland Ice Sheet (GrIS) surface that collect and drain surface meltwater through supraglacial stream-river networks to terminal moulins or lakes (Thomsen et al., 1989;Yang and Smith, 2016).IDC spatial and temporal characteristics and processes constrain the location, discharge, lag times, and total volume of surface meltwater penetrating into the ice sheet via moulins (Banwell et al., 2013;Yang and Smith, 2016;Smith et al., 2017), which in turn influences the timing of surface mass loss, subglacial hydrologic system evolution, and ice flow dynamics on seasonal and shorter-term timescales (Zwally et al., 2002;Sole et al., 2011;Banwell et al., 2013Banwell et al., , 2016;;Andrews et al., 2014;Arnold et al., 2014;Clason et al., 2015;Smith et al., 2015).
Previous studies have shown that planform IDC locations and shapes are largely induced by underlying bedrock controls on ice surface morphology, which is influenced by variations in bed roughness and slipperiness and the differing transmission of that variability to the ice surface (Lampkin and VanderBerg, 2011;Karlstrom and Yang, 2016;Crozier et al., 2018;Ignéczi et al., 2018).IDC areas generally increase with elevation due to lower moulin densities at high elevations (Poinar et al., 2015;Yang and Smith, 2016).As a result, high-elevation IDCs can drain nontrivial volumes of meltwater into the ice sheet even where overall melt rates are low (Yang and Smith, 2016;Smith et al., 2017).Analysis of satellite imagery suggests that by mid-July, over 95 % of IDCs in the southwest GrIS drain meltwater into moulins rather than lakes (Fitzpatrick et al., 2014;Smith et al., 2015).The variable sizes and shapes of IDCs influence the timing and magnitude of peak meltwater injection into the ice sheet through these moulins (Smith et al., 2017).
Few studies have examined meltwater surface routing processes through IDCs on the GrIS ablation zone surface.Therefore, our ability to simulate IDC moulin hydrographs and, by extension, realistic surface-to-bed connections remains limited.Most studies to date have used surface melt rates simulated by regional climate models (RCMs) to calculate runoff (hereafter called "moulin discharge") assumed to flow to the IDC's terminal outlet moulin, but without explicit treatment of surface meltwater routing processes and associated time lags (Smith et al., 2017).This can introduce large uncertainty in the timing and magnitudes of meltwater injection into the ice sheet, impacting the accuracy of subglacial hydrology and dynamical ice flow simulations.
Previous studies have addressed this problem in different ways.Clason et al. (2015) used a single-flow direction algorithm to route surface meltwater across the ice surface and accounted for factoring in runoff delays due to snowpack retention; Arnold et al. (1998) developed a distributed surface routing and lake-filling (SRLF) model to simulate moulin discharge.The SRLF model was designed for snowor bare-ice-covered IDCs and has since been used to simulate the effects of up-glacier snow line retreat on meltwater routing (Willis et al., 2002), drive subglacial hydrologic system evolution (Banwell et al., 2013(Banwell et al., , 2016)), fill supraglacial lakes (Banwell et al., 2012;Arnold et al., 2014), and illustrate ice flow patterns (de Fleurian et al., 2016).Leeson et al. (2012) developed a similar surface meltwater routing model based on Manning's equation for open-channel flow and Darcy's law for subsurface flow through a porous medium to simulate lake evolution on the GrIS.
As with terrestrial hydrologic models, supraglacial meltwater routing models are highly sensitive to the choices of surface routing scheme, some poorly quantified parameters (e.g., Manning's roughness coefficient n, mean flow velocity, lag time to peak, near-surface permeability), and input data (e.g., RCM model output, Automated Weather Stations).Karlstrom et al. (2014) measured and modeled subsurface porous flow in weathered ice to infer near-surface permeability and found it to be considerably smaller than commonly assumed parameterizations.Gleason et al. (2016) measured hydraulic geometries (width, depth, velocity, slope, Manning's roughness coefficient n, etc.) of nine supraglacial meltwater channels on the southwest GrIS and found n, in particular, to be variable (0.009-0.154).Smith et al. (2017) measured a 72 h in situ terminal outlet moulin hydrograph for the Rio Behar catchment, a moderately large instrumented IDC (∼ 63 km 2 ).These data were used to empirically calibrate a simple surface meltwater routing model, the Snyder synthetic unit hydrograph (SUH) for broader application across 799 surrounding IDCs remotely sensed across the GrIS ablation zone (Yang and Smith, 2016).However, SUH is a highly "lumped" routing model, meaning it does not distinguish among different physical flow process paths within a catchment (Singh et al., 2014) and thereby cannot differentiate among different meltwater routing processes.Some other more complex SUH methods have also been proposed for terrestrial hydrology but most of those methods cannot partition physical flow processes either (Singh et al., 2014).For example, the geomorphic instantaneous unit hydrograph (GIUH) method only focuses on open-channel flow but ignores hillslope flow (Moussa, 2008), which is not suitable for representing meltwater routing on the ice surface.
This study presents a spatially lumped, process-partitioned meltwater routing model to investigate surface meltwater routing parameters (meltwater travel distance, velocity, and time) in the Rio Behar catchment, using the in situ moulin hydrograph of Smith et al. (2017) for calibration.The lumped spatial domain is the moderate IDC scale (∼ 63 km 2 ).Two traditional terrestrial hydrograph analysis tools, the unit hydrograph (UH) and rescaled width function (RWF), are used to characterize IDC meltwater routing at a high spatial resolution (3 m) afforded by using a high-resolution remotely sensed digital elevation model (DEM) acquired simultaneously with the in situ measurements by the WorldView-1 satellite.The RWF offers particular advantages over SUH because it can partition between two different types of flow path (interfluve vs. open channel) and quantify their respective mean (or "bulk") meltwater travel velocities within the catchment.The performance of our RWF model is also compared with SRLF to better estimate the performance of RWF and to assess two different routing approaches.The resultant travel distance, velocity, and time are used to characterize two physically different meltwater routing processes on the ice surface, and the implications of seasonal routing evolutions are discussed.We conclude that a calibrated, highresolution RWF surface routing model offers good utility for obtaining important meltwater routing parameters and modeling IDC outlet moulin discharge hydrographs in the GrIS bare-ice ablation zone.

Study area
Our study area is the Rio Behar catchment (Smith et al., 2017), a moderately sized supraglacial IDC on the southwest GrIS surface (Fig. 1).It is located in the upper ablation zone, spanning 1207-1381 m a.s.l., near the long-term equilibrium line altitude (∼ 1550 m a.s.l.) of this area (van de Wal et al., 2015).In July 2015, the remotely sensed area of the Rio Behar catchment was 63.1 km 2 and the main stem length of its trunk supraglacial river was 13.8 km (Smith et al., 2017).Visual analysis of multi-temporal high-resolution satellite and UAV (drone) images shows that the supraglacial stream-river channel network was highly developed in the Rio Behar catchment by 21-23 July 2015, when 93.5 % of the surface was bare ice.A detailed description of the Rio Behar catchment and remotely sensed imagery, DEM, and catchment map is provided in Smith et al. (2017).Notably, the IDC exhibits a sub-grid scale with respect to RCMs and ice sheet models (ISMs) and is considered the dominating surface meltwater routing process on the southwest GrIS surface (Yang and Smith, 2016).

Data sources
For 72 continuous hours (11:00, 21 July to 10:00 western Greenland summer time (UTC − 2), 23 July 2015) measurements of meltwater discharge exiting the Rio Behar catchment were collected using an acoustic Doppler current profiler (ADCP) in the main stem supraglacial river (Smith et al., 2017).Hourly simulations of GrIS meltwater production M over the study period were generated using the MAR (Modèle Atmosphérique Régionale) 3.6 RCM (Fettweis et al., 2013).The 20 km MAR grid cells were reprojected to a common 5 km posting and map projection using nearest-neighbor resampling.Catchment-mean hourly melt (mm h −1 ) was obtained by clipping MAR grid cells with the remotely sensed Rio Behar boundary (Smith et al., 2017) and summing their corresponding melt values (Fig. 2) and was multiplied by the measured runoff coefficient for the catchment (0.69) to yield units of effective melt M (Smith et al., 2017).
Two catalogs of stereo WorldView-1 (WV1) panchromatic images (spatial resolution 0.5 m) acquired on 18 July 2015 were used for detailed mapping of supraglacial hydrologic features (rivers, lakes, and moulins) (Smith et al., 2017).A fixed-wing UAV (Ryan et al., 2015) acquiring aerial camera imagery (RGB bands) over the Rio Behar catchment from 20 to 22 July 2015 was used to validate the accuracy of supraglacial stream-river delineations derived from the WV1 image.
A concurrent high-resolution (3 m) DEM was derived from these WV1 stereo images using the open-source Ames Stereo Pipeline (ASP) toolkit method (Shean et al., 2016).We used the 30 m GIMP (Greenland Ice Mapping Project) DEM v2 (Howat et al., 2014) for comparison with the WV1 DEM.

Remote sensing of Rio Behar supraglacial river network
A supraglacial stream-river network integrates the hydrologic response of an IDC to surface melt from outside the supraglacial channel system (e.g., the "hillslope" in terrestrial hydrology; D' Odorico and Rigon, 2003) and openchannel flow (Montgomery and Foufoula-Georgiou, 1993).By late July the ice surface of the GrIS ablation zone becomes heavily dissected with very high drainage density (Smith et al., 2015;Yang et al., 2017) and very short distances between open channels.The term "interfluve" is borrowed from terrestrial fluvial geomorphology and refers to small areas of dissected terrain that slope toward rills or gullies.Interfluves are commonly referred to using the more general term hillslope; however in terrestrial geomorphology hillslopes can also refer to much larger features, whereas interfluves are a narrower term used for small upland areas with short runoff distances, typically found on heavily dissected surfaces.It is thus the more appropriate term for use on the ice sheet owing to the high observed stream density and correspondingly short distances between supraglacial open channels.Moreover, on the ice surface, there is no analog of soil creep, which is an important process in terrestrial hillslope geomorphology (Montgomery andDietrich, 1989, 1992).For these reasons we recommend use of the narrower term interfluve instead of the general term of hillslope, although the mathematical treatments are similar for both.
www.the-cryosphere.net/12/3791/2018/The Cryosphere, 12, 3791-3811, 2018 We delineated actively flowing supraglacial streams-rivers from the 0.5 m panchromatic WV1 image, following the automatic detection method of Yang et al. (2017) and Smith et al. (2017).The concurrent UAV image was used to validate the ability of this 0.5 m WV1 image to capture small streams.Two detection thresholds, one conservative (higher) and one nonconservative (lower) with values of 40 and 5 (out of 255), respectively, were applied separately to create two meltwater masks following Gabor filtering and path-opening processing of the WV1 image (Yang et al., 2017).The conservative threshold extracted linear features that are confidently classified as open-flow channels with clear channel banks and with high spectral contrast from the surrounding ice (Smith et al., 2017), while the nonconservative threshold extracted all the channel-like features in the image (Yang et al., 2015a).The two resultant river masks therefore represent upper and lower limits for the true distribution of the openchannel supraglacial stream-river network that was actively flowing on the Rio Behar catchment when the in situ hydrograph was collected.

High-resolution DEM processing
Extracting an IDC supraglacial stream-river network from a DEM requires assignment of a prescribed location for the catchment outlet (sink).For this study, the topographic depression containing the known location of the terminal outlet moulin was used as the sink; all other small depressions were filled as per Karlstrom and Yang (2016).This partially filled DEM was then used to calculate flow directions and a downstream flow-contributing area raster (Karlstrom and Yang, 2016).Finally, a global meltwater contribution area (A c ) threshold was used to simulate ice surface drainage networks.In practice, if A c is set too large (small), modeled drainage networks will underestimate (overestimate) realworld channel travel distances and overestimate (underestimate) actual interfluve travel distances (Montgomery and Foufoula-Georgiou, 1993;Yang and Smith, 2016).There-Figure 2. Within the Rio Behar catchment, MAR regional climate model inputs of melt (at top) are paired with in situ acoustic Doppler current profiler (ADCP) outputs of supraglacial river discharge (at bottom) to calibrate the unit hydrograph (UH) surface routing models assessed in this study."Effective melt" is the fraction of total MAR melt production that is transported all the way to the Rio Behar IDC terminal outlet moulin (Smith et al., 2017).The UH is calculated using effective melt as input and observed runoff as output as per Smith et al. (2017).
fore, by deliberately varying this parameter we are able to simulate the seasonal evolution of the supraglacial streamriver network, which tends to have the lowest drainage density early and late in the melt season (Yang et al., 2015b(Yang et al., , 2017;;King et al., 2016).This study used a DEM stream burning technique to force the DEM to produce a reliable actively flowing river network (Lindsay, 2016).To burn the WV DEM, elevations of DEM raster pixels that are spatially coincident with our remotely sensed supraglacial map were lowered ("burned") by 1.0 m, thereby forcing routed flow to pass through these accurately mapped supraglacial streamriver channels.

Unit hydrograph
The UH is a transfer function used to simulate the observed hydrograph (Q) at a catchment outlet for a unit input of water supply (e.g., 1 mm, 1 cm) applied uniformly across the catchment (Dingman, 2015).Effective M is the fraction of total MAR melt production M that is transported all the way to the terminal outlet moulin, i.e., the remainder after multiplying by the field measured runoff coefficient (0.69) (Smith et al., 2017).M is used as input to the transfer function, and the resulting simulation is called the "direct hydrograph" (as distinguished from the observed hydrograph).The UH was developed for terrestrial hydrologic applications in which precipitation (rain or snow) is the dominant hydrologic input (Singh et al., 2014) but is well suited for adaptation on ice sheets by substituting measured or modeled water equivalent from ice melt (Smith et al., 2017), as ice melt is the dominant hydrological input to IDCs on the southwest GrIS.Rainfall occasionally occurs during summer in our study area (Van As et al., 2017) but none occurred during our study period.
This allows MAR effective melt (M ) to be used as the input data for simulating (routing) the direct hydrograph at the IDC terminal outlet moulin as Q = M * UH, where * is a convolution operator (Dingman, 2015).Smith et al. (2017) used the Collins' method (Collins, 1939) to create a UH specific to the Rio Behar catchment (here called "MAR UH"), using MAR-produced M and ADCP-measured Q to calibrate the UH.
In the present study, we also derive other UH transfer functions, which can be built from a satellite image or DEM.To do this, the travel time (t) for each pixel within the IDC is required.Travel time represents the time needed for water to flow from each pixel to the catchment outlet, and the hourlybinned histogram of this travel time raster thus corresponds to a 1 h UH (Liu et al., 2003).Travel time t can be estimated as t = L/v, where L is meltwater flow distance and v is flow velocity.Flow distance L can be calculated from DEMs by assuming that meltwater flows from one pixel to the adjacent pixel with the steepest slope (Karlstrom and Yang, 2016).

Surface routing and lake-filling (SRLF) model
The SRLF model is a distributed, physically based model proposed by Arnold et al. (1998).Model input requirements include DEM elevations and a time series of meteorological forcing data.The model has been widely used for studies of surface meltwater routing in the GrIS ablation zone (Willis et al., 2002;Banwell et al., 2012Banwell et al., , 2013Banwell et al., , 2016;;Arnold et al., 2014;de Fleurian et al., 2016).The SRLF model uses Darcy's law to route meltwater flow through snow and Manning's equation to route meltwater flow over bare-ice surfaces.The present study focuses exclusively on the latter because satellite and UAV mapping revealed the surface of the Rio Behar catchment to be virtually all bare ice during the observational period (Smith et al., 2017).
The SRLF model can be used to calculate meltwater travel time and thus to create UH.In accordance with Arnold et al. (1998), we calculated a meltwater flow velocity (v) for each pixel in the Rio Behar catchment (see Appendix A).To compute meltwater travel time, a cost surface was created as 1/v, which was then used as one input raster and flow direction raster was used as another input to determine specific flow paths.The output raster gives travel time for each pixel.The SRLF UH was created by hourly binning the histogram of this travel time raster.

Rescaled width function (RWF)
The RWF is a conceptual runoff routing model that represents the total flow distance (L) in a catchment as a combination of an interfluve (hillslope in terrestrial settings), flow www.the-cryosphere.net/12/3791/2018/The Cryosphere, 12, 3791-3811, 2018 distance (L h ), and a channel flow distance (L c ), i.e., L = L h + L c .The RWF is an improved version of the width function (WF).The WF does not represent interfluve transport and therefore cannot be used for meltwater routing where interfluve flow is important (see Appendix B).If constant flow velocities are assumed for interfluve (v h ) and channel (v c ) zones, the travel time (t) for each pixel is the sum of interfluve travel time (t h ) and channel travel time (t c ) (Di Lazzaro, 2009).Consequently, catchment UH can be derived from the travel time, which is renamed the RWFUH (Singh et al., 2014).Determination of v h and v c is challenging because water flow velocities are difficult to measure on a catchment scale, especially for interfluve zones (Moussa, 2008).Although water flow velocities can be predicted theoretically for simple porous flow (Karlstrom et al., 2014) and our more general derivation in Appendix C, the field-measured IDC moulin hydrograph provides a good opportunity to directly calibrate v h and v c for a real-world melting ice sheet surface.To achieve this, different combinations of v h and v c were used to create RWFUHs.These RWFUHs were then used to simulate direct hydrographs at the IDC terminal outlet moulin for comparison with the corresponding in situ moulin hydrograph.To evaluate performance between the simulated vs. observed moulin hydrograph, we used the Nash-Sutcliffe model efficiency (NSE) cost function (Nash and Sutcliffe, 1970).The NSE was calculated for each RWFUH and compared to the field-measured moulin hydrograph and the "optimal" v h and v c defined as the combination that maximizes the NSE.

Seasonal evolution of the supraglacial stream-river network
By late July, shortwave radiation and air temperatures decline and meltwater production within the Rio Behar catchment also declines (Smith et al., 2017).In this study, we assume that low-order streams (i.e., very small tributaries as per Yang et al., 2016) stop flowing by late July (Yang et al., 2017), higher-order streams-rivers stop flowing by mid-August, and only very large high-order rivers flow into late August, an assumption supported by remotely sensed supraglacial streamriver maps in Smith et al. (2015) and Yang and Smith (2016).
In late July when the supraglacial drainage density is high, fast open-channel transport contributes heavily to meltwater flow routing exiting the IDC.As the supraglacial river network declines and drainage density decreases, interfluve zones expand and contribute more surface area to overland and/or porous media transport processes.
To test this idea, we simulated a temporal evolution of the supraglacial stream-river network within the Rio Behar IDC after late July to characterize the impact of drainage density decline on the RWFUH hydrograph.This test consisted of defining a series of A c thresholds (i.e., 250, 500, 1000, 2500, and 5000 pixels) that were used to create ar-tificial supraglacial drainage networks from the WorldView DEM.A c indicates the minimum meltwater contribution area required to form a supraglacial headwater stream.If a DEM grid cell exhibits a contributing area larger than A c , a supraglacial stream will form and thereby the grid cell belongs to the open-channel zone.In contrast, if a DEM grid cell exhibits a contributing area smaller than A c , a supraglacial stream will not form and thereby the grid cell belongs to the hillslope zone.Larger A c values will yield larger hillslope zones, whereas smaller A c values will yield larger open-channel zones.The minimum A c (250 pixels) was used to simulate a well-developed supraglacial stream-river network, while the maximum A c (5000 pixels) was used to simulate a poorly developed stream-river network.Variable A c values were used to simulate dynamic supraglacial streamriver networks and each A c value was assumed to last for 1 week.This sequence of A c thresholds reasonably mimics the seasonal contraction of the supraglacial stream-river networks after their maximum development in late July (Smith et al., 2015).The resultant drainage networks were then used to calculate moulin discharge (i.e., the direct hydrograph) based on optimal open-channel and hillslope velocities calibrated from RWFUH to demonstrate the influence of seasonal drainage network evolution on the shape and timing of the discharge hydrograph at the IDC terminal outlet moulin.

Supraglacial stream-river mapping
In total, 3381 km of actively flowing supraglacial streamriver lengths was confidently mapped (conservative threshold) within the Rio Behar catchment, yielding a drainage density of 53.6 km −1 with water bodies covering 8.2 % of the catchment surface.Applying the nonconservative detection threshold, 10 829 km of supraglacial stream-river lengths was mapped, yielding a higher drainage density of 164.3 km −1 with water bodies covering 24.1 % of the ice surface.These bounding estimates suggest that 76 %-92 % of the ice surface area consisted of interfluve zones with the remainder as supraglacial ponds, lakes, rivers, and streams.
Four test sites were selected to further illustrate the performance of conservative and nonconservative stream-river channel detections (Fig. 3).Sites 1 through 3 are located near the three main tributaries of the Rio Behar catchment where supraglacial streams-rivers are very well developed.Site 4 is located upstream, within an area where supraglacial streams are relatively sparse.Applying the conservative detection threshold delineated supraglacial streams-rivers clearly identifiable in the WV1 image, but narrower channellike (dark linear but not well-channelized) features among those streams-rivers were missed.In contrast, the nonconservative threshold captured these small features with the resultant streams-rivers being very well developed and hav- ing higher drainage density.Visual inspection of the 0.3 m UAV images (RGB bands) reveals that supraglacial channels mapped with 0.5 m resolution WV satellite imagery capture nearly all conservative channels that can be discerned in 0.3 m UAV camera imagery (Fig. 3).However, numerous smaller nonconservative supraglacial streams among large supraglacial rivers were not delineated.Therefore, we conclude that automatically mapped streams-rivers accurately estimate the minimum and maximum extents of channel and interfluve zones in the Rio Behar catchment.

Interfluve and open-channel travel distances
Meltwater travel distance rasters were calculated for each data source and processing approach (Table 1).Figure 4 shows the interfluve and channel travel distances with conservatively mapped rivers and burned WV DEM.The resultant mean channel travel distance is 7.1 ± 4.0 × 10 3 m, while the resultant mean interfluve travel distance is 19.7 ± 30.9 m.This signifies that in the Rio Behar catchment during our study period, meltwater travel distances through open channels were ∼ 3 orders of magnitude longer than travel distances through interfluves.
The mean interfluve distance we estimate for the Rio Behar catchment is larger than the 0.5-5 m values reported for the Juneau Icefield (Karlstrom et al., 2014) and the 9.0 ± 3.4 m interfluve distance reported for another low-elevation IDC of the southwest GrIS (McGrath et al., 2011).However, their interfluve distances were calculated as the nearest distance from a interfluve point to its adjacent channel rather than following the topographic flow direction calculated from DEMs, as in Karlstrom et al. (2014) and McGrath et al. (2011) (although for small slopes the two approaches should yield similar results).River detection thresholds significantly impact interfluve distances because higher-density Burned WV DEM, conservative threshold 7.1 ± 4.0 19.7 ± 30.9 0.4 0.9 -0.9443 4.9 ± 2.8 6.1 ± 9.5 0.4 0.7-1.2-> 0.9250 0.3-0.distributed meltwater channels lead to smaller interfluve distances.If the nonconservative river detection threshold is used, the mean interfluve distance calculated from our burned WV DEM is 6.7±15.0m, closer to these previously reported estimates.

Interfluve and open-channel travel velocities
The optimal RWF-calibrated mean open-channel velocity v c is on the order of 10 −1 m s −1 , while the optimal mean velocity v h for interfluves is on the order of 10 −3 -10 −4 m s −1 (Table 1 and Fig. 5).Both conservative and nonconservative approaches quantify open-channel velocities as v c = 0.3-0.5 m s −1 , which is similar to previous field-measured values of 0.25-0.5 m s −1 in small supraglacial streams (width < 0.5 m) at the Juneau Icefield (Karlstrom et al., 2014) and 0.35 m s −1 measured for a small supraglacial stream (0.2 m wide) at the southwest GrIS (Gleason et al., 2016).It is slower than faster velocities (0.2-9.4 m s −1 ) measured in large (> 10 m) supraglacial rivers (Smith et al., 2015;Gleason et al., 2016).The relatively low v c = 0.3-0.5 m s −1 quantified for an entire catchment suggests that small, relatively slow-flowing supraglacial streams that are vastly more numerous than large main stem supraglacial rivers dominated the mean RWF open-channel velocity, which is a bulk velocity averaged over the entire IDC (Appendix D).
The RWF-calibrated optimal interfluve velocity v h shows larger variation than v c .The range of v h is interpreted as v h = 0.2-1.5 × 10 −3 m s −1 if NSE = 0.9.If NSE = 0.925 is used as the calibration threshold, this optimal v h range is 0.3-1.2× 10 −3 m s −1 .Under the assumption that conservative (nonconservative) rivers overestimate (underestimate) interfluve distance, the nonconservative v h should be considered as the lower v h limit, while the conservative v h as the upper v h limit (Table 1).Although the specific v h may vary according to different calibration thresholds, Fig. 5 suggests that v h is on the order of 10 −3 -10 −4 m s −1 .This finding confirms that meltwater is routed through the Rio Behar catchment by slow interfluve flow (∼ 10 −3 -10 −4 m s −1 ) followed by fast open-channel flow (∼ 10 −1 m s −1 ).

Interfluve and open-channel travel time
The total supraglacial travel time (i.e., the combination of interfluve and channel travel time) for our study area and period was found to be ∼ 11 h.This suggests that, on average, a unit of application meltwater across the catchment takes 11 h to arrive at the terminal outlet moulin.However, this should be distinguished from the "time to peak", which describes the lag time between peak meltwater production and peak runoff entering the IDC terminal outlet moulin (Smith et al., 2017), using real-world inputs of melt production, which follow a strongly diurnal cycle.The 11 h reported here refers to a "unit" response, i.e., the average length of time for an instantaneous pulse of 1 cm of meltwater to drain from the Rio Behar IDC.
The optimal v h and v c combinations (which maximize NSE) were used to calculate mean meltwater travel time in interfluve (t h ) and open-channel (t c ) zones.The resultant mean t h is ∼ 6 h, while the mean t c is ∼ 5 h.This result differs from the results obtained for smaller bare-ice catchments (< 2 km 2 ), in which interfluve travel primarily controls or even dominates meltwater routing (Arnold et al., 1998;Karlstrom et al., 2014).For small catchments, meltwater channels are short and therefore fast travel time through open channels is less important than slow travel time through interfluve zones.

Moulin hydrograph simulations
UHs were created from the meltwater travel time maps (driven by the optimal v h and v c combinations in Table 1), allowing direct hydrographs at the Rio Behar IDC terminal outlet moulin to be simulated.Similarly, the SRLF model was applied to the WV DEM and the GIMP DEM to create UHs as well (Fig. 6), allowing direct comparisons with our RWF-based methods.Two contributing area thresholds (A c = 450 and 90 m 2 ) were applied to model supraglacial drainage networks with large (∼ 23 m) and small (∼ 9 m) interfluve distances, which were used for comparisons with the conservative and nonconservative image-mapped river networks, respectively (Table 1).
The UHs simulated by four RWF-based approaches, using two burned WV DEMs and two A c -based WV DEMs, generally capture the overall shape of the MAR UH (with a duration of 25 h; see Smith et al. (2017) for more details).All four RWF-based UHs smooth the MAR UH, signifying that surface routing of meltwater distributes runoff more uniformly over time (Dingman, 2015).The SRLF-based WV DEM UH also performs reasonably well, although its shape is more irregular, similar to the MAR UH.The UH simulated by the 30 m GIMP DEM is different from all the other UHs in that it distributes all of the input meltwater during the first 13 h, suggesting that this coarse-resolution DEM overly speeds up the surface meltwater transport (Fig. 6a).
Both RWF-based and SRLF-based UHs simulate the moulin hydrograph well (Fig. 6b).Except for the GIMP DEM SRLF approach, all the other RWF-based and SRLFbased approaches were able to accurately simulate the peak discharges of the observed moulin hydrographs.These approaches also captured the peak time of the first daily hywww.the-cryosphere.net/12/3791/2018/The Cryosphere, 12, 3791-3811, 2018 drograph, whereas a 2-4 h time lag was obtained for the second daily hydrograph.However, because SRLF lacks slow interfluve flow, it routes surface meltwater faster than RWF and distributes all of the input meltwater during the first 20 h (Fig. 6a).Ignoring the slow interfluve flow affects the performance of the SRLF method negatively, yielding an optimal NSE = 0.8742, smaller than for the RWF method (Table 1).Moreover, the SRLF-based GIMP DEM hydrograph is considerably different from the observed moulin hydrograph.This suggests that DEMs with a resolution exceeding ∼ 30 m may fail to accurately capture the velocities and time delays of surface meltwater routing.

Performance of conventional DEM-based simulations
The two conventional DEM-based simulations (Sect.4.2), assuming large and small values of the threshold A c , yielded results similar to those of the burned DEM approaches.The supraglacial drainage network simulated by a relatively large A c (450 m 2 , equivalent to 50 WV DEM pixels) is similar to the conservative image-mapped streams-rivers, whereas the drainage network simulated by the small A c (90 m 2 , equivalent to 10 WV DEM pixels) is similar to the nonconservative image-mapped streams-rivers (Table 1).Depending on this choice of A c threshold and also NSE threshold we find optimal v h values ranging from 0.3 to 1.2 × 10 −3 m s −1 , which is consistent with the optimal range obtained by using the burned WV DEM (v h = 0.2-1.5 × 10 −3 m s −1 ).However, v c shows large variations in optimal values, ranging from 0.4 to 2.0 m s −1 , because DEM-modeled supraglacial drainage networks do not match very well with remotely sensed river networks (Yang et al., 2015b), especially for small supraglacial streams (King et al., 2016).Consequently, the lower value of A c = 90 m 2 is recommended for use during the peak melting period if a high-resolution remotely sensed supraglacial stream-river map is not available.

Seasonal evolution of moulin discharge hydrographs
Seasonal changes in the relative proportion of open-channel vs. interfluve zones substantially alter the timing and magnitude of moulin discharge hydrographs (Fig. 7).If the supraglacial stream-river network is well developed (i.e., has high drainage density) and the interfluve zone is small, large diurnal variations in moulin discharge are simulated.This finding suggests that under well-developed conditions, openchannel travel is particularly important, similar to results reported for the SRLF model (Banwell et al., 2013).As the melt season progresses, smaller supraglacial streams dry up and their associated open-channel zone shrinks.Consequently, open-channel travel becomes secondary to interfluve travel.Under these conditions, meltwater delivery to the englacial system is further attenuated, yielding smaller diurnal variations (Fig. 7).This suggests that in absence of a well-developed supraglacial stream-river network, slow interfluve meltwater transport has a "smoothing" effect on terminal outlet moulin discharge.Similar behavior has been observed or simulated (Arnold et al., 1998;Karlstrom et al., 2014) and is explicitly modeled here using RWF.

Surface runoff delays on the Greenland Ice Sheet
The meltwater travel times quantified in this study confirm nontrivial-runoff delays are caused by at least two fluvial meltwater transport processes operating within the Rio Behar catchment.Such delays have previously been considered insignificant in studies of small ice surface catchments (Karlstrom et al., 2014).However, even for the moderately sized (∼ 63 km 2 ) Rio Behar catchment, supraglacial rivers are long (> 10 km long main stem), meaning meltwater can take several hours to pass through the open-channel network.In much larger (e.g., ∼ 245 km 2 reported in Yang and Smith, 2016) IDCs, channel routing delays are even longer.Therefore, the present study reinforces the importance of supraglacial stream-river networks in imparting nontrivial delays on surface meltwater transport as a function of IDC area, shape, and stream length (Smith et al., 2017), with a new contribution of considering slow interfluve flow as well as fast open-channel flow.
In contrast to these prior studies, our results suggest that both interfluve and open-channel processes control the timing and magnitude of meltwater transport on the ice sheet.This finding suggests that slow meltwater passage over short distances on interfluves is compensated for by fast meltwater transport over long distances through open channels, such that interfluve travel time was roughly equal to channel travel time during the time of the 2015 field experiment.It is possible, therefore, that supraglacial stream-river networks may mimic the classic graded river concept (Kesseli, 1941;Mackin, 1948), with the open-channel flow network developing into a sufficient density to convey available water supply generated on bare-ice interfluves.
Left untreated, surface routing delays degrade the utility of using RCM models to estimate inputs of meltwater to the subglacial environment and proglacial zone.Most current RCM models do not provide any surface routing functions to represent transport of runoff over the ice surface to moulins ( Van As et al., 2014;Cullather et al., 2016).To the best of our knowledge, MAR is the only RCM model integrating a runoff delay function to distribute runoff over time.This delay function was proposed by Zuo and Oerlemans (1996) and is based on the idea that surface meltwater reaches supraglacial channels sooner where the general ice surface slope is larger.The coefficients in the delay function were calibrated by albedo observations on the ice surface, and Lefebre et al. (2003) updated the coefficients to route meltwater more quickly.Applying the MAR delay function, the resultant runoff delay for the Rio Behar catchment is 8.6 days based on Zuo and Oerlemans (1996) and 7.5 days based on Lefebre et al. (2003).In contrast, the runoff delay obtained using RWF in our study is only ∼ 11 h, much shorter than these lumped delays.Van As et al. (2017) built a statistically based supraglacialto-proglacial delay function to optimally match modeled runoff with observed proglacial river discharge measurements collected in the Watson River, near Kangerlussuaq.Applying this delay function, the runoff delay for the entire glacial and proglacial system is 3.9 days for meltwater generated in the Rio Behar catchment.Van As et al. (2017) used a 10 h smoothing per 100 m elevation bin to represent supraglacial routing delays, comparable to our 11 h travel delay calculated for the Rio Behar catchment.Because this purely statistical approach of Van As et al. ( 2017) is calibrated using time series of in situ proglacial discharge measurements rather than formulas applied to DEMs as per Zuo and Oerlemans (1996) and Lefebre et al. (2003), its close agreement with our RWF routing model lends confidence in its more physically realistic routing scheme.Note that the catchment areas of most IDCs are commonly smaller than one MAR cell (Yang and Smith, 2016;Smith et al., 2017).Therefore, for all but the largest IDCs many of the surface routing delays www.the-cryosphere.net/12/3791/2018/The Cryosphere, 12, 3791-3811, 2018 modeled explicitly here could plausibly be parameterized at the scale of a single large RCM grid cell.

Seasonal evolution of the supraglacial drainage network
Supraglacial stream-river networks undergo a dramatic seasonal evolution from low to high to low drainage density within just 3-4 months (Lampkin and VanderBerg, 2014;Smith et al., 2015), modifying the shape of IDC terminal outlet moulin hydrographs.Figure 7 suggests that the moulin hydrograph of the Rio Behar catchment will show small diurnal variations at the beginning and end of a melt season and large diurnal variations during peak melt season.Note that this differs from the classic signal of alpine glaciers, which tend to display a steadily intensifying diurnal cycle throughout the summer as seasonal snow lines climb to higher elevations (Elliston, 1973).Because this seasonal variation may significantly modulate subglacial water pressure and consequently ice flow velocities, this effect warrants further study for other IDCs and using ice dynamical models.Very interestingly, moulin discharges simulated by Fig. 7 are qualitatively similar to subglacial water pressure variation measured in a borehole ∼ 20 km away from our catchment (lat 67.201, long −49.289;Fig. 7 in Wright et al., 2016).Although these field-measured subglacial water pressures were obtained during 2011, they show similarly large diurnal variations during late July, smaller diurnal variations during early August, and very small diurnal variations around late August.This may indicate a direct control of seasonally varying surface meltwater routing on subglacial water pressure, which in turn impacts subglacial pathway evolution and ice flow dynamics (Banwell et al., 2013;Wright et al., 2016).Meanwhile, the subglacial drainage network is well developed in late summer (August) (Andrews et al., 2014) so this may also contribute to the small diurnal variation in subglacial water pressure.
The SRLF model is the most commonly used model for simulating meltwater delivery to moulins and our results suggest it performs well for surface hydrologic conditions similar to those in our study IDC (Banwell et al., 2013(Banwell et al., , 2016)).More recently, the SUH has been advanced as a simple way to model the magnitude and timing of moulin runoff based on remotely sensed IDC area, shape, and main-stem stream length (Smith et al., 2017).However, neither SRLF nor SUH considers the seasonal evolution of supraglacial stream-river networks.Owing to its flexibility for partitioning interfluve and open-channel zones and their changing ratios over time (Mutzner et al., 2016), the RWF provides a good opportunity for improved representation of both interfluve and openchannel processes and their evolution over space and time.
The RWF model is calibrated here using a single in situ supraglacial river discharge hydrograph, so the obtained optimal meltwater routing velocities can thus only be confidently attributed to one particular IDC (Rio Behar catch-ment) during one particular study period (21-23 July 2015).However, other bare-ice IDCs surrounding the Rio Behar catchment have similar surface conditions during bare-ice conditions, suggesting some transferability of our results.For example, the calibrated open-channel velocity, i.e., v c = 0.3-0.5 m s −1 , can be used to estimate Manning's n (n = R 2/3 S 1/2 /v c ); if hydraulic radius R is set to 0.035 m (Arnold et al., 1998) and slope S is set to the mean ice surface slope of the Rio Behar catchment, i.e., 0.024, the resultant n is 0.033-0.054,which matches up well with n = 0.050 used in Arnold et al. (1998) and n = 0.035 ± 0.027 estimated by Gleason et al. (2016).As such, we submit that our estimated bulk averaged value of n = 0.03-0.05may be a reasonable estimate for supraglacial streams-rivers under bare-ice conditions for use in surface meltwater routing models.

Is interfluve meltwater dominated by overland flow
or subsurface flow?
Our results suggest that during late July 2015, surface meltwater in the Rio Behar catchment was routed by a combination of slow interfluve flow (∼ 10 −3 -10 −4 m s −1 ) and fast open-channel flow (∼ 10 −1 m s −1 ).The latter RWF-inferred open-channel flow velocities correspond closely with field measurements (Karlstrom et al., 2014;Gleason et al., 2016), lending confidence that our bulk catchment-averaged values are grounded in a real-world process.Less clear, however, is the physical process governing interfluve flow.The RWF method we implement assumes that surface meltwater is routed by two distinct flow processes, i.e., open-channel flow and interfluve flow.The inferred mean interfluve velocity is ∼ 10 −3 -10 −4 m s −1 , which is 2-3 orders of magnitude smaller than the open-channel velocities (∼ 10 −1 m s −1 ).While RWF partitioning does not prescribe a physical transport process operating within interfluves, two likely candidate mechanisms are surface sheet flow and subsurface porous flow.
Sheet flow is overland (over ice in our case) flow taking the form of a thin, continuous film over relatively smooth surfaces and not concentrated into rills or channels (Mays, 2010).Manning's kinematic solution is generally used to analyze the sheet flow (Mays, 2010) and the meltwater flow velocity can be calculated as v s = f (n, M , S, L h ), where n is set to 0.05, M = 1.7 cm is daily effective melt in the Rio Behar catchment, S = 0.024 is the average catchment slope, and L h is set to 1-100 m.The resultant v s is ∼ 3-8 × 10 −2 m s −1 , which is similar to the terrestrial interfluve velocities (Moussa, 2008;Di Lazzaro, 2009;Singh et al., 2014;Mutzner et al., 2016) but still 1-2 orders faster than the ice surface v h quantified in this study.
This discrepancy suggests that interfluve transport is most likely controlled by subsurface meltwater flow, i.e., porous media flow.A fully saturated Darcy's law has been used in Arnold et al. (1998) and Leeson et al. (2012) (among many others) to describe meltwater routing on firn and snow sur- faces.However, to our knowledge, all models of meltwater routing over bare ice assume ice is impermeable and that Darcy's law is therefore not applicable (Arnold et al., 1998;Banwell et al., 2013;de Fleurian et al., 2016).Field studies do reveal that the bare-ice surface of ablating glaciers is often characterized by a porous ice layer termed "weathering crust" (Müller and Keeler, 1969;Fountain and Walder, 1998;Irvine-Fynn et al., 2011;Stevens et al., 2018), and lowdensity well-developed weathering crust has been observed in bare ice of the Rio Behar catchment (Cooper et al., 2018).Our results suggest that in contrast to current practice, principles of porous-media flow may be applied even in the bareice ablation zone if conditions of weathering crust and porous low-density bare ice are found.
The classic treatment for water transport through porous media is Darcy's law.Darcy's velocity (v d ) is defined as v d = kS/ϕ, where k is hydraulic conductivity (m s −1 ), S is slope (m m −1 ), and ϕ is weathering crust ice porosity; and hydraulic conductivity k is calculated as k = Kρ w g/µ, where K is absolute permeability (m 2 ), ρ w is water density (kg m −3 ), and µ is water viscosity (kg m s −1 ) (Arnold et al., 1998;Leeson et al., 2012).We followed Arnold et al. (1998) to set µ = 1.8 × 10 −3 kg m s −1 and Karlstrom et al. (2014) to set ϕ = 0.1 for weathering crust ice.Slope S is set as the mean slope (0.024) of the Rio Behar catchment.Near-surface ice permeability is highly uncertain, but applying the 10 −10 -10 −11 m 2 range estimated in Karlstrom et al. (2014), we estimate Darcy's velocity v d as 1.3 × 10 −4 -10 −5 m s −1 , 1 order smaller than the interfluve velocity v h we quantified.This implies that interfluve flow was not fully saturated in our study area, at least during the time when the ADCP supraglacial river discharge measurements were collected.This is consistent with the idea that subsurface flow through permeable weathering crust ice is only partially saturated, except perhaps in regions near channel heads that exhibit many interconnected small lakes and fully saturated slush.Partially saturated subsurface flow can be described by the Boussinesq equation (Bear, 1972), obtained by com-bining Darcy's law for porous flow with continuity of water, forced by meltwater recharge due to melting (Karlstrom et al., 2014).The resultant partially saturated (unconfined aquifer) velocity will be similar to or lower than the fully saturated velocity v d -as shown in the Appendix C and Fig. 8, reasonable values result in v h ∼ 10 −4 -10 −5 m s −1 .Because neither of these simple models for porous flow matches the inferred meltwater velocity v h = 10 −3 -10 −4 m s −1 in interfluve zones of the Rio Behar catchment, we suspect that multiple physical processes are involved in v h .For example, the combination of a relatively fast overland flow (∼ 10 −2 m s −1 ) and a slower porous subsurface flow (< 10 −4 m s −1 ), such as might occur for ephemeral channels on a variably saturated substrate, could explain the larger velocities.We leave mechanistic study of such issues to future work.

Advantages and limitations of RWF
The SRLF model is the first to attempt routing of surface meltwater downslope (Arnold et al., 1998).More recently, the Snyder SUH was used to derive moulin hydrographs (Smith et al., 2017).Both methods simulate observed moulin hydrographs reasonably well, but they cannot insightfully reveal the physical process of surface meltwater routing.Recently, permeable weathering crust was found on the Greenland bare-ice surface (Cooper et al., 2018), rather than impermeable bare ice as previously assumed (Arnold et al., 1998).
For this reason, it may not be appropriate to apply principles of supraglacial open-channel flow everywhere on the ice surface, i.e., subsurface flows may be more suitable for describing meltwater transport in the interfluve (hillslope) areas of higher-elevation ice separating meltwater channels.This reality calls for an easy-to-use, straightforward method to partition ice surface into channel vs. non-channel (i.e., interfluve) flow with each experiencing different physical flow processes.The RWF is our proposed solution for this partitioning.
We selected RWF over other SUH methods for the following reasons: (1) most SUH methods do not include interfluve (hillslope) transport and consider only the open-channel network on water routing (Singh et al., 2014), whereas RWF includes both hillslope and open-channel flows; (2) RWF is straightforward to implement and couple with remote sensing, requiring only hillslope and open-channel zones as inputs; (3) although RWF is a spatially lumped model, it can provide catchment-scale meltwater routing velocities, which are crucial for broad-scale understanding of ice surface hydrology.The derived mean open-channel velocity is comparable to field-measured velocities in small supraglacial streams, and the derived hillslope velocity is comparable to simulations of a partially saturated subsurface hydrological model.Therefore, RWF appears to be a simple and useful tool for modeling meltwater routing across broad-scale areas of melting ice.
The central hypothesis of UH theory is that catchment response to rainfall (here, melt production) is linear, i.e., variations in input rainfall and melt change only the ordinates, not duration, of the direct hydrograph (Dingman, 2015).Meltwater routing is also assumed to be fully determinable by the morphometric characteristics of the catchment (here, IDC drainage networks, shape, area, etc.) (Singh et al., 2014).These assumptions of linearity and fixed basin response simplify routing models but also create some limitations.Particularly, hydraulic geometries of meltwater channels are not independent of IDC characteristics but instead vary nonlinearly with channel discharge (i.e., Q = wdv, w = aQ b , d = cQ f , v = kQ m ).Gleason et al. (2016) suggest that supraglacial meltwater channels primarily accommodate greater discharges by increasing v c (m = 0.63-0.95),which is also supported by Brykała (1999) (m = 0.49).Thus, at higher discharges open-channel flow velocities will increase nonlinearly.Moreover, v h is also affected by surface melt dynamics (Leeson et al., 2012;Karlstrom et al., 2014;Karlstrom and Yang, 2016), and spatiotemporal patterns in surface meltwater production surely influence IDC hydrographic responses, just as spatiotemporal variations in rainfall pattern impact terrestrial hydrographic responses (Nicótina et al., 2008).Therefore, supraglacial IDCs may not respond to surface melt linearly.Future studies should consider varying v h and v c based on different surface melt patterns, which should generate variable RWFUHs during a melt season.
In addition to spatiotemporal variations in meltwater input, v h and v c surely vary spatially as well.RWF is a spatially lumped hydrologic model yielding fixed constants of v h and v c averaged across the catchment (Table 1).However, in reality, v h and v c vary diversely within a catchment (Maidment et al., 1996;D'Odorico and Rigon, 2003).The SRLF model (Arnold et al., 1998) offers spatially varying surface meltwater routing velocities based on ice surface topography much like distributed terrestrial hydrologic models (Maidment et al., 1996;Liu et al., 2003), thus providing a physically based approach to investigate spatially varying meltwater routing velocities.Therefore, combining RWF and SRLF would be a promising future direction for producing a physically based, spatially distributed surface routing model for use on ice sheets.As a starting point, RWF-derived interfluve travel velocities could be included in the SRLF model to parameterize meltwater flows through bare, porous low-density bare ice (Cooper et al., 2018), especially for partially saturated subsurface conditions.We leave spatially distributed routing models for future studies because of two reasons: first, these models need more data inputs and parameters (which are difficult to estimate) than RWF; second, we need to determine what additional scientific value would be gained from more complex models.
Although the four RWF-based approaches presented here (nonconservative mapped rivers, conservative mapped rivers, high A c threshold, low A c threshold; see Fig. 6) all simu-late the Rio Behar moulin hydrograph reasonably well, it is important to interpret each method from a physically based standpoint and "get the right answers for the right reasons" (Kirchner, 2006).The four RWF-based approaches all perform well because we calibrated v h and v c from field and remote-sensing observations and then used the optimal velocity combination to recreate a measured hydrograph.However, these bulk calibrated v h and v c values may or may not be reasonable estimates for channel and interfluve velocities more broadly.They were collected using field and remotesensing observations collected during a narrow window of time, from 21 to 23 July 2015 when snow was gone and the Rio Behar IDC supraglacial stream-river network was highly developed on fully bare ice.We are encouraged that our derived values agree broadly with other field studies (McGrath et al., 2011;Chandler et al., 2013), but additional field investigations are needed to confirm that the mean values of v h and v c derived here may be usefully applied to other bare-ice locations on the ice sheet.

Field site and observation recommendation
Selection of an IDC for field study is logistically challenging and requires careful planning and design.We selected the Rio Behar catchment by considering surface melt intensity, distance to ice edge, distance to automatic weather stations, catchment size and shape, catchment outlet (moulin) conditions, and safety conditions (Smith et al., 2017).Two types of field measurements will be crucial for better understanding of surface meltwater routing process: supraglacial river discharge and subglacial water pressure.Supraglacial river discharge hydrographs can be used to validate the performance of surface meltwater routing methods, while subglacial water pressure can be used to estimate the hydrological responses of subglacial environments to different supraglacial meltwater inputs (moulin discharge).
In situ investigation is also necessary to characterize interfluve conditions.Cooper et al. (2018) analyzed the density and hydrological properties of bare, ablating ice away from open channels by drilling holes into wet bare ice and measuring the subsurface porosity and water infilling rate, properties that cannot be measured from remote sensing.Satellite images are certainly useful for providing preliminary observations for ice surface conditions.For example, Smith et al. (2017) partitioned bare ice and snowpack zones using high-resolution satellite imagery and Ryan et al. (2018) investigated ice surface albedo, surface impurities, and cryoconite holes using higher-resolution UAV images.That said, we are unaware of any remote-sensing solution to confirm the presence or absence of saturated subsurface weathering crust and its hydraulic conductivity, so field measurements remain essential at present.

Future research directions
It is crucial to couple surface meltwater routing models with subglacial hydrological models to build a complete understanding of surface-to-bed meltwater connections.One path forward would be to use SUH, SRLF, and/or RWF to calculate moulin hydrographs using DEMs of different sources and spatial resolutions, then coupling this output to the Subglacial Hydrology and Kinetic, Transient Interactions (SHAKTI) subglacial hydrology model (Sommers et al., 2018).Doing so would allow derivation of hourly changes in subglacial water pressure in response to different moulin discharge inputs.A logical next step would be to then analyze the potential impact of these varying subglacial water pressures on subglacial hydrologic system evolution and ice flow dynamics.An ultimate objective should be to model the complete surface-to-bed meltwater transfer process by using RCMs to generate surface melt, surface routing to generate moulin discharge hydrographs, and subglacial models to track basal water pressure, subglacial hydrological system evolution, and ice flow.
Crucial ice surface topographic characteristics, such as slope, flow direction, flow length, drainage area, and drainage networks, are scale dependent.Zhang and Montgomery (1994) illustrated DEM resolution significantly impacts hydrological responses of terrestrial catchment to rainfall, using 2, 4, 10, 30, and 90 m DEMs.We suggest that DEM source and catchment morphometry both affect a DEM's capability of simulating meltwater routing on the ice surface.In general, a 100 m or coarser-resolution DEM may yield larger offsets in simulating moulin hydrographs compared to a 30 m resolution DEM but the specific offsets need further estimation.Moreover, high-resolution Arc-ticDEM (Noh andHowat, 2015, 2017) raises prospects for studying meltwater routing in unprecedented detail and it covers the entire Greenland Ice Sheet.The ArcticDEM products are now released at 2, 10, 32, 100, 500, and 1000 m resolution (Release 7, September 2018).Therefore, we recommend using ArcticDEM products in future meltwater routing studies.
A particularly promising area for future work will be incorporating surface routing delays in studies of proglacial discharge, in order to remove the effects of supraglacial delays before interpreting subglacial delays and/or storages from proglacial river discharge hydrographs.For example, in southwestern Greenland, numerous supraglacial IDCs form on the ice surface each summer and route meltwater into moulins (Thomsen et al., 1989;Banwell et al., 2012Banwell et al., , 2013;;Arnold et al., 2014;Yang and Smith, 2016).By integrating RWF surface routing with the lumped supra-, en-, and subglacial delays obtained from proglacial river studies (Rennermalm et al., 2013;Van As et al., 2014, 2017), subglacial runoff delays can be better separated from supraglacial delays.The corrected subglacial delays could then be used to better interpret the coupling of surface melt-runoff with sub-glacial water pressures (van de Wal et al., 2015;Wright et al., 2016) and used more generally to investigate the evolution of the subglacial hydrological system over short timescales (Banwell et al., 2013;Andrews et al., 2014).As a simple example, the supraglacial delay we obtain here for the Rio Behar catchment using RWF is ∼ 0.5 days (11 h), whereas the lumped delay obtained from Van As et al. (2017) is 3.9 days.This suggests subglacial delays contributed ∼ 90 % of the total delay between runoff generation on the ice sheet and its appearance in proglacial river discharge at the ice edge.

Conclusions
Numerous internally drained catchments (IDCs) are distributed on the southwest Greenland Ice Sheet.These catchments collect and drain meltwater through supraglacial stream-river networks to large, terminal outlet moulins, consequently constraining the timing and discharge of meltwater flowing over the ice surface to moulins that deliver meltwater to discrete locations on the bed.A growing literature is recognizing the nontrivial influence of supraglacial meltwater transport processes on meltwater received at the bed and proglacial zone.This study has investigated surface meltwater routing processes in the Rio Behar catchment, a moderately sized IDC near Kangerlussuaq, using high-resolution satellite and UAV observations and an in situ field-measured supraglacial river discharge hydrograph collected in late July 2015.Key meltwater routing parameters were quantified using a rescaled width function (RWF) surface routing model, which distinguishes fast meltwater transport through supraglacial stream-river channels from slow transport over interfluves, the latter likely involving partially saturated, near-surface porous-media flow.Our main contribution is thus to partition interfluve vs. open-channel flow in a surface routing model adapted for use on ice surfaces, using remote sensing and sparse in situ measurements.The model includes two terrestrial hydrologic techniques (unit hydrograph and rescaled width function) in a simple and flexible approach.Its scientific utility includes quantifying runoff delays on the Greenland Ice Sheet, improved understanding of open-channel versus interfluve water transport on bare melting ice, improved interpretation of proglacial river discharge hydrographs, and improved coupling of RCM datasets with subglacial borehole studies and models of subglacial hydrology and ice dynamics.the resultant mean channel width is 2.5 ± 2.0 m, supporting our conclusion that numerous small supraglacial streams control bulk-catchment channel velocity v c .We defined large supraglacial rivers as the features that can be identified by moderate-resolution (10-30 m) satellites (e.g., Sentinel-2 and Landsat 8), while small supraglacial streams are the features that can only be identified by high-resolution (0.5-2.0 m) satellites (e.g., WorldView-1/2/3/4).It is subjective to determine a threshold width but if required, we recommend 10 m.

Figure 1 .
Figure 1.The Rio Behar catchment is a moderately sized internally drained catchment (IDC) on the southwest Greenland Ice Sheet.A highly developed supraglacial stream-river network was mapped from a high-resolution (0.5 m) panchromatic WorldView-1 (WV1) image acquired on 18 July 2015.The catchment outlet moulin is under the black star and 5 m MAR grid cells are shown in grey rectangles.A derivative 3 m resolution DEM is used to delineate the topographic catchment boundary (black line).An in situ hourly hydrograph of supraglacial river discharge collected by Smith et al. (2017) (black star) was used to calibrate a Rio Behar catchment unit hydrograph (UH) for the rescaled width function (RWF) and surface routing and lake-filling (SRLF) models.

Figure 4 .
Figure 4. Meltwater travel distances for (a) open channels and (b) interfluves of the Rio Behar internally drained catchment (IDC) as obtained by burning the 18 July 2015 WorldView DEM with the conservative remotely sensed active supraglacial stream-river surface drainage network.

Figure 5 .
Figure 5. Calibration of channel and interfluve velocities using (a) conservative and (b) nonconservative supraglacial river delineations.Contour labels show Nash-Sutcliffe model efficiencies (NSEs) obtained by applying different combinations of mean interfluve (v h ) and open-channel velocities (v c ) and comparing with the field-measured hydrograph just upstream of the Rio Behar catchment terminal outlet moulin.

Figure 6 .
Figure 6.(a) Unit hydrographs (UHs) and (b) simulated direct hydrographs at the terminal moulin of the Rio Behar catchment, as modeled by different data sources and methods.The MAR UH is calibrated from effective MAR melt and measured supraglacial river discharge (Fig. 2).The rescaled width function (RWF) and surface routing and lake-filling (SRLF) UHs are obtained by calibrating open-channel and interfluve velocities to optimally match the MAR UH.

Figure 7 .
Figure 7. Simulated Rio Behar catchment discharge at the terminal outlet moulin for 20 July to 31 August 2015 assuming a seasonally evolving supraglacial stream-river network.A series of contributing area thresholds (A c ) of 250, 500, 1000, 2500, and 5000 pixels is used to mimic an evolving supraglacial stream-river drainage network from the 18 July 2015 WorldView DEM.The minimum A c (250 pixels) is used to simulate well-developed river networks, while the maximum A c (5000 pixels) is used to simulate poorly developed river networks.Variable A c values are used to simulate dynamic river networks, which are considered to best capture realistic seasonal evolution.Drainage density is lowest in the early and late season and highest in July.The diurnal variations in moulin discharge are strongly controlled by supraglacial stream-river network patterns.

Figure 8 .
Figure8.Schematic diagram of subsurface meltwater porous flow through weathering crust and permeable near-surface low-density ice to estimate interfluve transport speed v h .Melt generated at the near surface at rate M percolates through porous ice, supplying a water table that transports melt downhill towards streams at heights h 0 and h L above a layer of impermeable ice (dashed line).

Figure D1 .
Figure D1.Channel width histogram of conservatively mapped supraglacial stream-river networks in the Rio Behar catchment.Most supraglacial meltwater channels are narrower than 4 m with a mean width of 2.5 ± 2.0 m, confirming that numerous small supraglacial streams dominate the bulk-catchment average channel velocity v c .

Table 1 .
Surface meltwater travel distances, velocities, and times in the Rio Behar catchment for July 2015.