Direct measurements of meltwater runoff on the Greenland ice sheet surface.

Meltwater runoff from the Greenland ice sheet surface influences surface mass balance (SMB), ice dynamics, and global sea level rise, but is estimated with climate models and thus difficult to validate. We present a way to measure ice surface runoff directly, from hourly in situ supraglacial river discharge measurements and simultaneous high-resolution satellite/drone remote sensing of upstream fluvial catchment area. A first 72-h trial for a 63.1-km2 moulin-terminating internally drained catchment (IDC) on Greenland's midelevation (1,207-1,381 m above sea level) ablation zone is compared with melt and runoff simulations from HIRHAM5, MAR3.6, RACMO2.3, MERRA-2, and SEB climate/SMB models. Current models cannot reproduce peak discharges or timing of runoff entering moulins but are improved using synthetic unit hydrograph (SUH) theory. Retroactive SUH applications to two older field studies reproduce their findings, signifying that remotely sensed IDC area, shape, and supraglacial river length are useful for predicting delays in peak runoff delivery to moulins. Applying SUH to HIRHAM5, MAR3.6, and RACMO2.3 gridded melt products for 799 surrounding IDCs suggests their terminal moulins receive lower peak discharges, less diurnal variability, and asynchronous runoff timing relative to climate/SMB model output alone. Conversely, large IDCs produce high moulin discharges, even at high elevations where melt rates are low. During this particular field experiment, models overestimated runoff by +21 to +58%, linked to overestimated surface ablation and possible meltwater retention in bare, porous, low-density ice. Direct measurements of ice surface runoff will improve climate/SMB models, and incorporating remotely sensed IDCs will aid coupling of SMB with ice dynamics and subglacial systems.

The production and transport of meltwater (runoff) is an important hydrological process operating on the surface of the Greenland ice sheet (GrIS). Total GrIS mass loss from runoff and solid ice dynamics (glacier calving) now exceeds ∼260 Gt/y, contributing >0.7 mm annually to global mean sea level rise (1–3). Since 2009, approximately two-thirds of this total mass loss has been driven by negative ice sheet surface mass balance (SMB) and associated runoff increases, as calculated from climate/SMB models (3, 4). This runoff passes through supraglacial stream/river networks entering moulins (englacial conduits) and crevasses that connect to the bed (5–9), temporarily influencing basal water pressures and/or ice motion (10–13) and forming a dynamic subglacial drainage system that expels water toward the ice edge and global ocean. The new dominance of runoff as a driver of GrIS total mass loss will likely persist into the future, because of further increases in surface melting (14), reduced meltwater storage in firn due to formation of near-surface ice layers (15), and possibly a waning importance of dynamical mass losses as ice sheets retreat from their marine-terminating margins (16). Therefore, the hydrological process of ice surface runoff warrants study, both for basic scientific understanding and to improve representation and/or parameterization of runoff processes in climate/SMB models.

A key uncertainty in climate/SMB projections of future GrIS runoff contributions to global sea level is that estimating runoff requires partitioning of SMB among some poorly constrained processes, with the modeled “runoff” (R) simply an error-sensitive residual of the sum of modeled meltwater production (M), rainfall and condensation, minus modeled retention, refreezing, and sublimation in snow and firn. Representation of these various elements varies by model (), but in all cases R is an error-sensitive residual that is not independently validated with in situ field measurements collected on the ice surface. Previous efforts to validate R have used proglacial river discharge (outflow) emerging from the ice edge (7, 17–20), but outflow fundamentally differs from R because it incorporates complex en- and subglacial processes that can delay, remove, or add water, including cavity storage/release, reservoir constrictions, conduit pressurization, subdaily variations in hydraulic potential gradient, basal melting, and subglacial aquifers (5, 10, 17, 21–23). Furthermore, basin delineations for proglacial river outlets have high uncertainty (7, 17, 24), are keenly sensitive to user choice of a hydraulic potential parameter [i.e., the k-value (25, 26)], and are vulnerable to water piracy between adjacent basins (27, 28). Proglacial river discharge measurements can suffer large uncertainty due to heavy sediment loads, braided channels, and mobile beds (29). In short, proglacial outflow does not confidently reflect the timing of SMB and runoff processes operating on the GrIS surface, especially at diurnal time scales.

