Intercomparison 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 (GrIS) 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 pressures. We compare hourly supraglacial moulin discharge simulations from three surface meltwater routing models – the synthetic unit hydrograph (SUH), the bare-ice component of surface routing and lake filling (SRLF), and the 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 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 produce significant differences in peak moulin discharge and time to peak. In particular, the RWF yields later, smaller peak moulin discharges than the SUH or SRLF due to its representation of slow interfluve flow between supraglacial meltwater channels, and it can readily accommodate the seasonal evolution of supraglacial stream and river networks. Differences among the three routing models are reflected in a 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.

1 Introduction

Large volumes of meltwater are routed through supraglacial stream and river networks on the Greenland ice sheet (GrIS) 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 (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 uncertainties in hydrodynamic coupling with ice motion.

In most studies that investigate surface-to-bed meltwater connections and the behavior of the subglacial hydrologic system, surface meltwater routing is either simplified or simply ignored, i.e., moulin hydrographs are estimated directly from regional climate model (RCM) output. In these instances, internally drained catchments (IDCs) drain surface meltwater through supraglacial stream and river networks to terminal moulins or lakes (Yang and Smith, 2016). Then, RCM instantaneous grid cell runoff is simply summed over each IDC to estimate its corresponding supraglacial moulin 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 longer-term studies; however, surface meltwater routing has been found to substantially modify ice surface runoff (Banwell et al., 2012, 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 subglacial hydrologic system evolution and hydrodynamics, especially for large IDCs and short (i.e., diurnal) timescales. 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 (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 (Yang et al., 2018). Alternatively, empirical unit hydrograph (UH) methods can be used to estimate the hydrologic response based on catchment shape and area alone. The UH is “a transfer function that is widely used for modeling catchment runoff response to rainfall events for some unit duration and unit depth of effective water input” (Smith et al., 2017) and thus enables 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 digital elevation model (DEM)-based surface routing and lake filling (SRLF) model. The 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 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 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.

Most recently, Yang et al. (2018) utilized the 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 ($\sim {\mathrm{10}}^{-\mathrm{3}}$–10⁻⁴ m s⁻¹) and fast open-channel flow ($\sim {\mathrm{10}}^{-\mathrm{1}}$ m s⁻¹) combined to simulate the downstream moulin hydrograph.

These three surface meltwater routing models (SRLF, SUH, and RWF) use different assumptions and have different data requirements (Table 1). The 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, the RWF requires a partitioning of slow (interfluve) and fast (open-channel) flow that can substantially impact the performance of the RWF without careful spatial and temporal calibration. Moreover, the SRLF and RWF rely on surface DEMs to calculate meltwater flow paths (Yang et al., 2018), and the 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 high-resolution DEM datasets (e.g., 2 m ArcticDEM) of the GrIS surface (Noh and Howat, 2015, 2017, 2018).

Figure 1 Four supraglacial internally drained catchments (IDCs) on the southwestern Greenland ice sheet were selected for intercomparison of surface meltwater routing models. IDC topographic boundaries (black lines) are extracted from high-resolution (2 m) ArcticDEMs. Supraglacial river networks and lakes (blue lines) are mapped from a Landsat-8 panchromatic image (Yang and Smith, 2016). The base map is 10 m Sentinel-2 image acquired on 30 July 2018.

Table 1 A brief summarization of surface meltwater routing models.

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. The routing models' ability to simulate meltwater transport velocity, flow length, transport time, UH, peak moulin discharge, and time to peak is compared, and the impact of DEM spatial resolution on SRLF and RWF is examined. The differences in using these models to drive subglacial effective-pressure variations are qualitatively investigated with an idealized case of surface 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).

Notably, this study can, in good faith, neither focus on method comparison nor determine the “best” (i.e., best able to reproduce a real-world moulin hydrograph) model 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. We conclude with a general discussion of surface meltwater routing on the GrIS and include some recommendations for best practices and future research.

2 Methods and data

2.1 Study area and data sources

We selected four moderate-size supraglacial IDCs in the Russell Glacier region, southwestern GrIS, to explore the impact of the SUH, SRLF, and RWF meltwater routing models on moulin discharge and subglacial effective pressure (Fig. 1, Table 2). They are distributed at approximately 200 m elevation intervals in order to span the elevational range of most well-developed IDCs found in the Russell Glacier region and the variable surface melt conditions of this region (Fitzpatrick et al., 2014; Yang and Smith, 2016). Large supraglacial lakes can store considerable volumes of surface meltwater and thus impact moulin discharge at the IDC outlet (Yang et al., 2019). Therefore, these four IDCs were deliberately chosen to avoid the presence of large supraglacial lakes, and all surface meltwater is routed to the terminal moulin at the catchment outlet (Fig. 1). As such, surface runoff produced in each IDC should be equal to the moulin discharge (Smith et al., 2017). The four IDC areas range from 53.0 to 66.9 km²; surface elevations vary from 1086 to 1672 m and ice thicknesses from 854 to 1432 m (Table 2). A moulin discharge hydrograph collected at Rio Behar catchment (IDC2 in this study; 67.049346^∘ N, 49.025809^∘ W) for 72 h, from 20 to 23 July 2015, was used to calibrate parameters of the SUH and RWF models (see Sect. 2.5 and 2.6). It is problematic to apply these empirically derived parameters over large spaces and long times (Yang et al., 2018). 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).

Table 2 Summary of four study catchments.

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 the Polar Geospatial Center (PGC) and provide an unprecedented opportunity to investigate ice surface landscape and hydrologic processes (Noh and Howat, 2015, 2017, 2018). Supraglacial IDC boundaries were derived from ArcticDEM following the method of Karlstrom and Yang (2016). 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. Catchment-averaged hourly runoff (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 runoff is termed RCM instantaneous runoff.

2.2 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 contrast, a routing model can be lumped and only determines how surface meltwater produced in the catchment is temporally routed to the catchment outlet. The 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):

$$\begin{array}{}\text{(1)}& Q=R*\mathrm{UH},\end{array}$$

where * is the convolution operator. For each IDC, convolution of the UH with hourly R yields a discharge hydrograph at the 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 the UH can be considered routing models.

2.3 Snyder synthetic unit hydrograph

The Snyder SUH (Snyder, 1938) assumes that catchment morphometry controls the shape of the UH and utilizes the catchment main-stem stream length (L, in kilometers), the distance from the catchment outlet (here, the terminal moulin) to the point on the main channel nearest to the catchment centroid (L_ca, in kilometers), and two parameters, C_p and C_t, to determine peak discharge (h_p, per hour) and time to peak (t_p, in hours) of the SUH as

$$\begin{array}{}\text{(2)}& {\displaystyle }{t}_{\mathrm{p}}=C_{\mathrm{t}}(L{L}_{\mathrm{ca}}{)}^{\mathrm{0.3}},\text{(3)}& {\displaystyle }{h}_{\mathrm{p}}=C_{\mathrm{p}}/t_{\mathrm{p}},\end{array}$$

C_p and C_t 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–23 July 2015 to calibrate C_p and C_t. Following Smith et al. (2017), we use C_p=0.72 and C_t=1.61 to determine h_p and t_p 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, C_p and C_t, 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 for the Russell Glacier region, southwestern GrIS, during July 2015.

2.4 Surface routing and lake filling

The SRLF has been used to model filling of 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 flow over bare-ice surfaces (Arnold et al., 1998). The 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, 2013, 2016; Arnold et al., 2014). 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. The SRLF assumes universal presence of supraglacial open-channel flow everywhere on the bare-ice surface by applying Manning's open-channel flow equation (Manning, 1891) to calculate channelized meltwater velocity v for every pixel along a DEM-derived flow path:

$$\begin{array}{}\text{(4)}& v=R_{\mathrm{H}}^{\mathrm{2}/\mathrm{3}}{S}^{\mathrm{1}/\mathrm{2}}/n,\end{array}$$

where R_H is the hydraulic radius of the supraglacial meltwater channel, S is water surface slope calculated from the DEM, and n is the Manning roughness parameter. Following Arnold et al. (1998), we set R_H=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 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 the SUH and RWF, a 1 h UH was calculated by hourly binning the transport time raster (Yang et al., 2018), with the resultant UH termed the SRLF UH.

2.5 Rescaled width function

The RWF routing model partitions surface meltwater transport pathways into interfluve (i.e., hillslope) and open-channel distances (L_h and L_c) to determine two catchment-averaged meltwater velocities (v_h and v_c) to represent routing efficiency over these regions (D'Odorico and Rigon, 2003). Then, the transport time for each pixel in a catchment can be calculated as

$$\begin{array}{}\text{(5)}& t=t_{\mathrm{h}}+t_{\mathrm{c}}=L_{\mathrm{h}}/v_{\mathrm{h}}+L_{\mathrm{c}}/v_{\mathrm{c}},\end{array}$$

where t_h is the interfluve travel time, and t_c is the channel travel time. Similar to the 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 (Yang et al., 2018). By applying the constant interfluve velocity v_h to interfluve pixels and channel velocity v_c 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 the RWF UH. The primary difference between the RWF and SRLF is the way they determine meltwater flow velocity for each catchment cell. For the RWF, using the in situ moulin hydrograph of Smith et al. (2017), Yang et al. (2018) determined the velocity combination of v_h=0.0006 and v_c=0.4 m s⁻¹ to be the optimal match to field observations. Although v_h and v_c may change notably during the melt season because of variable surface melt intensity and ice surface topography, it is reasonable to assume that v_h and v_c are persistent during a relatively short time period for similar hydrological and glaciological environments. Therefore, in this study, the optimal velocity combination obtained in Yang et al. (2018) is assumed to be transferable and used to create RWF UHs for the four IDCs.

2.6 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, 10, 30, and 90 m resolutions, similar to Zhang and Montgomery (1994), and the corresponding meltwater flow paths, velocities, transport time, UHs, and moulin hydrographs for the SRLF and RWF subsequently calculated. For the SRLF, all five resampled DEMs were used. For the RWF, only 2, 5, and 10 m resolutions were used to prevent DEM resolution from exceeding the typical interfluve transport distance, ∼10 m for southwestern GrIS (Yang et al., 2018); otherwise, interfluve transport distance would be overestimated, and the resultant hydrograph would be inappropriate (Hancock et al., 2006).

2.7 Temporal evolution of surface meltwater routing

The spatial pattern of supraglacial stream and 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 open-channel and interfluve transport distances change (Montgomery and Foufoula-Georgiou, 1993; van Meerveld et al., 2019). The RWF can mimic this effect by varying the partitioning of interfluve versus open-channel zones (Yang et al., 2018). The bare-ice component of the 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. The SUH relies on C_p and C_t to build the UH. If a moulin hydrograph is available during the entire melt season, we can calculate time-evolving C_p and C_t using variable time to peak (t_p), peak discharge (h_p), and main-stem channel length (L), thus creating multitemporal UHs. Without such direct measurements, the parameters cannot be realistically varied in time, and the SUH cannot mimic variable hydrologic response of a catchment to surface melt.

Capturing the rapid seasonal evolution of supraglacial stream and river networks is challenging (Lampkin and VanderBerg, 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 initiation threshold, A_c) values (A_c=100, 250, 500, 1000, 2500, 5000 m²) to simulate the hypothetical seasonal evolution of supraglacial stream and river networks by evolving the partitioning of interfluve vs. open channels and resultant surface meltwater routing simulations (Yang et al., 2018). A_c defines the surface area needed to initiate open-channel flow. Smaller values of A_c signify periods of intense surface melt and actively flowing, well-developed supraglacial stream and river networks (Smith et al., 2017), whereas larger values signify poorly developed network characteristics of the late melt season (Yang et al., 2018). A total of six 5 d RWF UHs were created for each A_c value, broadly consistent with a decreasing RCM runoff trend during July 2015. Each RWF UH was used to calculate moulin discharge for 5 d, and the resultant moulin discharge is termed as dynamic A_c discharge.

2.8 Subglacial hydrology modeling

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 modeled effective pressure resulting from the different surface meltwater routing models considered in this study, we treat their respective moulin discharges as direct meltwater inputs to the bed via a single moulin and simulate subsequent evolution of the local subglacial drainage system in an idealized setting using the SHAKTI model (Sommers et al., 2018).

SHAKTI is built within the Ice-sheet and Sea-level System Model (ISSM; 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 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).

To explore and compare the influence of different surface meltwater routing models on idealized subglacial hydrology in a small experiment, we apply surface 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 Neumann boundary conditions on the other three sides (i.e., enforcing $\mathrm{d}h/\mathrm{d}x=\mathrm{0}$ at the upstream right boundary and $\mathrm{d}h/\mathrm{d}y=\mathrm{0}$ at the top and bottom boundaries, where h is hydraulic head, which is directly related to water pressure). For each IDC, the mean surface slope and mean ice thickness are used (Morlighem et al., 2017) for the domain geometry, with the bed slope set equal to the surface slope to create a uniformly thick tilted slab of ice, and constant sliding velocity is prescribed as 100 % of the annual mean observed surface velocity (Table 3). Note that two-way coupling is not currently enabled between SHAKTI and the ISSM; the prescribed sliding velocity influences the subglacial system, but the feedback of evolving effective pressure on sliding velocity is not implemented. Constants and parameter values used in the simulations are summarized in Table 3. Certain parameters are highly uncertain (for example, bed bump height and bed bump spacing), so we select values used in earlier work with SHAKTI (Sommers et al., 2018) and in the Subglacial Hydrology Model Intercomparison Project (de Fleurian et al., 2018). There is a great need to better constrain uncertain parameters used in subglacial hydrology models, but that endeavor is beyond the scope of this work.

Table 3 Constants and parameters used in subglacial hydrology simulations.

3 Results

3.1 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 and 15:00 LT (local time) for all the four studied IDCs (Table 2). Inclusion of routing processes, however, introduces considerable delay between the 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 LT (Table 2, Figs. 2 and S1–S3 in the Supplement).

Figure 2 Presentation of unit hydrographs (UHs; column 1) and moulin discharges (column 2) of IDC1 during July 2015, as simulated by three surface meltwater routing models (SUH, RWF, and SRLF). Simultaneous RCM runoff (gray line) is shown to indicate the effect of surface meltwater routing processes on moulin discharge.

SUH-routed hydrographs simulate peak moulin discharge between approximately 19:00 and 21:00 LT, a delay of 6 h compared to the RCM instantaneous peak discharge, while SRLF-routed delays vary from 3 to 10 h, and RWF-routed delay varies from 7 to 10 h (Table 2). These delays are longer than the ∼2 h measured for small IDCs near ice margins (Shepherd et al., 2009; McGrath et al., 2011; Bartholomew et al., 2012) but are comparable to the ∼5–7 h measured for the ∼60 km² Rio Behar catchment, consistent with findings indicating that catchment size partially controls 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 UHs for each routing model. The UH isolates the impact of basin characteristics on moulin discharge by providing the basin-integrated runoff resulting from 1 cm of surface melt over 1 h. The resultant patterns of moulin hydrographs for the four IDCs are similar, so we use IDC1 to illustrate our results here. The three UHs (i.e., SUH, SRLF UH, and RWF UH) and their routed moulin discharges for the other three IDCs are presented in the Supplement (Figs. S1–S3). For IDC1, the peak values of RWF UHs are lower than SUHs and 2 m SRLF UHs. Therefore, RWF UHs temporally distribute surface meltwater most smoothly (Fig. 2a, c, and e). The peak RWF UH is ∼0.10 (Fig. 2e), indicating that ∼10 % surface meltwater is routed to the moulin during the peak runoff hour (Singh et al., 2014), while the peak values are ∼0.13 for the SUH (Fig. 2a) and ∼0.15 for the 2 m SRLF UH (Fig. 2c).

UHs with a lower amplitude 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 (Figs. 2b, d, and f). For IDC1, the peak moulin discharges simulated by the three routing models are >25 % lower than the RCM instantaneous peak moulin discharge, while the diurnal discharge ranges are >27 % lower. Among the three routing models, the RWF introduces slightly smoother moulin hydrographs than the SUH and 2 m SRLF (Fig. 2b, d, and f), 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.

3.2 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 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 temporal pattern of effective-pressure simulation is similar among the four IDCs, although the magnitude and amplitude differ substantially due to differences in ice thickness, slope, and sliding velocity of each catchment. We focus on IDC1 to illustrate our results here (Fig. 3), with results for the other three IDCs presented in the Supplement (Figs. S4–S6). For IDC1, the temporal mean of the spatial mean effective pressure is relatively consistent between routing models, varying by <4 % (ranging from 3.36 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 3 Subglacial effective pressures for IDC1 as simulated by SHAKTI, with inputs to the subglacial system via a single moulin prescribed by the moulin discharges (shown in Fig. 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.

The relative magnitudes of diurnal variations seen in the meltwater inputs are reflected in the subglacial response. The effective-pressure amplitude scales relative to the meltwater input amplitude, with the largest amplitude produced by the RCM instantaneous runoff inputs (Fig. 3). While all surface meltwater routing models damp 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 (i.e., a low “channel”) develops from the moulin location to the outflow at the left edge of the domain, as expected, based on the topography and boundary conditions of the experiment design (Fig. 4). This efficient drainage pathway evolves with the inputs through time and is characterized by higher gap height and effective pressure (lower hydraulic head and water pressure) than the surrounding bed perpendicular to flow (Fig. 4, and the animation in Fig. S7 in the Supplement).

Figure 4 Snapshots of dynamic subglacial hydrology fields on day 23 in IDC1 using the SUH routing method to drive moulin input (see full animation of channel evolution and fluctuation in the Supplement). An efficient channelized drainage pathway develops from the moulin location at the center of the domain to the outflow at the left, characterized by higher gap height, water flux, effective pressure, and lower hydraulic head than its surroundings perpendicular to flow.

For all surface meltwater routing models, the daily minimum moulin input generally corresponds to daily maximum effective pressure (i.e., minimum water pressure) in these experiments. 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 high meltwater input required to maintain a larger gap height. A magnified example of this timing is seen in Fig. 5, which presents the average 2 d cycle of moulin discharge input using the three routing models overlaid with effective pressure in IDC1 (results for the other three IDCs are presented in Figs. S8–S10). All three routing models achieve minimum moulin discharges around 09:00–11:00 LT and minimum effective pressures around 17:00–19:00 LT, yielding a time lag of 8–9 h. In contrast, the RCM instantaneous runoff achieves minimum moulin discharge around 00:00 LT and minimum effective pressure around 10:00 LT. The timing of effective pressure 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 to that of maximum effective pressure simulated by the routing models (∼1 h difference).

Figure 5 The average 2 d cycle of moulin discharge (Q) for IDC1 during July 2015 derived from Figs. 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 h 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. As the new water inputs are accommodated, efficient pathways reform, and effective pressure increases (subglacial water pressure decreases).

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 qualitative channelization behavior and temporal mean effective pressure over the 31 d simulations are relatively consistent between models (Figs. 3 and 5).

Note that the small, idealized subglacial model domain was designed to compare the impact of each surface routing model on local subglacial hydrology behavior in the vicinity of a moulin input at the bed. These numerical experiments explore differences resulting from forcing by the different surface routing models under the same conditions and are explicitly not a representation of realistic subglacial evolution in a given catchment. We acknowledge that uncertain parameters and the particular boundary conditions certainly play a role in the specific behavior and pressure magnitudes shown here, and these results are intended solely to illustrate and compare differences between methods in an idealized setting. To investigate the influence of the upstream boundary condition in our domain, we have compared sensitivity experiments (not shown) using a steady moulin input of 15 m³ s⁻¹ (the mean input of the SUH routing model), run to a steady subglacial drainage configuration, using a no-flux ($\mathrm{d}h/\mathrm{d}x=\mathrm{0}$) Neumann boundary as used in the experiments presented here vs. a Dirichlet boundary, setting the head equivalent to 50 % of the ice overburden pressure (identical to the downstream boundary condition). The spatial mean effective pressure (our quantity of interest) in the case with Dirichlet conditions at both upstream and downstream edges is 5 % higher than with the Neumann upstream condition. The prescribed pressure of the downstream condition also influences the magnitude of effective pressure: setting this to 40 % of overburden results in a mean effective pressure 17 % higher than with the 50 % overburden condition; increasing the downstream boundary to 60 % overburden decreases the resulting mean effective pressure by 15 %. We find that imposing high-pressure boundary conditions of 90 % overburden at the upstream and downstream edges appears to somewhat inhibit the development of the subglacial system in response to the moulin inputs as the entire domain is highly pressurized (and is therefore not as helpful for our purpose of comparing the behavior with inputs derived from the different routing models).

As such, our results should not 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 surface meltwater routing and the associated modification of the 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 in this experimental setup, depending on the surface routing model used (Fig. 3). We hope this will serve as inspiration for future work to pursue broader-scope simulations to thoroughly investigate the larger- and longer-scale effects of surface meltwater input variations on effective pressure and ice dynamics.

3.3 Temporal evolution of moulin discharge simulations

To illustrate the strong potential dynamism of supraglacial stream and river network extent on bare ice, Fig. 6 presents dynamic supraglacial stream and river networks of IDC1 under different channel initiation A_c thresholds (proxy for time). Numerous small supraglacial streams are generated when A_c is set to 5000 m² (Fig. 6a), while only large main-stem channels are obtained when A_c is set to 10⁶ m² (Fig. 6d), yielding drainage density varying from 19.3 to 0.8 km⁻¹. Such strong, seasonal variation in supraglacial stream and river network extent is well supported by previous high-resolution remote-sensing studies, which report similar expansion and contraction of actively flowing supraglacial stream and river networks during the melt season (Smith et al., 2015, 2017; Yang et al., 2018).

Figure 6 Variable supraglacial stream and river network for IDC1, as simulated by applying variable accumulative area threshold (channel initiation threshold, A_c) values to ArcticDEM.

Well-developed supraglacial stream and river networks (e.g., simulated using A_c=100 m²) route surface meltwater efficiently (Smith et al., 2015), yielding high UH and moulin discharge peaks, whereas relatively poorly developed supraglacial stream and river networks indicate relatively slow routing of meltwater in expanding interfluve zones (Yang et al., 2018), yielding smaller and delayed UH and moulin discharge peaks (Fig. 7a and b). As such, temporal evolution of supraglacial stream and river networks substantially alter surface meltwater routing processes and yield dynamic unit hydrographs and moulin discharges. This result is consistent with Montgomery and Foufoula-Georgiou (1993) and van Meerveld et al. (2019), who found that river network extent controls the response of catchments to precipitation (melt in our case). During the first 10 d (240 h), the channel initiation thresholds A_c are small, so numerous supraglacial streams/rivers are developed; thereby, the dynamic A_c RWF routing performs similarly with the SUH, 2 m SRLF, and the original RWF. As supraglacial stream and river networks shrink (i.e., A_c becomes larger), our simple modeling experiment (by increasing A_c thresholds every 5 d for the month of July 2015) suggests a gradual damping of diurnal variations until, by the end of the month, such variations are ∼50 % smaller than the same routing method using a smaller A_c value (Fig. 7c). Moreover, in the dynamic A_c RWF routing, the damping effect displayed in the moulin discharge inputs carries over into the resulting effective-pressure cycles as well (Fig. 7c).

Figure 7 Presentation of (a) unit hydrographs (UHs), (b) moulin discharges, and (c) subglacial effective pressure of IDC1 during July 2015, as simulated by dynamic A_c values. A_c is the cumulative contributing area required to initiate a supraglacial meltwater channel, and dynamic A_c values are used as a proxy for time to simulate the temporal evolution of supraglacial stream and river networks. Simultaneous RCM runoff (gray line) is shown to indicate the effect of surface meltwater routing processes on moulin discharge.

3.4 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 (Fig. 2e and f) but progressively lower-amplitude SRLF UHs (Fig. 2c and d). The RWF relies on high-resolution DEMs to calculate interfluve transport distance and cannot use DEMs of coarser resolution (>10 m; Yang et al., 2018). At 90 m resolution, SRLF UH shapes exhibit diminished and delayed peaks (Fig. 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 DEMs yield a ∼0.15 UH peak at hour 5, with all meltwater evacuated within 15 h. The SRLF applied on higher-resolution DEMs yields larger peak moulin discharges and diurnal discharge signals (Fig. 2d). For IDC1, SRLF-routed peak moulin discharge derived from 2 m DEMs is 52.4 % larger than that from 90 m DEMs, while the corresponding value for SRLF-routed diurnal discharge range is 179.0 %. This finding indicates that SRLF performance is strongly controlled by DEM spatial resolution, with coarser resolutions resulting in increased damping of simulated moulin hydrographs. The damping effect of SRLF moulin discharge induced by DEM spatial resolutions carries over into the resulting effective-pressure cycles as well (Fig. 3b), urging caution when using coarse-resolution DEMs to route surface meltwater and to simulate subglacial effective pressure.

4 Discussion

4.1 Intercomparison of surface meltwater routing models

This study conducts an intercomparison 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 (Figs. 2 and S1–S3) can provide some insightful information to qualitatively evaluate the performance of certain routing methods. First, inclusion of 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 LT (Chandler et al., 2013; 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.

Moreover, the SUH and RWF both rely on several empirical parameters (C_p and C_t for the SUH and v_h and v_c for the RWF) calibrated from a moulin hydrograph measured at the Rio Behar catchment, southwestern GrIS, during a very short time period (72 h), in July 2015 (Smith et al., 2017). The SRLF is applicable over large spaces and long times because it only relies on DEMs to calculate meltwater flow velocities (Banwell et al., 2012). In this study, we assume these empirical parameters are transferable over space and time, but this assumption needs further validation. It may hold for ice sheet surface with 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.

4.2 Influence on diurnal subglacial pressure variations

Our results demonstrate that the routing models and DEM resolution can modulate the diurnal variability of both surface meltwater transport and associated moulin inputs (Fig. 2). Surface meltwater routing alters hourly inputs, with only slight changes to total integrated daily inputs. However, hourly variations in moulin discharge can influence ice dynamics (Fig. 3). Diurnal variations in subglacial pressure, for example, can trigger variations in ice velocity and 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 without sustained input to maintain efficient drainage channels, and 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 modeling surface-to-bed meltwater connections (Figs. 3–5 and S4–S10). As we see reflected in both the magnitude of effective-pressure diurnal amplitudes (Fig. 3) and in the timing of minimum and maximum effective pressure relative to moulin inputs (Fig. 5), different surface routing models yield a range of diurnal behavior. However, as described in Sect. 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 even though the amplitude varies greatly. This suggests that depending on the timescale of interest, diurnal inputs may not actually be vital to capture the relevant ice dynamics, and this should be investigated in future work with detailed catchment-scale (or larger) modeling.

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 the subglacial system in the absence of a 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 present that the supraglacial hydrologic system acts as short-term storage for surface meltwater, as exhibited by the time lag of moulin inputs between models; therefore, application of an appropriate surface meltwater routing model may reduce the 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). It is conceivable, for example, that routing models could help to parameterize the storage term built into subglacial hydrology models or that some portion of this term should in fact be apportioned to supraglacial routing delays. 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 (Figs. 3–5 and S4–S10). While there is limited evidence that, within a fully coupled ice flow and subglacial hydrology model, diurnal variations 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., Fig. 7c). Subglacial channel development is key to terminating accelerated sliding in the early melt season (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 model may, in addition to providing the best representation 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 channelization behavior and temporal mean effective pressure over the 31 d simulations are quite similar across routing models (Fig. 5). Further work is needed to explore the implications of using subdaily inputs versus time-averaged inputs on long-term subglacial hydrology and ice dynamics in realistic, large domains. The case may be that the peaks and troughs of large-amplitude diurnal input variations effectively balance out to yield equivalent results as time-averaged inputs.

4.3 Influence of seasonal supraglacial drainage evolution on meltwater routing

Seasonal evolution of supraglacial drainage pattern can substantially alter surface meltwater routing and moulin discharge characteristics (Figs. 2 and S1–S3). Such changes also involve the removal and deposition of a winter snowpack (Nienow et al., 2017) and the development of a weathering crust on base ice (Cooper et al., 2018), both of which can reduce the prevalence of open-channel flow (Yang et al., 2018). 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 and river networks is a distinct advantage of the RWF (Fig. 7a and b).

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 modeling the evolution of the subglacial drainage system and the development of subglacial channels; a limited supraglacial stream and river extent will damp the diurnal range of meltwater inputs into the subglacial system and maintain a more constant subglacial pressure (e.g., Fig. 7c), which in turn results in less variation in ice motion. Furthermore, limiting diurnal variations may result in more rapid growth of subglacial channels early in the melt season (Schoof, 2010; Hewitt, 2013) and thus potentially produce a more realistic transition to primarily channelized drainage (Banwell et al., 2013, 2016).

4.4 Impact of DEM resolution on supraglacial meltwater routing

Analyzing the influence of DEM spatial resolution on the hydrological response of catchments 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 catchments to surface melt. We found that RWF routing is not affected by DEM spatial resolution, whereas SRLF routing is (Table 2, Figs. 2 and S1–S3). Notably, the RWF can only operate on high-resolution (<10 m) DEMs, whereas the SRLF can be applied to coarser-resolution DEMs as well. The 5 m SRLF UHs and resultant moulin discharges best match with RWF and SUH simulations (Figs. 2 and S1–S3). 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 behavior 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 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 (Fig. 2c and d), whereas stable flow length contributes to the RWF's better performances under different DEM resolutions (Fig. 2e and f).

4.5 Future research directions of surface-to-bed meltwater connection

A challenge remaining unsolved is mapping or simulating the seasonal evolution of supraglacial stream and river networks. The spatial extent of supraglacial stream and river networks determines the partition of interfluve and open-channel zones and thereby controls IDC hydrologic response to surface melt. The dynamic A_c thresholds provide a simple simulation of supraglacial stream and river network evolution but may not represent the realistic evolution (Yang et al., 2018). In future, a more solid relationship between RCM surface runoff and supraglacial drainage pattern should be built.

Moreover, all three routing models described in this study have limited capacity to represent complex processes of ice surface degradation (“weathering crust”) from penetration of solar radiation and associated subsurface melting (Karlstrom et al., 2014; Cooper et al., 2018). The SUH relies on catchment and stream and river network drainage pattern (Smith et al., 2017) and does not consider the possibility of shallow subsurface flow through weathering crust. While the snow component of the SRLF can handle seasonal snow cover, its bare-ice component also cannot accommodate a permeable bare-ice surface (Arnold et al., 1998; Banwell et al., 2012). The RWF, in contrast, can accommodate shallow subsurface flow through weathering crust using field-calibrated interfluve flow velocities but does not support a third option of seasonal snow cover (Yang et al., 2018). As a result, the RWF cannot be used to simulate moulin discharge for snow-covered IDCs during the early melt season. We suggest that future work should pursue an appropriate parameterization of shallow porous media flow through both seasonal snow cover and weathering crust in order to accurately describe meltwater transport in this area throughout the duration of the melt season.

Finally, the subglacial hydrology simulations in this study should be interpreted as an experimental demonstration of how different meltwater routing models influence local subglacial hydrology in the vicinity of a moulin through variations in moulin inputs. The specific effective pressures presented here are admittedly influenced by uncertain modeling parameters and boundary conditions and should be viewed as a comparison between methods, not as a representation of real behavior at the bed in these catchments. Regional simulations covering multiple IDCs in a fully coupled ice flow and 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 modeled with multiple moulins distributed in areas where borehole observations are available, or relevant quantities may be inferred from radar products (Chu et al., 2016).

5 Conclusions

Surface meltwater routing is crucial for understanding the surface-to-bed meltwater connection of the Greenland ice sheet but remains poorly studied to date. This study presents a first intercomparison of three different meltwater routing models and employs a subglacial hydrologic model to explore the impacts of their differences 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 DEMs and seasonal evolution of the supraglacial stream and 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 behavior 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 and 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 modeling.