Inter-comparison of surface meltwater routing models for the Greenland Ice Sheet and influence on subglacial effective pressures

. Each summer, large volumes of surface meltwater drain off the Greenland Ice Sheet surface through moulins to the bed, impacting subglacial hydrology and ice flow dynamics. Supraglacial surface routing delays may propagate to englacial and subglacial hydrologic systems, requiring accurate assessment to correctly estimate subglacial effective 15 pressures. We compare hourly supraglacial moulin discharge simulations from three surface meltwater routing models: synthetic unit hydrograph (SUH), bare ice component of surface routing and lake filling (SRLF), and rescaled width function (RWF), for four internally drained catchments on the southwestern Greenland ice sheet surface. The routing models are forced identically using surface runoff from the Modèle Atmosphérique Régionale regional climate model (RCM). For each catchment, simulated moulin hydrographs are input to the SHAKTI subglacial hydrologic model to simulate diurnally 20 varying subglacial effective pressure variations in the vicinity of a single moulin. Overall, all three routing models produce more realistic moulin discharges than simply using RCM runoff outputs without surface routing, but with significant differences in peak moulin discharge and time-to-peak. In particular, RWF yields later, smaller peak moulin discharges than SUH or SRLF due to its representation of slow interfluve flow between supraglacial channels, and can readily accommodate the seasonal evolution of supraglacial stream/river networks. Differences among the three routing models are reflected in a 25 series of simple idealized subglacial hydrology simulations that yield different diurnal effective pressure amplitudes; however, the supraglacial hydrologic system acts as short-term storage for surface meltwater and the temporal mean effective pressure is relatively consistent across routing models. on subglacial effective pressure (ice overburden pressure – water pressure) in the vicinity of a moulin in an idealized subglacial setting. Results show that inclusion of surface meltwater routing introduces significantly smaller peak moulin discharges and delayed time-to-peak relative to the RCM instantaneous runoff. Different surface meltwater routing models, as well as different spatial-resolution 25 DEMs and seasonal evolution of the supraglacial stream/river networks, induce variable diurnal moulin discharges and effective pressures, which in reality would influence ice sliding velocity. While the different routing models produce different diurnal amplitudes in effective pressure variation, the qualitative channelization behaviour and temporal mean effective pressures are relatively consistent across surface meltwater routing models. Together, these findings urge caution for better representations of surface meltwater routing and moulin discharge simulations to drive subglacial hydrology and 30 ice dynamics models, as well as highlight the need for further research to investigate the cumulative effects of diurnal inputs to the subglacial drainage system and the relevant impacts on ice dynamics through detailed, large-scale modelling.

substantially different diurnal effective pressure amplitudes depending on the applied surface meltwater routing model; however, display relatively consistent mean effective pressure across routing models.