At the present time, climate/SMB models contain little or no provision for retention and/or refreezing of runoff in bare ice (i.e., either on or below the ice surface), or for flow routing (lateral transport) of runoff over the ice surface to moulins. Instead, residual M converts instantly to R and is assumed to depart the ice surface. This is acceptable for estimating net SMB but not for estimating the timing and volume of runoff delivered to moulins, the dominant pathway linking supraglacial with subglacial hydrological systems (7, 21, 24). This, in turn, clouds understanding of the interplay between SMB and ice dynamics, especially at short time scales. Moulins inject surface runoff into a transient, subglacial hydrological system exerting primary control on diurnal to multiday changes in ice sheet basal motion and water pressure (11, 12, 30–33). Subdaily delays or lags between the timing of surface melt and basal water pressures are often used to infer capacity of the subglacial drainage system, yet supraglacial routing delays receive little or simplified treatment (9, 10, 31, 34).

Finally, solar radiation supplies most energy for melting ice on the GrIS margins and bare-ice ablation zone, followed by the turbulent flux of sensible heat (35, 36). As a result, temporal scales governing energy and mass exchange between the atmosphere and ice surface range from seconds (for turbulent eddies) to daily and monthly for net radiative surface energy balance. Because solar radiation dominates melting, it is imperative to resolve the effect of diurnal cycles in the surface energy balance on surface runoff. The diurnal time scale is especially important for runoff generation in the midelevation ablation zone, where daytime melting is interrupted by nighttime freezing (37), causing heat loss from the ice surface and potential refreezing of meltwater. Diurnal variations in runoff also influence ice dynamics, because ice motion accelerations are driven by variability in meltwater input (10, 12). Meltwater alternatively flows from subglacial channels into the distributed basal system during intervals of high supply/high pressure, and from the distributed system into channels during intervals of low supply/low pressure (33, 38, 39). This diurnal pressurization of the distributed system drives diurnal variations in ice velocity. Numerical modeling shows increases in diurnal ice motion and a slight increase in annual mean velocity when diurnal variations in surface runoff input are considered (40).

In sum, climate/SMB models are essential tools for simulating SMB runoff inputs to subglacial systems and to the global ocean (41, 42), but they currently lack validating field measurements of runoff timing and quantity, especially over diurnal time scales. To address these challenges, we present a field-based approach to measure R directly on the ice sheet surface—before en- and subglacial interferences—at the scale of a supraglacial internally drained catchment (IDC). IDCs are defined by fluvial supraglacial stream/river networks, which dominate surface drainage patterns on the southwestern GrIS (43). They have areas of order ∼101–102 km2, a geographic scale comparable to the grid cells of most regional climate/SMB models. The field procedure is demonstrated for a representative IDC having an area of 63.1 km2 (our best estimate of catchment area, with upper and lower uncertainty bounds of 69.1 and 51.4 km2, respectively), hereafter called the Rio Behar catchment in honor of the late Dr. Alberto E. Behar (Fig. 1). Spanning an elevation range of 1,207–1,381 m, Rio Behar catchment is located just below the long-term equilibrium line [∼1,500 m above sea level (a.s.l.) in this area (34)], experiences seasonal melting from June through August of each year, and is centrally located in one of the highest runoff-producing regions of the GrIS (3, 14). Our field trial was conducted in late July 2015, near the end of the peak runoff season when the region’s supraglacial stream/river networks are fully developed, yet before the onset of reduced melting in August (Fig. 2).

Fig. 1.

WorldView-1/2 satellite-derived map of Rio Behar catchment, a moderately sized (63.1 km2) internally drained catchment (IDC) centrally located in a melt-intensive area of the GrIS (Inset). From 20 to 23 July 2015, we collected 72 h of continuous in situ ADCP discharge measurements in the main-stem supraglacial river (Rio Behar) at our base camp (black star; 67.049346N, 49.025809W), ∼300 m upstream of the catchment’s terminal moulin. Measurements of ice surface ablation were collected at base camp and by the PROMICE KAN_M automated weather station. Four GPS-surveyed red tarpaulins visible in satellite and drone imagery were used as ground control points (GCP) to aid image geolocation and georectification. Eight years of topographic Rio Behar catchment boundaries, delineated from WorldView satellite stereo-photogrammetric DEMs (multishaded gray lines), establish overall catchment stability from 2008 to 2015. The 18 July 2015 DEM boundary, adjusted for small areas of stream piracy, was used for calculations presented in this study (thick black line; 63.1 km2). Manually identified stream channel heads (headwater channel incision points) mapped in the 18 July 2015 satellite image constrain minimum (green circles, inner) and maximum (red circles, outer) plausible catchment boundaries, respectively. The minimum boundary eliminates crevasse fields in the southeast catchment headwater area. Polygons bound small confirmed (red polygons) and potential (purple polygons) internally drained subareas (i.e., internal moulins) not draining to the large terminal moulin. Four small, nondraining supraglacial lakes were fully integrated into the stream/river network with no impoundment of flow. This map was created in part using DigitalGlobe, Inc., satellite imagery.

Fig. 2.

The 20–23 July 2015 field experiment (dashed lines) was timed for late July near the end of the peak runoff season, when Rio Behar catchment was bare ice, its seasonal surface drainage pattern was fully developed, and before the onset of cooler temperatures and reduced melting in August. Colored lines show daily melt rates (M) from the HIRHAM5, RACMO2.3, MAR3.6, and Point SEB climate/SMB models; melt rate is not supplied by MERRA-2.

Conceptually, our approach is simple, requiring only hourly measurements of discharge in an IDC main-stem supraglacial river (i.e., to measure the volume of runoff physically departing the source catchment) and high-resolution mapping of the IDC’s contributing upstream catchment area. Note that “runoff” has units of depth per model time step in gridded climate model output (L T−1, typically mm ⋅ d−1 or mm ⋅ h−1) but units of discharge when obtained from in situ measurements (L3 T−1, typically m3 ⋅ s−1). Remotely sensed catchment area (L2, typically km2) is required for conversion between the two units of runoff.

We measured discharge hourly in the main-stem supraglacial river of Rio Behar catchment for 72 h from 20 to 23 July 2015 by deploying an RTK GPS SonTek RiverSurveyor Acoustic Doppler Current Profiler (ADCP) from a bank-operated cableway suspended across the river immediately upstream of its descent to the catchment’s terminal moulin (). During the same period, we obtained high-resolution images from the DigitalGlobe WorldView-1 and WorldView-2 satellites (resolution 0.5 m panchromatic, 2.0 m multispectral) and from a custom-made fixed-wing drone [unmanned aerial vehicle (UAV); visible band, resolution 0.3 m]. These acquired images were used to map Rio Behar catchment boundaries, surface drainage pattern, and snow cover. Topographic divides of the catchment were delineated from a high-resolution digital elevation model (DEM) of the ice surface, derived stereo-photogrammetrically from a WorldView-1 image pair acquired 18 July 2015. The long-term stability of this catchment was established from older WorldView image pairs beginning in 2008 (Fig. 1). The 2015 topographic boundary was later manually adjusted for small areas lost (2.7 km2) or gained (0.8 km2) due to stream piracy (breaching) across divides, and for small internal subareas draining to minor internal moulins (1.6 km2). Intersection of this corrected catchment area (63.1 km2) and its maximum plausible extent (69.1 km2, identified by mapping outer channel heads; Fig. 1) and minimum plausible extent (51.4 km2, identified by mapping inner channel heads and removal of 4.1 km2 of crevasse fields; see Fig. 1 and ) with gridded outputs from the HIRHAM5, MAR3.6, RACMO2.3, MERRA-2, and Point SEB climate/SMB models enables a first direct comparison between modeled and measured on-ice R for the Rio Behar catchment (Figs. 3 and 7).

Fig. 3.