Introduction
Large volumes of meltwater are routed through supraglacial stream/river networks on the Greenland ice sheet each summer (Smith et al., 2015). In temperate areas of the ice sheet, most of this surface meltwater is injected to the bed via moulins 5 (Catania et al., 2008;Lampkin and VanderBerg, 2014;Smith et al., 2015;Yang and Smith, 2016;Koziol and Arnold, 2018), where it can modulate ice flow (Bartholomew et al., 2011;Palmer et al., 2011;Banwell et al., 2013;Hewitt, 2013;Andrews et al., 2014;de Fleurian et al., 2016). However, the role of the supraglacial system in controlling subglacial hydrology remains poorly studied to date (Flowers, 2018). In particular, there are limited constraints on the efficiency of surface meltwater routing (Smith et al., 2017), resulting in large uncertainties in hydrodynamic coupling with ice motion. 10 In most studies that investigate surface-to-bed meltwater connections and the behaviour of the subglacial hydrologic system, surface meltwater routing is either simplified or simply ignored, i.e., moulin hydrographs are estimated directly from output of regional climate models (RCMs). In these instances, RCM instantaneous grid cell runoff is simply summed over internally drained catchments (IDCs, "hydrologic units on the GrIS surface that collect and drain meltwater through supraglacial stream/river networks to terminal moulins or lakes") (Yang and Smith, 2016) to estimate supraglacial moulin 15 discharge or to fill supraglacial lakes (McGrath et al., 2011;Fitzpatrick et al., 2014). This summed value is used to drive models of subglacial hydrologic system evolution and/or ice flow dynamics. Such a simplification may be appropriate for small IDCs or for very long-term studies; however, surface meltwater routing has been found to substantially modify ice surface runoff (Banwell et al., 2012;Banwell et al., 2013;Arnold et al., 2014), by altering the magnitude and timing of peak moulin discharge (Smith et al., 2017;Yang et al., 2018). As such, surface meltwater routing has strong potential to affect 20 subglacial hydrologic system evolution and hydrodynamics, especially for large IDCs and short (i.e., diurnal) time scales.
The current lack of representation of the surface meltwater routing leads to an insufficient understanding of surface-to-bed meltwater connections and ice dynamics, particularly on diurnal timescales. Therefore, constraints on IDC discharge can provide critical boundary conditions for studies of the subglacial hydrologic system.
In terrestrial systems, including ice sheet and glacier surfaces, water routing may be characterized by transport velocity 25 (v), flow length (L), and transport time (t). Meltwater transport velocity and flow length can be estimated from ice surface topography (Arnold et al., 1998), with the distribution of transport times representing the hydrologic response of a particular supraglacial catchment to surface melt . Alternately, empirical unit hydrograph (UH, "a transfer function that is widely used for modelling catchment runoff response to rainfall events for some unit duration and unit depth of effective water input" (Smith et al., 2017)) methods can be used to estimate the hydrologic response based on catchment 30 shape and area alone, enabling derivation of moulin hydrographs through convolution of input RCM surface runoff with remotely-sensed UH parameters.
In order to route supraglacial meltwater on ice surfaces, Arnold et al. (1998) and Banwell et al. (2012) developed the DEM-based Surface Routing and Lake Filling (SRLF) model. SRLF assumes all meltwater is transported downslope, with Manning's open-channel flow equation used to calculate spatially variable meltwater transport velocities for every bare-ice DEM pixel.
More recently, the Snyder Synthetic Unit Hydrograph (SUH) (Snyder, 1938) was adapted for use on ice sheets to 5 estimate moulin hydrographs without the need for a DEM (Smith et al., 2017). To do this, Smith et al. (2017) calibrated two key SUH parameters using a field-measured moulin hydrograph, and a Gamma function was then used to build individual SUHs for hundreds of remotely-sensed internally drained catchments (IDCs) (Smith et al., 2017). Because the Snyder SUH is a simple lumped model that does not rely on DEMs, it provides a straightforward but static approach to route surface meltwater to moulins. 10 Most recently, Yang et al. (2018) utilized Rescaled Width Function (RWF) (D'Odorico and Rigon, 2003) to further partition the bare ice surface into either interfluve (i.e., hillslope) or open-channel zones using high-resolution satellite imagery to discern between the two. Catchment-averaged meltwater transport velocities for each zone were then calibrated using a field-measured moulin hydrograph (Smith et al., 2017), with slow interfluve flow (~10 -3 -10 -4 m/s) and fast openchannel flow (~10 -1 m/s) combined to simulate the downstream moulin hydrograph. 15 These three surface meltwater routing approaches (SRLF, SUH, and RWF) use different assumptions and data requirements. SUH is the simplest model; it requires catchment shape, area, and several empirically-derived parameters to estimate surface meltwater transport time (Smith et al., 2017). In contrast, RWF requires a partitioning of slow (interfluve) and fast (open channel) flow that can substantially impact the performance of RWF without careful spatial and temporal calibration. Moreover, SRLF and RWF rely on surface DEMs to calculate meltwater flow paths  and 20 SRLF also relies on DEMs to calculate meltwater routing velocities (Arnold et al., 1998). Characteristics such as slope, flow direction, flow length, drainage area, and drainage networks are readily extracted from DEMs but are scale-dependent (Montgomery and Foufoula-Georgiou, 1993), signifying that DEM spatial resolutions influence their computation (Zhang and Montgomery, 1994;Hancock et al., 2006). This DEM resolution dependence can now be tested with the advance of highresolution DEM datasets (e.g., 2 m ArcticDEM) of the GrIS surface (Noh and Howat, 2015, 2018. 25 To assess supraglacial discharge variability among current surface meltwater routing models, this study simulates surface meltwater runoff through four moderate-size IDCs using the SUH, SRLF, and RWF models described above ( Figure   1).Routing model ability to simulate meltwater transport velocity, flow length, transport time, UH, peak moulin discharge, and time to peak are intra-compared, and the impact of DEM spatial resolution on SRLF and RWF is examined. For all of these models, the consequent effect of using them to drive subglacial effective pressure variations is investigated for the 30 illustrative case of meltwater input to the bed via a single moulin, using the Subglacial Hydrology and Kinetic, Transient Interactions (SHAKTI) subglacial hydrologic model (Sommers et al., 2018). From this, we demonstrate the impact of surface meltwater routing on subglacial effective pressure fluctuations.
Notably, this study cannot, in good faith, focus on method comparison, nor determine the "best" (i.e., best able to reproduce a real-world moulin hydrograph) due to the lack of calibration and validation data. Owing to this limitation, the goal of this study is to assess differences among the three surface meltwater routing models, rather than to determine which model most realistically simulates surface meltwater routing on the ice surface. By using the outputs from all three routing models as meltwater inputs to drive the SHAKTI subglacial model, we seek to characterize the impact of differences in 5 surface routing on subglacial pressures and evolution, particularly over diurnal timescales. Nevertheless, these results demonstrate the extent to which the choice of surface meltwater routing models can alter modelled subglacial conditions. We conclude with a general discussion of surface meltwater routing on the Greenland ice surface and include some recommendations for best practices and future research.

Study area and data sources
We selected four moderate-size supraglacial IDCs in the Russell Glacier region, southwestern GrIS to explore the impact of SUH, SRLF, and RWF meltwater routing models on moulin discharge and subglacial effective pressure ( Figure 1, Table 2). They are distributed at approximately 200 m elevation intervals in order to span the elevational range of most welldeveloped IDCs found in the Russell Glacier region and the variable surface melt conditions of this region (Fitzpatrick et al., 15 2014;Yang and Smith, 2016). Large supraglacial lakes are absent in these four IDCs and surface meltwater is all routed to the moulin at the catchment outlet ( Figure 1). As such, surface runoff produced in each IDC should equal to the moulin discharge (Smith et al., 2017). The four IDC areas range from 53.0 km 2 to 66.9 km 2 ; surface elevations vary from 1086 m to 1672 m and ice thicknesses from 854 m to 1432 m (Table 2). A moulin discharge hydrograph collected at Rio Behar catchment (IDC2 in this study; 67.049346N, 49.025809W) for 72 h from 20 to 23 July 2015 was used to calibrate parameters 20 of SUH and RWF models (see Sections 2.5 and 2.6). It is problematic to apply these empirically-derived parameters over large spaces and long times . Therefore, the selected IDCs are distributed in a relatively small region and the areas of IDC1, IDC3, and IDC4 are similar to the Rio Behar catchment (IDC2).
The high-resolution (2 m) ArcticDEM is the primary data source for our study of surface meltwater routing in the four supraglacial IDCs. ArcticDEM products are created from high-resolution (0.5 m) WorldView-1/2/3 stereo images by Polar 25 Geospatial Center (PGC) and provide an unprecedented opportunity to investigate ice surface landscape and hydrologic processes (Noh and Howat, 2015, 2018. MAR (Modè le Atmosphé rique Ré gionale) 3.6, a regional coupled atmosphere-land climate model with 20 km native horizontal resolution (Fettweis et al., 2013), provides hourly surface runoff simulations.2.2 Regional climate model runoff simulations Hourly simulations of surface runoff (R) during July 2015 were obtained from MAR 3.6 (Smith et al., 2017) in order to 30 be to be coincident with the field-collection of the moulin discharge hydrograph. Supraglacial IDC boundaries were derived from ArcticDEM following the method of Karlstrom and Yang (2016). Catchment-averaged hourly R [mm/h] was obtained by clipping MAR grid pixels with the supraglacial IDC boundaries and summing their corresponding runoff values (Smith et al., 2017). Absent any routing, the resulting area-integrated R is termed RCM instantaneous runoff.

Unit Hydrograph
A surface meltwater routing model does not need to explicitly determine how meltwater produced on a cell is routed to its downstream grid cell(s) in a catchment, which is the aim of spatially-distributed routing models (Singh et al., 2014). In 5 contrast, a routing model can be lumped and only determines how surface meltwater produced in the catchment is temporally routed to the catchment outlet. Unit hydrograph (UH) is designed for this purpose and provides a transfer function between surface runoff R and moulin discharge Q (i.e., direct hydrograph): where * is the convolution operator. For each IDC, convolution of UH with hourly R yields a discharge hydrograph at the 10 catchment outlet (i.e., moulin) and thereby determines how surface meltwater produced in the catchment is temporally routed to the catchment outlet. Therefore, the three surface meltwater routing methods (SUH, SRLF, and RWF) based on UH can be considered routing models.

Snyder Synthetic Unit Hydrograph
Snyder Synthetic Unit Hydrograph (SUH) (Snyder, 1938) assumes catchment morphometry controls the shape of UH 15 and utilizes the catchment main-stem stream length (L, in km), the distance from the catchment outlet (here, the terminal moulin) to the point on the main channel nearest to the catchment centroid (Lca, in km), and two parameters, Cp and Ct, to determine peak discharge (hp, in hr -1 ) and time-to-peak (tp, in hours) of SUH as tp = Ct(LLca) 0.3 , hp = Cp/tp. Cp and Ct are two parameters depending on "units and drainage-basin characteristics" (Snyder, 1938). Smith et al. (2017) used a moulin discharge hydrograph collected at the Rio Behar catchment during 20 to 23 July 2015 to calibrate Cp and Ct . Following 20 Smith et al. (2017), we use Cp = 0.72 and Ct = 1.61 to determine hp and tp of SUHs and employ a Gamma distribution to determine the catchment-specific SUH for each IDC.
Although Smith et al. (2017) demonstrated parameter transferability using two other independently field-measured moulin hydrographs, the two calibrated parameters, Cp and Ct, are collected at one moulin for one moment in time so parameter transferability over large regions and long timescales may still be limited. However, it is reasonable to apply them 25 for the Russell Glacier region, southwestern GrIS during July 2015.

Surface Routing and Lake Filling
SRLF has been used to fill supraglacial lakes within supraglacial catchments (Banwell et al., 2012). Here we use the model to route surface meltwater downslope to a moulin and subsequently to the ice sheet bed. The SRLF model employs Darcy's law to route surface meltwater flow through snow and Manning's open-channel flow equation to route meltwater 30 flow over bare-ice surfaces (Arnold et al., 1998). SRLF is the first surface meltwater routing model for the GrIS and has been successfully applied to simulate supraglacial lake growth and to drive subglacial hydrological evolution during several complete melt seasons (Banwell et al., 2012;Banwell et al., 2013;Arnold et al., 2014;Banwell et al., 2016;de Fleurian et al., 2016). In this study, we use its bare-ice component to make it comparable with the SUH and RWF (Smith et al., 2017;Yang et al., 2018) and representative of mid melt season conditions in the ablation zone. SRLF assumes universal presence of supraglacial open-channel flow everywhere on bare ice surface, by applying Manning's open-channel flow equation 5 (Manning, 1891) to calculate channelized meltwater velocity v for every pixel along a DEM-derived flow path: where RH is the hydraulic radius of the supraglacial meltwater channel, S is water surface slope calculated from DEM, and n is the Manning roughness parameter (Arnold et al., 1998). Following Arnold et al. (1998), we set RH = 0.035 m and n = 0.05.
We calculated S per pixel using ArcticDEM. The shortest meltwater flow length (L) was determined from ArcticDEM using 10 the D8 flow routing algorithm following Karlstrom and Yang (2016). Meltwater transport time (t) for each DEM pixel was then determined as t = L/v. To compare with synthetic unit hydrograph and rescaled width function, a one-hour unit hydrograph was calculated by hourly binning the transport time raster , with the resultant UH termed SRLF UH.

Rescaled Width Function 15
The rescaled width function (RWF) routing model partitions surface meltwater transport pathways into interfluve (i.e., hillslope) and open-channel distances (Lh and Lc), to determine two catchment-averaged meltwater velocities (vh and vc) to represent routing efficiency over these regions (D'Odorico and Rigon, 2003). Then, the transport time for each pixel in a catchment may be calculated as: (2) 20 where th is the interfluve travel time and tc is the channel travel time. Similar to SRLF, the shortest meltwater flow path for interfluve and channel transport distances is determined using a DEM (Karlstrom and Yang, 2016), with the relative coverage determined by the extent of supraglacial streams and rivers . By applying the constant interfluve velocity vh to interfluve pixels and channel velocity vc to open channel pixels, the meltwater transport time from every catchment pixel to the moulin can be determined and a consequent catchment UH generated, hereafter referred to as RWF 25 UH. The primary difference between RWF and SRLF is the way they determine meltwater flow velocity for each catchment cell. For RWF, using the in situ moulin hydrograph of Smith et al. (2017), Yang et al. (2018) determined the velocity combination of vh = 0.0006 m/s and vc = 0.4 m/s to be the optimal match to field observations. Although vh and vc may change notably during the melt season because of variable surface melt intensity and ice surface topography, it is reasonable to assume that vh and vc are persistent during a relatively short time period for similar hydrologic and glaciological 30 environments. Therefore, in this study, the optimal velocity combination obtained in  is assumed to be transferable and used to create RWF UHs for the four IDCs.

Controls of DEM resolution on surface meltwater routing
DEM spatial resolution has important impacts on the hydrologic response of terrestrial catchments to precipitation (Zhang and Montgomery, 1994;Hancock et al., 2006). To estimate the control of DEM spatial resolution on ice surface meltwater routing, 2 m resolution ArcticDEM data were resampled to 5 m, 10 m, 30 m, and 90 m resolutions, similar to Zhang and Montgomery (1994), and the corresponding meltwater flow paths, velocities, transport time, UHs, and moulin 5 hydrographs for SRLF and RWF subsequently calculated. For SRLF, all five resampled DEMs were used. For RWF, only 2 m, 5 m, and 10 m resolutions were used to prevent DEM resolution from exceeding the typical interfluve transport distance, ~10 m for southwestern GrIS ; otherwise, interfluve transport distance would be overestimated and the resultant hydrograph would be inappropriate (Hancock et al., 2006).

Temporal evolution of surface meltwater routing 10
The spatial pattern of supraglacial stream/river networks is known to vary significantly during a melt season (Lampkin and VanderBerg, 2014). This affects the hydrologic response (and UH) of a supraglacial IDC as transport distances change (Montgomery and Foufoula-Georgiou, 1993; van Meerveld et al., 2019). RWF can mimic this effect by varying the partitioning of interfluve versus open-channel zones . SRLF, in contrast, assumes the bare ice surface has a stable response to the surface melt and uses static meltwater routing velocities to build the UH. SUH relies on Cp and Ct to 15 build the UH. If a moulin hydrograph is available during the entire melt season, we can calculate time evolving Cp and Ct using variable time-to-peak (tp), peak discharge (hp), and main-stem channel length (L), thus creating multi-temporal UHs.
Without such direct measurements, the parameters cannot be realistically varied in time and SUH cannot mimic variable hydrologic response of a catchment to surface melt.
Capturing the rapid seasonal evolution of supraglacial stream/river networks is challenging (Lampkin and VanderBerg, 20 2014;Yang et al., 2017) and remains largely unperformed (Flowers, 2018), but a positive relationship exists between surface melt and supraglacial drainage density (Yang et al., 2017). As a simplified simulation, this study uses a series of cumulative contributing area (or channel initial threshold, Ac) values (Ac = 100 m 2 , 250 m 2 , 500 m 2 , 1000 m 2 , 2500 m 2 , 5000 m 2 ) to simulate the hypothetical seasonal evolution of supraglacial stream/river networks by evolving the partitioning of interfluve vs. open channels and resultant surface meltwater routing simulations . Ac defines the surface area needed 25 to initiate open channel flow. Larger Ac values will yield smaller open-channel zones because larger contributing interfluve areas are required to form open channels (King et al., 2016). Smaller values of Ac signify periods of intense surface melt and actively flowing, well-developed supraglacial stream/river networks (Smith et al., 2017), whereas larger values signify poorly developed network characteristics of the late melt season . A total of six 5-day RWF UHs were created for each Ac value, broadly consistent with a decreasing RCM runoff trend during July 2015. Each RWF UH was 30 conducted to calculate moulin discharge for five days and the resultant moulin discharge is termed as dynamic Ac discharge.

Subglacial hydrology modelling
A variety of numerical models have been formulated to simulate different aspects of the subglacial drainage system and compute subglacial effective pressure (defined as the difference between ice overburden pressure and subglacial water pressure) (Flowers, 2015;de Fleurian et al., 2018). To examine the variations in modelled effective pressure resulting from the different surface meltwater routing models considered in this study, we treat their respective moulin discharges as direct 5 meltwater inputs to the bed via a single moulin, and simulate subsequent evolution of the local subglacial drainage system using the SHAKTI model (Sommers et al., 2018).
SHAKTI is built within the Ice Sheet System Model (Larour et al., 2012). Using finite elements, it applies a single set of equations over the entire model domain, conserving mass, energy, and momentum, with a subglacial water flux formulation allowing for local development of both turbulent and laminar flow regimes, as well as regimes corresponding to 10 the broad transition between laminar and turbulent. SHAKTI calculates transient effective pressure, hydraulic head, subglacial gap height, subglacial water flux, and "degree of channelization" (defined as the ratio of the rate of opening by melt to the rate of opening by sliding over bumps in the bed). Subglacial effective pressure and geometry evolve naturally to produce continuous configurations ranging from sheet-like to channelized drainage. A thorough description of the SHAKTI model equations, methods, and features, including limitations, can be found in Sommers et al. (2018). 15 To examine the influence of different surface meltwater routing models on local subglacial hydrology, we apply meltwater input at a single moulin at the center of a 1 km square domain. The domain is discretized with an unstructured triangular mesh consisting of 614 elements and a typical element edge length of 50 m. Outflow is to the left edge of the square domain, with a Dirichlet boundary condition setting the subglacial water pressure to 50 % of the ice overburden pressure (a reasonable assumption for our drainage catchments far from the ice sheet margin), and zero-flux pressure 20 boundaries on the other three sides. For each IDC, the mean surface slope and mean ice thickness are used (Morlighem et al., 2017), with the bed slope set equal to the surface slope to create a uniformly thick tilted slab of ice, and sliding velocity is prescribed as 100% of the annual mean observed surface velocity (Table 3). Constants and parameter values used in the simulations are summarized in Table 3.

Simulations of supraglacial moulin discharge
Inclusion of surface meltwater routing clearly modifies the timing of peak moulin discharge relative to RCM runoff obtained from MAR, with each of the three routing models performing somewhat differently. When RCM instantaneous runoff alone is used to simulate moulin discharge (i.e., no routing), moulin discharges peak between 13:00-15:00 local time for all the four study IDCs ( Table 2). Inclusion of routing processes, however, introduces considerable delay between the 30 timing of peak RCM runoff generation and its arrival at the moulins, with all routing models displaying peak moulin discharges occurring from 17:00 to 23:00 local time (  Figures 2 and S1-S4). SUH-routed hydrographs simulate peak moulin discharge between approximately 19:00 and 21:00, a delay of 6 hours compared to the RCM instantaneous peak discharge, while SRLF-routed delays vary from 3 hours to 10 hours and RWFrouted delay varies from 7 hours to 10 hours (Table 2). These delays are longer than ~2 hours measured for small IDCs near 5 ice margins (Shepherd et al., 2009;McGrath et al., 2011;Bartholomew et al., 2012) but are comparable to ~5-7 hours measured for the ~60 km 2 Rio Behar catchment, confirming the control of catchment size on surface meltwater routing delays (Smith et al., 2017).
Surface meltwater routing also controls peak moulin discharge magnitude and diurnal discharge range. This observation is clearly demonstrated when examining the derived unit hydrographs (UHs) for each routing model. The UH isolates the 10 impact of basin characteristics on moulin discharge by providing the basin integrated runoff resulting from one cm of surface melt over one hour. The three UHs (i.e., SUH, SRLF UH, and RWF UH) and their routed moulin discharges for the four IDCs are presented in the supplement ( Figures S1-S4); the resultant patterns of moulin hydrographs are similar so we use IDC1 to illustrate our results here. For IDC1, the peak values of RWF UHs are lower than SUHs and 2 m SRLF UHs.
Smoother UHs yield smaller peak moulin discharge magnitudes and diurnal discharge ranges. Therefore, inclusion of surface meltwater routing introduces smaller peak moulin discharge relative to the RCM instantaneous peak discharge (Figures 2b, 2d, and 2f). For IDC1, the peak moulin discharges simulated by the three routing models are >25 % lower than 20 the RCM instantaneous peak moulin discharge, while the diurnal discharge ranges are >27 % lower. Among the three routing models, RWF introduces slightly smoother moulin hydrographs than SUH and 2 m SRLF (Figures 2b, 2d, and 2f), yielding small differences in peak moulin discharge (<7%) and diurnal discharge range (<12%). This finding suggests that all three models are more representative of the observed processes than using RCM instantaneous runoff without routing, but that each model collects and distributes surface meltwater with different efficiency. 25

Simulations of subglacial effective pressure
When applied to provide meltwater input to the bed via a single moulin, the SRLF, SUH, and RWF meltwater routing models all introduce differences in diurnal and long-term fluctuations in subglacial effective pressure (ice pressurewater pressure) relative to effective pressure produced using RCM instantaneous runoff. The effective pressures presented are spatial mean effective pressures over the entire 1-km square domain of the bed surrounding the moulin input. The resultant 30 pattern of effective pressure simulation is similar among the four IDCs so we focus on IDC1 to illustrate our results here, with results for the other three IDCs presented in the supplement. For IDC1, the temporal mean of the spatial mean effective pressure is relatively consistent between routing models, varying by < 4% (ranging from 3.36 MPa to 3.48 MPa). RCM instantaneous runoff produces a slightly higher temporal mean effective pressure of 3.53 MPa. The amplitude of the fluctuations around the mean, however, differs substantially between different routing models (Figure 2 and Figures S1-S4 in the supplement).
The relative magnitudes of diurnal signals seen in the meltwater inputs are reflected in the subglacial response. The 5 effective pressure amplitude scales relative to the meltwater input amplitude, with the largest amplitude produced by the RCM instantaneous runoff inputs (Figures 3). While all surface meltwater routing models dampen these amplitudes as compared to RCM instantaneous runoff, the models that produce relatively high-amplitude moulin inputs correspond to high-amplitude effective pressure variations. In all cases a clear preferential pathway develops from the moulin location to the outflow at the left edge of the domain. This efficient channelized drainage pathway evolves and fluctuates with the inputs 10 through time and is characterized by higher gap height and effective pressure (i.e. lower hydraulic head and water pressure) than the surrounding bed perpendicular to flow (Figure 4, and supplementary animation).
For all surface meltwater routing models, the daily minimum moulin input generally corresponds to daily maximum effective pressure (i.e. minimum water pressure)., The minimum effective pressure (maximum water pressure) occurs several hours later as moulin inputs increase again into a subglacial system that has shut down due to ice creep in the absence of 15 high meltwater input required to maintain a larger gap height. A magnified example of this timing is seen in Figure 5, which presents the average two-day cycle of moulin discharge input using the three routing models overlaid with effective pressure in IDC1. All three routing models achieve minimum moulin discharges around 09:00-11:00 and minimum effective pressures around 17:00-19:00, yielding a time lag of 8-9 hours. In contrast, the RCM instantaneous runoff achieves minimum moulin discharge around 00:00 and minimum effective pressure around 10:00. The timing of effective pressure 20 produced using RCM instantaneous runoff is visibly different than with the routing methods. Interestingly, the timing of minimum effective pressure simulated by the RCM instantaneous runoff is very close (~ 1 hour) to that of maximum effective pressure simulated by the routing models.
In general, we find that while the moulin inputs generated by the different surface routing models produce variations in diurnal range of effective pressure, the overall channelization behaviour and temporal mean effective pressure over the 31-25 day simulations are relatively consistent between models. The subglacial model domain and duration were chosen to illustrate the impact of the chosen supraglacial routing model on local subglacial hydrology in the vicinity of a moulin input at the bed. As such, our results cannot necessarily be extrapolated to infer large-scale or seasonal evolution of the subglacial hydrologic system in response to different surface forcings; however, the results do provide insight into the potential diurnal sensitivity of the subglacial system to changes in supraglacial meltwater routing and the associated modification of the 30 discharge hydrograph. For example, the amplitude of diurnal effective pressure variation for a particular day may range from < 0.5 MPa to > 3.0 MPa, depending on the surface method used (Figure 3). While future work should pursue broader-scope simulations to more fully investigate the larger-and longer-scale effects of surface input variations on effective pressure and ice dynamics, our results suggest that in a fully coupled ice dynamics/subglacial hydrology model, different routing methods may not produce significantly different cumulative or time-averaged effects in effective pressure for simulation time scales longer than daily.

Temporal evolution of moulin discharge simulations
To illustrate the strong potential dynamism of supraglacial stream/river network extent on bare ice, Figure 6 presents dynamic supraglacial stream/river networks of IDC1 under different channel initial Ac thresholds (proxy for time). Numerous 5 small supraglacial streams are generated when Ac is set to 5000 m 2 (Figure 6a), while only large main stem channels are obtained when Ac is set to 10 6 m 2 (Figure 6d), yielding drainage density varying from 19.3 m -1 to 0.8 m -1 . Such strong, seasonal variation in supraglacial stream/river network extent is well-supported by previous high-resolution remote sensing studies, which report similar expansion and contraction of actively flowing supraglacial stream/river networks during the melt season (Smith et al., 2015;Smith et al., 2017;Yang et al., 2018). 10 Well-developed supraglacial stream/river networks (e.g., simulated using Ac = 100 m 2 ) route surface meltwater efficiently (Smith et al., 2015), yielding high UH and moulin discharge peaks, whereas relatively poorly-developed supraglacial stream/river networks indicate relatively slow routing of meltwater in expanding interfluve zones , yielding smaller and delayed UH and moulin discharge peaks (Figures 2g and 2h). As such, temporal evolution of supraglacial stream/river networks substantially alter surface meltwater routing processes and yield dynamic unit 15 hydrographs and moulin discharges. This result is consistent with Montgomery and Foufoula-Georgiou (1993) and van Meerveld et al. (2019), which found that river network extent controls the response of catchments to precipitation (melt in our case). During the first 10 days (240 hours), the channel initiation thresholds Ac are small so numerous supraglacial streams/rivers are developed; thereby, the dynamic Ac RWF routing performs similarly with SUH, 2 m SRLF, and the original RWF. As supraglacial stream/river networks shrink (i.e., Ac becomes larger), our simple modelling experiment (by 20 increasing Ac thresholds every 5 days for the month of July 2015) suggests a gradual dampening of diurnal variations until, by the end of the month, such variations are ~50 % smaller than the same routing method using a smaller Ac value ( Figure   2h). Moreover, the average two-day cycle of moulin discharge simulated by the dynamic Ac RWF routing is considerably dampened compared to those simulated by SUH, 2 m SRLF, and the original RWF ( Figure 6). In the dynamic Ac RWF routing, the dampening effect displayed in the moulin discharge inputs carries over into the resulting effective pressure 25 cycles as well (Figures 3 and 6).

Effects of DEM spatial resolution on surface meltwater routing
RWF routing is largely unaffected by DEM spatial resolution, but SRLF routing is significantly affected. Resampling 2 m ArcticDEM data to 5 m and 10 m yields similar RWF UHs and moulin discharges (Figures 2e and 2f), but progressively smoother versions using SRLF UHs (Figures 2c and 2d). At 90 m resolution, SRLF UH shapes exhibit diminished and 30 delayed peaks (Figure 2c). For IDC1, for example, use of a 90 m DEM yields a ~0.10 UH peak at hour 9 and a ~0.04 UH peak at hour 17, whereas 2 m DEM yield a ~0.15 UH peak at hour 5 with all meltwater evacuated within 15 hours. SRLF applied on higher-resolution DEMs yields larger peak moulin discharges and diurnal discharge signals (Figure 2d).
For IDC1, SRLF-routed peak moulin discharge derived from 2 m DEM is 52.4 % larger than that from 90 m DEM, while the corresponding value for SRLF-routed diurnal discharge range is 179.0 %. This finding indicates that SRLF performance is 5 strongly controlled by DEM spatial resolution, with coarser resolutions resulting increased damping of simulated moulin hydrographs. The dampening effect of SRLF moulin discharge induced by DEM spatial resolutions carries over into the resulting effective pressure cycles as well (Figure 3b), urging caution when using coarse-resolution DEMs to route surface meltwater and to simulate subglacial effective pressure.

Implications of surface meltwater routing method inter-comparison
This study conducts an inter-comparison of three surface meltwater routing models (SUH, SRLF, and RWF). Due to the lack of field-measured moulin hydrographs, it is difficult to determine which surface meltwater routing model performs most realistically over large areas and long times. However, our simulations of moulin discharges (Figures 2, and S1-S4) can provide some insightful information to qualitatively evaluate the performance of certain routing methods. First, inclusion of 15 surface meltwater routing introduces lower peak moulin discharges and delayed time to peak relative to the RCM instantaneous runoff. Second, for IDCs with similar areas and elevations as those examined here, simulated peak moulin discharge consistently occurs between 19:30 and 22:00 Smith et al., 2017). The timing of peak discharge suggests that several routing scenarios are unlikely to be realistic, including any technique that uses unmodified RCM outputs as these values qualitatively underestimate peak discharge time. 20 Additionally, we can use observations of channelized subglacial pressure to strengthen our qualitative comparison.
While exact field measurements are not available for direct comparison to our modelled effective pressures, limited field observations indicate that hydraulic head within subglacial channels varies diurnally by at least 40 m and up to 150 m in regions with slightly thinner ice Meierbachtol et al., 2013;Andrews et al., 2014). Furthermore, these reported pressure variations in channelized regions do not fall below 40% of overburden at the observed points. While 25 somewhat influenced by the particular geometry, moulin inputs, and other parameters used, our model results for effective pressure in the vicinity of a moulin show modelled mean effective pressure variations ranging from ~2 MPa to ~4 MPa for IDC1 (~26% to 50% of overburden pressure). The RCM runoff produces effective pressures with the largest amplitude fluctuations, and is unlikely to provide physically realistic ice surface runoff and the resolution of SRLF needs to be carefully chosen to generate realistic diurnal meltwater inputs to the subglacial drainage system. 30 Moreover, SUH and RWF both rely on several empirically parameters (Cp and Ct for SUH, and vh and vc for RWF) calibrated from a moulin hydrograph measured at the Rio Behar catchment, southwestern GrIS during a very short time period (72 hours), July 2015 (Smith et al., 2017). SRLF is applicable over large spaces and long times because it only relies on DEM to calculate meltwater flow velocities (Banwell et al., 2012). In this study, we assume these empirically parameters are transferable over space and time but this assumption needs further validation. It may hold for ice sheet surface with 5 similar hydrologic and glaciological environments but is problematic to apply over large spaces and long times due to evolving ice surface characteristics. Multiple independent, long-term moulin hydrographs will help eliminate the need for this assumption.

Influence on diurnal subglacial pressure variations
Our results demonstrate that the routing models and DEM resolution can modulate the diurnal variability of both 10 surface meltwater transport and associated moulin inputs (Figures 3 and 4). Surface meltwater routing alters hourly inputs, with only slight changes to total integrated daily moulin input. However, hourly variations in moulin discharge can influence ice dynamics. Diurnal variations in subglacial pressure, for example, are directly associated with diurnal variations in ice velocity, which have the capacity to induce additional surface-to-bed connections (Carmichael et al., 2015;Hoffman et al., 2018). The diurnal amplitude of moulin discharge also affects the transient evolution of the subglacial drainage system 15 without sustained input to maintain efficient drainage channels, effective pressure decreases, resulting in increased sliding velocities, potentially influencing seasonal ice displacement (Hewitt, 2013; van de Wal et al., 2015). As such, including the diurnal signal of moulin discharge may be important for accurately modelling surface-to-bed meltwater connections ( Figure   3 & S1-S3). As we see reflected in both the magnitude of effective pressure diurnal amplitudes ( Figure 3) and in the timing of minimum/maximum effective pressure relative to moulin inputs ( Figure 5), different surface routing methods yield a 20 range of diurnal behaviour. However, as described in Section 3.2, the temporal mean effective pressure produced in our subglacial hydrology simulations with inputs from the various routing models is relatively consistent between routing models. Depending on the time scale of interest, diurnal inputs may not be vital to capture the relevant ice dynamics.
Many subglacial hydrology models commonly invoke a numerical term (the "englacial void ratio") to represent englacial storage in order to provide short term storage and release of meltwater that cannot be accommodated rapidly within 25 the subglacial system, in the absence of more realistic representation of supraglacial and englacial storage (Hewitt, 2013;Werder et al., 2013;Hoffman et al., 2016). In SHAKTI, this englacial void ratio is also included as an option (Sommers et al., 2018), but our simulations considered here do not employ this term in the equations. Our results demonstrate that the supraglacial hydrologic system acts as short-term storage for surface-derived meltwater, as exhibited by the time lag of moulin inputs between models; therefore, application of an appropriate surface meltwater routing scheme may reduce the 30 dependence of some subglacial models on a somewhat arbitrary englacial storage term to produce realistic diurnal effective pressure variations and timing lags (Werder et al., 2013;Hoffman et al., 2016). The assumption of surface inputs being instantaneously delivered to the subglacial system is an approximation that ignores the largely uncertain complex flow paths through the englacial system from surface to bed, but it is reasonable to assume that variability in surface inputs should influence variability in inputs to the subglacial system. Supraglacial routing delays can affect the amplitude and timing of subglacial pressure variations (Figures 3 and 6).
While there is limited evidence that within a fully coupled ice dynamics/subglacial hydrology model that diurnal variations 5 in effective pressure can contribute to changes in ice dynamics (Hewitt, 2013), the evolution of the supraglacial hydrologic system can substantially alter effective pressure on diurnal timescales (e.g. Figure 3d). Subglacial channel development is key to terminating early melt season accelerated sliding (Hoffman and Price, 2014) and subglacial channelization occurs more readily and persistently under constant supraglacial meltwater inputs (Schoof, 2010;Poinar et al., 2019). The appropriate choice of a time evolving surface meltwater routing scheme may, in addition to providing the best representation 10 of diurnal effective pressure variations, result in an improved representation of subglacial evolution and improve quantitative representation of seasonal variations in ice flow in coupled models. Our results, however, show that while meltwater inputs with highly variable diurnal amplitudes yield effective pressures with highly variable diurnal amplitudes, the general channelization behaviour and mean effective pressure over the 31-day simulations are quite similar across routing models.
Further work is needed to explore the implications of using sub-daily inputs versus time-averaged inputs on long-term 15 subglacial hydrology and ice dynamics in realistic, large domains. The case may be that the peaks and troughs of largeamplitude diurnal input variations effectively balance out to yield equivalent results as time-averaged inputs.

Influence of seasonal supraglacial drainage evolution on meltwater routing
Seasonal evolution of ice surface drainage pattern can substantially alter surface meltwater routing and moulin discharge characteristics (Figures 2 and S1-S4). Such changes also involve the removal and deposition of a winter snowpack 20 (Nienow et al., 2017) and the development of a weathering crust on base ice , both of which can reduce the prevalence of open-channel flow . While these processes cannot be explicitly represented with the meltwater routing parameterizations described here, time evolving UHs can integrate the impact of seasonal changes in the relative proportions of porous and open-channel flow. As such, the relative ease of creating dynamic UHs to mimic seasonal evolution of supraglacial stream/river networks is a distinctive advantage of RWF (Figures 2g and 2h). 25 The development of time varying UHs provides an opportunity to simulate dynamic moulin discharges (Smith et al., 2017;Yang et al., 2018). Such time-varying UHs could be important to realistically modelling the evolution of the subglacial drainage system and the development of subglacial channelsa limited supraglacial stream/river extent will dampen the diurnal range of meltwater inputs into the subglacial system and maintain a more constant subglacial pressure (e.g. Figure   3d), which in turn results in less variation in ice motion. Furthermore, limiting diurnal variations may result in more rapid 30 growth of subglacial channels early in the melt season (Schoof, 2010;Hewitt, 2013), thus, potentially produce a more realistic transition to primarily channelized drainage (Banwell et al., 2013;Banwell et al., 2016).

Impact of DEM resolution on supraglacial meltwater routing
Analysing the influence of DEM spatial resolution on the hydrological response of catchment to precipitation is an important research topic in terrestrial hydrology (Zhang and Montgomery, 1994;Hancock et al., 2006). We expand it to ice sheet hydrology and attempt to investigate the impact of DEM spatial resolution on the hydrological response of supraglacial catchment to surface melt. We found that RWF routing is not affected by DEM spatial resolution, whereas SRLF routing is 5 (Table 2 and Figure 2). The 5 m SRLF UHs and resultant moulin discharges best match with RWF and SUH simulations (Figures 2). This is consistent with the finding that "the most appropriate DEM grid size for topographically driven hydrologic models is somewhat finer than the interfluve scale identifiable in the field" (Zhang and Montgomery, 1994). This behaviour occurs because a coarser-resolution DEM represents bare ice surface topography more smoothly and yields lower bare ice surface slopes, similar to terrestrial topography (Zhang and Montgomery, 1994); lower topographical slopes yield 10 smaller velocities via Manning's open-channel equation.
Different meltwater routing velocities control the meltwater transport time because meltwater flow length is minimally impacted by DEM resolutions. As a result, lower SRLF meltwater routing velocities induce longer meltwater transport time and distribute diurnal surface runoff more smoothly over time (Figures 2c and 2d), whereas stable flow length contributes to RWF's better performances under different DEM resolutions (Figures 2e and 2f). 15

Future research directions of surface-to-bed meltwater connection
A challenge remaining unsolved is to map or simulate seasonal evolution of supraglacial stream/river networks. Spatial extent of supraglacial stream/river networks determine partition of interfluve and open-channel zones and thereby controls IDC hydrologic response to surface melt. Moreover, all three routing models described in this study do not explicitly include complex ice surface composition and transformation, including the retreat of seasonal snow cover (Hubbard and Nienow, 20 1997) and the modification of the ice surface by incoming solar radiation and other processes in the interfluve zone (Karlstrom et al., 2014;Cooper et al., 2018). This may lead to uncertainties in unit hydrograph and moulin discharge simulations. Recently, weathering crust has been found to be widely distributed on the Greenland ice surface , rather than impermeable bare ice layer as previously assumed (Arnold et al., 1998). Therefore, an appropriate parameterization of porous media flow may be important to accurately describe meltwater transport in the interfluve zone of 25 ice surface (Karlstrom et al., 2014;Cooper et al., 2018;Yang et al., 2018). Additionally, the hourly MAR runoff is a cumulative runoff on one hour and may induce a potential delay of one hour in the peaks; therefore, 10-15 minute runoff outputs could be used to better capture runoff delays in future.
Furthermore, the subglacial hydrology simulations in this study should be interpreted as a focused view into how different meltwater routing methods influence local subglacial hydrology in the vicinity of moulins through variations in 30 moulin inputs. Regional simulations covering multiple IDCs in a fully coupled ice dynamics-subglacial hydrology model remain for future work, an important next step in understanding the surface-to-bed meltwater influence on regional and large-scale ice dynamics. A more direct comparison of model results with observations will be appropriate for this future work, in which the entire drainage catchments are modelled with multiple moulin distributed in areas where borehole observations are available or relevant quantities may be inferred from radar products (Chu et al., 2016).

Conclusions
Surface meltwater routing is crucial for understanding the surface-to-bed meltwater connection of the Greenland Ice Sheet 5 but remains poorly studied to date. This study presents a first inter-comparison of three different meltwater routing models and employs a subglacial hydrologic model to explore the impacts of their differences on subglacial effective pressure (ice pressurewater pressure) in the vicinity of moulins. Results show that inclusion of surface meltwater routing introduces significantly small peak moulin discharges and delayed time to peak relative to the RCM instantaneous runoff. Different surface meltwater routing models, as well as different spatial-resolution DEMs and seasonal evolution of the supraglacial 10 stream/river networks, induce variable diurnal moulin discharges and effective pressures, which influence ice sliding velocity. While the different routing models produce different diurnal amplitudes in effective pressure variation, the overall channelization behaviour and temporal mean effective pressures are relatively consistent across surface meltwater routing models. Together, these findings urge caution for better representations of surface meltwater routing and moulin discharge simulations to drive subglacial hydrology and ice dynamics models, as well as highlight the need for further research to 15 investigate the cumulative effects of diurnal inputs to the subglacial drainage system and the relevant impacts on ice dynamics.

Author contributions 25
KY and AS designed the study. KY and AS performed the data analysis. KY wrote the paper with contributions from all authors.   Figure 2) generated by the various surface meltwater routing models. The effective pressure shown is the spatial mean over the entire 1-km square domain which contains the moulin input at its center.   2 and 3. The daily minimum input in supraglacial moulin discharge (solid lines) corresponds generally to maximum effective pressure (dashed lines), and is followed within 8-9 hours by the daily minimum effective pressure (maximum subglacial water pressure). This suggests that the system shuts down due to creep with low meltwater input, and becomes highly pressurized as meltwater input increases again.