Hourly supraglacial runoff R from the Rio Behar catchment obtained from in situ ADCP discharge measurements (red) and as estimated by five climate/SMB models (color-shaded envelopes) during the 20–23 July 2015 field experiment. Observed runoff is attenuated and delayed relative to modeled runoff due to nonrepresentation of fluvial transport (routing) in current models. An exception is MAR3.6 (green), which uses a simple delay-to-ice-edge assumption, thus greatly smoothing the diurnal runoff signal. Units of R in climate/SMB models (mm ⋅ h−1) are converted to discharge (m3 ⋅ s−1) by multiplication with remotely sensed catchment area (Fig. 1), enabling direct comparison with ADCP measurements. The uncertainty bounds shown for modeled R thus reflect Rio Behar catchment area uncertainty, with centerlines denoting the optimal catchment area estimate of 63.1 km2 and upper and lower uncertainty reflecting the maximum and minimum plausible catchment area estimates of 69.1 and 51.4 km2, respectively. Error bars for in situ data are SDs calculated from multiple ADCP profiles collected within each measurement hour. Local time for Rio Behar catchment is Coordinated Universal Time (UTC) minus 2 h.

Fig. 7.

Climate/SMB model simulations compared with field measurements of (A) runoff and (B) ice surface lowering (ablation) during the 20–23 July 2015 field experiment. (A) Cumulative hourly supraglacial runoff R from the Rio Behar catchment as measured from in situ ADCP measurements (in red) and as estimated by five climate/SMB models (color-shaded envelopes). Note that values of cumulative modeled R (m3) derive from summation of hourly discharges (m3 ⋅ s−1), which are obtained by multiplying climate/SMB model outputs with the remotely sensed catchment area(s) of Fig. 1. Upper and lower uncertainty bounds in modeled R thus reflect Rio Behar catchment area uncertainty, with centerlines denoting the optimal catchment area estimate of 63.1 km2 and upper and lower uncertainty bounds reflecting maximum and minimum plausible catchment area estimates of 69.1 and 51.4 km2, respectively. Error bars (red) for in situ measurements denote the following: (A) Cumulative SDs calculated from multiple ADCP supraglacial river discharge measurements collected within each measurement hour; and (B) cumulative ice surface–lowering measurements as measured manually at 15 ablation stakes in our Rio Behar base camp (mean values also shown) and by the KAN_M AWS. Upper and lower uncertainty bounds in modeled ice ablation reflect assumptions of either solid ice (0.918 g ⋅ cm−3) or lower observed (0.688 g ⋅ cm−3) (50) bare-ice density to convert model outputs of M from units of liquid water equivalent to solid ice equivalent. The vertical dashed line in B indicates time of cessation of our ADCP discharge experiment in A. MERRA-2 is not shown in B because M is not supplied by MERRA-2. Local time for Rio Behar catchment is Coordinated Universal Time (UTC) minus 2 h.

Results

Comparison of our hourly discharge measurements () with hourly climate/SMB model outputs of catchment R quantifies the attenuation and delay of observed R delivered to the Rio Behar catchment terminal moulin (Fig. 3). Because evacuation of runoff requires physical passage through the IDC’s fluvial drainage pattern, some duration of time must pass between the timing of peak R generated across the IDC and the timing of peak R (i.e., peak discharge) received by the moulin. This duration is called “time-to-peak” (t, in hours) in traditional terrestrial hydrograph analysis (44, 45). In general, time-to-peak delays will increase for catchments having larger and/or more elongate areas and lower stream densities, with soil porosity, topographic slope, and land cover being contributing factors (45). For a given uniform depth of R generated across the catchment, larger catchments produce greater total discharge and peak discharge (Q) than do smaller catchments, due to their larger source areas. Applied to southwest Greenland, where most IDCs have areas of tens of square kilometers (43), these fluvial catchment processes are thus intrinsic to the scale of a climate/SMB model grid cell.

To demonstrate how influential fluvial supraglacial catchments are to the timing (t) and peak discharge (Q) of GrIS meltwater runoff delivery to moulins, we use our Rio Behar discharge measurements to calibrate a simple lumped (catchment-scale) morphometric routing model for use on the ablating ice surface, the synthetic unit hydrograph (SUH; ). Three advantages of the SUH routing model are that it isolates the impact of basic IDC properties (area, shape, and stream length) on t and Q delivered to the catchment outlet (here, the terminal moulin), which can all be obtained with remote sensing; it does not require use of DEMs [which are acutely sensitive to choice of a depression-filling threshold and do not always reflect true surface drainage patterns (46)]; and it is designed to be transferable to ungauged catchments.

Extension of our field-calibrated SUH to a broad-scale (13,563 km2), remotely sensed map of 799 surrounding IDCs (43) quantifies temporal and spatial heterogeneities in runoff delivery to terminal moulins due solely to differences in IDC areas, shapes, and river lengths (Fig. 4). For a theoretical unit runoff depth of 1 cm (i.e., a 1-cm-deep layer of water assumed to materialize uniformly across the ice sheet surface in 1 h), catchment-induced time-to-peak delays would range from as low as 0.4 h to as high as 9.5 h, due solely to varying IDC areas, shapes, and river lengths (Fig. 4). Peak discharges entering moulins would range from as low as 0.7 m3 ⋅ s−1 to as high as 53.0 m3 ⋅ s−1 (Fig. 4), again due solely to these basic fluvial catchment properties that are not currently represented in climate/SMB models.

Fig. 4.

Application of our field-calibrated Synthetic Unit Hydrograph (SUH) routing model to 799 remotely sensed IDCs on the southwest GrIS [gray borders; mapped previously from a 19 August 2013 panchromatic Landsat-8 image (43)] illustrate how fluvial, supraglacial IDCs impart spatially heterogeneous modifications to meltwater runoff delivered to terminal moulins and hence the bed. Each IDC contains a remotely sensed supraglacial river (not shown for visual clarity) terminating in a major, catchment-terminating moulin. These theoretical SUH maps assume a spatially uniform, 1-cm-deep layer of meltwater released over a duration of 1 h and isolate the influence of remotely sensed IDC area, shape, and river length on (A) time-to-peak delays of peak runoff arrival at each catchment’s terminal moulin (t, in hours) and (B) magnitude of peak discharge received at each catchment’s terminal moulin (Q, m3 ⋅ s−1). More realistic maps, forced by climate/SMB models, are shown in Fig. 5 and .

Fig. 5.

Supraglacial IDCs modify the timing and magnitude of runoff delivered to terminal moulins, as demonstrated here at 1400 local western Greenland time on 21 July 2015 using (A) MAR3.6, RACMO2.3, and HIRHAM5 climate/SMB model outputs of corrected meltwater production (M′; see ) to estimate (B) instantaneous area-integrated runoff and (C) more realistic, SUH-routed runoff. MERRA-2 is not shown because it does supply M. Point SEB is not shown because its output is not gridded. The boundaries of 799 IDCs (gray borders) were mapped previously from a 19 August 2013 panchromatic Landsat-8 image (43). Each IDC contains a remotely sensed, moulin-terminating supraglacial river (not shown for visual clarity). Climate/SMB model output M′ has units of water depth equivalent (mm ⋅ h−1), which converts to runoff in discharge units (m3 ⋅ s−1) following multiplication with intersected IDC catchment boundaries (B and C). The black star at ∼67N, 49W denotes the Rio Behar IDC. In both B and C large IDCs enable large moulin discharges above 1,500 m a.s.l. elevation, despite lower overall melt rates. SUH routing (C) yields lower peak moulin discharges at this time of day than instantaneous area-integrated runoff (B), because SUH requires more time for runoff to travel through fluvial supraglacial stream/river networks. A companion nighttime version of this figure 10 h later (00:00 on 22 July; see ) shows the opposite effect, with shutdowns in A and B but high moulin discharges in C.

A more realistic scenario, using climate/SMB model outputs of melt production M and a Gamma function to synthesize each IDC’s unique SUH (47) (), yields similarly heterogeneous spatial patterns not present in gridded climate model output (Fig. 5 and ). These heterogeneities include large discharges (>20 m3 ⋅ s−1) entering moulins at high elevations on the ice sheet (>1,500 m a.s.l.) despite low melt rates there, due to the presence of large IDCs (43, 48). Importantly, peak moulin discharges are significantly reduced if climate/SMB model output is subjected to unit hydrograph theory (Fig. 5), rather than the practice of instantaneously aggregating model output within each IDC (33, 49) (Fig. 5). The opposite is true at night, when modeled melt and instantaneous area-aggregated runoff shut down but SUH-routed runoff is high (). Averaging across all 799 IDCs (including many small catchments) reduces Q by 13.5 ± 10.0% if climate/SMB model output is subjected to SUH routing (Fig. 6). Diurnal variability in Q is reduced by 15.1 ± 12.5%, and the mean timing delay between peak melt production and peak moulin discharge lengthens by 2.9 ± 2.8 h. For the larger IDCs (>30 km2, n = 122), for which routing delays are greatest, these averages increase to 30.4 ± 9.1%, 37.0 ± 12.0%, and 5.1 ± 4.6 h, respectively.

Fig. 6.

Comparison of SUH-routed runoff (Fig. 5) with instantaneous area-integrated runoff (Fig. 5) for all 799 IDCs: (A) peak moulin discharge; (B) diurnal difference between maximum and minimum moulin discharge; and (C) time delay between peak melt production across the catchment and peak discharge received at the terminal moulin. Applying SUH routing to climate/SMB model output yields lower peak discharges, suppressed diurnal variability, and delayed, asynchronous timing of peak runoff delivered to catchment-terminating moulins.

Although these numbers should be viewed cautiously because our SUH model depends, in part, on parameters calibrated only at the Rio Behar catchment, a successful retroactive application of SUH to the IDCs of two older field studies (8, 32) is encouraging (). Depending on choice of input climate/SMB model, SUH-estimated peak runoff times for a 1.1-km2 IDC nearly 300 km distant from the Rio Behar catchment range from 16:00 to 20:00 (local Greenland time), comparable to 16:30–17:00 observed in field observations acquired in August 2009 (8, 32) (). For an 18.2-km2 IDC ∼14 km distant from Rio Behar, SUH-estimated peak runoff times range from 17:00 to 22:00, comparable to field measurements of 18:00–20:00 acquired in late June/early July 2011 (32). Such independent reproducibility of runoff timing delays measured at other times and sites on the ice sheet suggest utility of SUH elsewhere on the southwest GrIS ablation zone. However, collection of additional supraglacial discharge datasets, especially from large IDCs and colder regions, is needed for further calibration and validation of the SUH approach.

With regard to the absolute magnitudes of measured versus modeled runoff, comparison of our cumulative ADCP discharge measurements with cumulative modeled R over our 72-h field experiment found that climate/SMB models overestimated R by +21 to +58% for this particular location and time on the ice sheet (for a five-model average, assuming lower and upper constraints on watershed extent, respectively). Taken separately, four of five models overestimated R (Fig. 7). Similarly, four of four models (for which melt M is available) overestimated ice surface lowering (ablation), if their outputs of M are compared with in situ ice surface–lowering measurements collected from 15 ablation stakes at our base camp and sonic surface–lowering data from the nearby PROMICE KAN_M automated weather station (AWS) (Fig. 7 and ). This conclusion holds regardless of whether the density of solid ice (0.918 g ⋅ cm−3) is used to convert M to units of ice thickness equivalent, or a lower, near-surface ice density (0.688 g ⋅ cm−3) averaged from 10 shallow cores drilled at our base camp (50). Point-based ablation measurements have known limitations (51), but both field datasets display less ice surface lowering than modeled M (Fig. 7), similarly to how the models overestimate R (Fig. 7).

One interpretation of Fig. 7 is that the models overestimated M, and hence R. However, our comparison of modeled versus in situ AWS surface energy balance () reveals that modeled energy balance components closely matched in situ AWS measurements. In general, RACMO2.3 albedo, radiation, and turbulent fluxes track AWS observations too well to advance model overestimation of surface energy receipt as the leading explanation for model overestimations of ice ablation and R (). For example, the radiative effects of clouds (52) may have contributed slightly to model overestimation of R during the third day of the field experiment, but not the first 2 d when the sky was clear (). Importantly, the Point SEB model is driven purely by AWS measurements, yet similarly overestimates observed surface ablation and R like the other, reanalysis-driven models (Fig. 7).

All of this suggests some meltwater loss or retention process that is external to the “skin” surface energy balance allocated to the top of the ice surface by most models. We hypothesize that subsurface melting (53) and subsequent retention and/or refreezing of meltwater in porous, low-density bare ice [called “weathering crust” (6, 50, 54)] may contribute to or explain the observed discrepancies between modeled M and R and measured ice surface lowering and supraglacial river discharge, respectively. Runoff infiltration into crevasses (8) cannot explain the observed runoff deficit, because crevassed areas are eliminated from our minimum bounding catchment map (51.4 km2) and are thus already included in the lower model uncertainty bounds of Fig. 7 (and Fig. 3). Although the possibility of additional, missed leakage cannot be fully ruled out, there is no evidence for this in our high-resolution UAV imagery (). Missed meltwater retention in seasonal snow also seems unlikely: the climate/SMB models indicate bare ice, and snow classifications from our UAV mapping and two WorldView-2 images confirm that Rio Behar catchment had <6.5% snow cover at the time of our field experiment, and perhaps as little as 0.9% (). Remotely sensed retrievals of lake volume storage rule out the possibility of runoff impoundment in four supraglacial lakes contained within the Rio Behar catchment (). The remaining hypothesis (i.e., of water retention/refreezing in the bare-ice weathering crust) is explored further in Discussion and Conclusions and in the .

Regardless of mechanism, a first-order, empirical correction for any missed retention processes and/or model overestimations of M for the Rio Behar catchment during our field experiment is supplied by a set of empirical, model-specific runoff coefficients relating observed runoff R to modeled melt production M for HIRHAM5, MAR3.6, RACMO2.3, and Point SEB (). No values are supplied for MERRA-2 because M is not an output of this model. Although these coefficients are computationally identical to how runoff coefficients are calculated for terrestrial catchments (i.e., river discharge divided by catchment water input), they also include any missed model over- or underestimation of M and are more properly treated as correction factors for climate/SMB models instead of traditional runoff coefficients. For bare-ice surface conditions similar to those observed at the Rio Behar catchment during our field experiment, these correction factors may be multiplied by M to obtain alternate, lower estimates of M (here termed effective melt M′) in addition to standard model output.

Discussion and Conclusions

Although the field protocol presented here is currently logistically impractical for sustained monitoring or deployment at numerous sites, it offers a useful, and perhaps only, direct way to independently measure supraglacial R for validating climate/SMB models used to simulate ice sheet runoff and associated inputs to subglacial and marine systems. Our provision of field-calibrated, model-specific runoff coefficients and SUH parameters offers an initial step in this direction, enabling generation of SUHs, peak moulin discharges (Q), and runoff time-to-peak delays to moulins (t) from standard climate/SMB model outputs of melt production M () at a time of maximum drainage efficiency on the ice sheet surface.

Our successful retroactive testing of SUH runoff timing delays against in situ observations of two earlier field studies (8, 32) conducted in different years, locations, and elevations than Rio Behar catchment suggests plausible transferability of SUH to other areas of the GrIS ablation zone. One reason for this success may be that only three SUH parameters (C, C, and m) require in situ calibration; the others (t and h) derive purely from remotely sensed catchment characteristics and are thus adjusted individually for each IDC. That said, further field experiments are needed at other locations and times on the ice sheet to derive additional runoff coefficients and SUH parameters for differing surface conditions. Hydrological measurements from Haut Glacier d’Arolla, Switzerland, for example, suggest that earlier in the runoff season the presence of snow also suppresses diurnal contrasts and introduces delays between peak melt production and peak moulin discharge (55). Similarly, the seasonal evolution of supraglacial stream/river drainage networks may influence early-season runoff coefficients and the values of C, C, and m presented here due to lower stream density and/or temporary retention of runoff in slush and seasonal snow (43). Note that the most likely outcome of these processes would be to further delay runoff delivery to moulins (Fig. 6), further suppress diurnal variability (Fig. 6), and further suppress peak moulin discharges (Fig. 6), rendering conservative our scientific conclusions about the influence of fluvial supraglacial catchments on meltwater delivery to moulins and the bed.

The field measurements and SUH calculations presented here illustrate the critical importance of IDCs in modulating the timing and magnitude of runoff evacuated off the ice surface to moulins (Figs. 4 and 5 and ). Previous studies have shown the importance of filling and draining supraglacial lake basins (20, 30, 56), but even in the absence of lake basins, runoff becomes unevenly redistributed over space and time due to water collection and transport through fluvial supraglacial stream/river catchments. Based on our observed values of C and C, these catchment-scale processes on ice are not unlike those on land (), despite known hydraulic differences between supraglacial and terrestrial channels (57).

Because IDC areas vary greatly and moulins convey meltwater quickly to the bed (32), the timing and volume of surface runoff received at the bed are thus arrhythmic in time and heterogeneous in space, unlike outputs from gridded climate/SMB models (Fig. 5 and ). Subdaily time lags between surface climatology and melt-induced ice motion are established first on the ice surface, which might otherwise be attributed to en- and/or subglacial delays or modes in melt-induced ice motion (8, 11, 12, 25, 30, 56), basal pressure (23), or subglacial drainage capacity (34). Diurnal variability in moulin discharge is lower than that of climate/SMB-modeled runoff fields, potentially reducing rates of inferred subglacial channelization (40). Where the diurnal variability of meltwater delivered to the bed is dampened by surface routing delays, there should be an impact on ice-sliding velocities, especially at higher elevations on the ice sheet. Using climate/SMB runoff to drive ice dynamics models in such areas could thus overestimate diurnal subglacial pressure variability, leading to small overestimations in the diurnal range of ice velocities and perhaps annual mean velocity as well. Conversely, large IDCs have the capacity to amplify moulin discharge, including at high elevations where melt rates are low but IDCs are large (43, 48), especially if moulins are first initiated through hydrofracturing and drainage of interior-advancing supraglacial lakes (21, 24, 58) then subjected to extreme and/or sustained melt events (59). In sum, the supraglacial drainage pattern on the GrIS surface influences a host of important subglacial processes, especially at short time scales.

Our finding that modeled and observed surface energy balances largely agree (), yet both overestimate observed ice surface lowering and runoff (Fig. 7), leads us to hypothesize that subsurface melting and delay/retention/refreezing of meltwater in porous, low-density weathering crust may be an important bare-ice physical process not represented in the climate/SMB model simulations presented here. Shortwave radiation penetration and subsurface melting of bare ice certainly promotes the development of weathering crust (6, 53, 54) at our study site (), which is characterized by abundant cryoconite holes and porous, water-saturated, low-density bare ice at least 1.1 m deep (50). Ablating weathering crust typically experiences less surface lowering than expected from skin surface energy balance calculations alone, due to internal melt within the subsurface ice matrix (50, 60, 61). Any meltwater retained within this porous medium—for example, due to deepening of the crust, enlargement of cryoconite holes, or enlargement of pore space volume—would result in model overestimation of R because current modeling schemes do not permit water retention in bare ice. Moreover, any refreezing of this meltwater (which we observed nightly during the field experiment) requires that it remelt to become true runoff, consuming additional melt energy not currently allocated in energy balance models for the bare-ice zone. Any model that correctly quantifies surface melt energy but does not simulate these processes will overestimate both ice surface lowering and runoff ().

Although mismatched scale and timing preclude direct comparison of our field results with GRACE (Gravity Recovery and Climate Experiment) satellite gravity data, we note in the that two previously published, sector-aggregated GRACE observations similarly show less actual mass loss than simulated by climate/SMB models (SMB-D) in some key melt-intensive sectors, including ours in southwest Greenland (62, 63) (). However, we are reluctant to draw general conclusions about climate/SMB model performance at other times or locations on the GrIS due to the short duration and small geographic area (relative to model domains) of our field experiment. The observed spread in modeled runoff estimates for the Rio Behar catchment (Fig. 7) is consistent with a broader intercomparison of modeled outputs across the GrIS, including heightened model uncertainty in the ablation zone (64). New field experiments are needed to determine how to refine climate/SMB model simulations of ice surface ablation and runoff in the bare-ice zone, as well as remote-sensing SMB estimates that use satellite/airborne altimetry measurements of ice surface lowering.

Regardless of absolute magnitudes of R, the timing and amplitude of meltwater runoff are clearly modified by fluvial geomorphology and fluvial catchment processes operating on the GrIS surface. Lateral flow routing through internally drained catchments predictably delays the arrival, reduces the peak discharge, and suppresses the diurnal variability of R entering moulins. Large catchments yield high moulin discharges, even at high elevations where overall melt rates are low. These realities, together with possible delays/retention/refreezing of runoff in bare-ice ablation zone weathering crust, signify that supraglacial drainage processes critically preconfigure the timing and flux of meltwater delivered to the bed. Incorporating fluvial catchments, hydrological theory, and field calibrations into ice sheet models should improve coupling of SMB with subglacial and marine systems.