Rapid basal melting of the Greenland Ice Sheet from surface meltwater drainage.

Subglacial hydrologic systems regulate ice sheet flow, causing acceleration or deceleration, depending on hydraulic efficiency and the rate at which surface meltwater is delivered to the bed. Because these systems are rarely observed, ice sheet basal drainage represents a poorly integrated and uncertain component of models used to predict sea level changes. Here, we report radar-derived basal melt rates and unexpectedly warm subglacial conditions beneath a large Greenlandic outlet glacier. The basal melt rates averaged 14 mm ⋅d-1 over 4 months, peaking at 57 mm ⋅d-1 when basal water temperature reached +0.88 ∘C in a nearby borehole. We attribute both observations to the conversion of potential energy of surface water to heat in the basal drainage system, which peaked during a period of rainfall and intense surface melting. Our findings reveal limitations in the theory of channel formation, and we show that viscous dissipation far surpasses other basal heat sources, even in a distributed, high-pressure system.

The flow of ice sheets and glaciers is controlled by basal motion, which takes place through some combination of hard bed sliding (1–3), sliding-induced cavity formation (4, 5), and deformation of subglacial sediment (6, 7). All forms of basal motion require a thawed thermal state in order to be substantial (8), with more heat produced at (or delivered to) the bed than lost through conduction into the colder ice above (9). Basal motion is also strongly influenced by the way in which hydrologic systems evacuate meltwater (10), which is produced basally as well as at the surface. In settings where surface meltwater is transferred to the bed, drainage is often expected to occur through large channels, which become increasingly efficient in terms of discharge when they grow in size (11, 12). The resulting decrease in water pressure produces arborescent networks in which larger channels capture water from their less efficient surroundings, including smaller channels as well as water stored in small cavities (13), thin films (14), or porous sheets (15). In Greenland, channelized basal drainage has been observed as far as 30 km inland from the land-terminating southwest margin (16), and recent studies show that channels may also form beneath marine-terminating glaciers (17, 18), which drain 88% of the ice sheet (19). However, the evolution of basal drainage system efficiency, and channels’ ability to form under thick ice, remain highly uncertain (20, 21).

The central process in channel formation is energy dissipation through turbulence and viscous resistance in the water flow, which should make small cavities or sheets unstable (22) and result in channel growth until wall melting balances creep closure (11). In the classic theory of steady-state water flow in subglacial channels, Röthlisberger (11) assumed that heat transfer occurs instantaneously and that the temperature of water is fixed at the pressure-dependent melting temperature of the ice. Nye (23), followed by Spring and Hutter (24), extended this theory to consider transient water flow with temperature-dependent heat transfer in Icelandic subglacial outburst floods (Jökulhlaup) (25). However, with a paucity of data to confirm how energy is dissipated in basal drainage systems more broadly, Röthlisberger’s simpler theory has become a cornerstone in hydrologic glacier models today (26).

Here, we report a time series of radar-derived basal melt rates (BMRs) together with contemporaneous, colocated borehole records showing basal water pressure and temperature of water beneath a large Greenlandic outlet glacier. The reported BMR is unprecedented because it is two orders of magnitude higher than previous estimates for an ice sheet (27, 28) and comparable to the rate of meltwater generation at the surface. The high magnitude is corroborated by independent borehole records, which capture the temperature dependency of heat transfer and viscous energy dissipation in the basal drainage system of an ice sheet.

Results

To quantify basal melting, we used an autonomous phase-sensitive radio echo sounding (ApRES) instrument, which has millimeter range precision (29), to track the vertical displacement of internal layers and the ice–bed interface of Sermeq Kujalleq (Store Glacier), West Greenland ( and ). The ApRES instrument was installed 30 km inland from the glacier front, at site S30 (Fig. 1), where the local ice thickness was estimated to be 604 m to 606 m, and ice properties and basal conditions are well constrained (30–32). The ApRES instrument was configured to obtain 4-hourly measurements in unattended mode and operated continuously from 3 August to 4 December in 2014.

Fig. 1.

Vertical ice deformation and BMRs observed by radar beneath Store Glacier, Greenland. (A) Location and LandSat-8 image (acquired 1 July 2014) showing Store Glacier and the S30 study site in Greenland. Black line shows the central flowline. (B) Daily vertical deformation rate (VDR) of the ice column at S30 derived by tracking internal layers’ displacements over time (positive when ice column thickens). Light orange shading highlights periods with surface melt. Dark orange shading shows a cyclonic rainfall event with intensified surface melting due to warm atmospheric conditions. Red bars represent the SE. (C) Daily BMR obtained by subtracting total vertical ice deformation shown in B from phase-sensitive measurements of the ice column thickness. Red bars indicate SE calculated as the square-root sum of error terms. Days with insufficient samples (green star) were excluded from the time series.

Radar-Derived BMRs.

Daily BMR was calculated following the same approach used in studies of Antarctic ice shelves (33, 34). In a two-step approach (32), we started by fitting a vertical velocity model to the observed displacement of internal reflectors (). The best fit throughout the entire 4 month period was a linear regression model (32) (), resulting in positive vertical deformation rates (thickening) averaging mm ⋅d−1 during the observational period. In the second step, we subtracted the strain-induced thickening from the observed displacement of the ice–bed interface, which was identified clearly at a depth of 604 m to 606 m (). The resulting BMR was positive and persistently high, especially during summer, when the average rate was mm ⋅d−1 (Fig. 1). We also recorded a distinct peak in basal melting ( mm ⋅d−1) on 18 August, coincident with high surface melt rates of 56 mm w.e. d−1 during a rainfall event that brought 80 mm of precipitation over 6 d (Fig. 2). In winter, BMR was notably lower (9.8 ±0.9 mm ⋅d−1) and less variable (Fig. 1).

Fig. 2.

Borehole records from site S30 on Store Glacier. (A) Basal water pressure (p) recorded in a borehole drilled to the bed at site S30. The ice overburden pressure (p) is derived from precise colocated ApRES measurements of ice column thickness. The difference between p and p is the effective pressure (N) used to estimate the basal shear stress and frictional heat produced at the bed. The p fraction of overburden pressure () is calculated from the mean radar-derived ice thickness. Vertical gray lines denote approximate time of sensor installation. Inset shows dampened diurnal variations in p (blue line) together with strong diurnal fluctuations in surface air temperature (green line), for period marked by black box. (B) Borehole-installed temperature records from hydrological system at thawed glacier bed (T1) and sensors which froze into basal ice immediately above the bed (M1) and ∼3 m (T2) and 7 m (T3) higher. The horizontal dashed line indicates pressure-dependent water–ice phase transition temperature ( –0.40 ∘C). (C) Ice surface velocity (U, right axis) recorded from GPS installed at drill site together with rates of basal motion (U, right axis) obtained by subtracting ice deformation recorded as tilt in the borehole. is the mean ice velocity after the melt season has ended. Stacked bar plot (left axis) shows surface melt recorded by an automatic weather station at the drill site (dark blue) and additional precipitation (light blue) derived from NCAP/NCAR reanalysis data. Borehole records shown in A and B and glacier velocity shown in C were adapted from Doyle et al. (30) published under CC-BY license (https://doi.org/10.1002/2017JF004529).

Quantifying Sources and Sinks of Heat.

The BMR of a grounded ice sheet has not previously been observed or calculated at the precision and daily resolution presented here. We find that the BMR beneath our field location on Store Glacier is two orders of magnitude higher than previous estimates of 0.10 m ⋅y−1 derived from airborne radio echo sounding profiles and attributed to high geothermal heat flux in the central Greenland interior (27). To understand why our BMR is so much higher than previous estimates, we quantified the most widely recognized sources and sinks of basal heat in our study area. As an initial analysis (see breakdown in and Eq. ), we included heat sourced from the geothermal heat flux (0.06 W ⋅m−2) and frictional heat (0.9 W ⋅m−2 to 2.6 W ⋅m−2, depending on sliding speed and basal shear stress) and an upward conductive heat loss into the ice base (-0.060 W ⋅m−2) (Fig. 3). Contemporaneous in situ measurements of basal conditions in boreholes drilled to the ice base next to the ApRES (32) enabled all of these contributions to be constrained by direct observations (), with the exception of geothermal heat flux, which was inferred from crustal thickness (35). This initial heat budget analysis accounts for basal melting of 0.12 mm ⋅d−1 to 0.30 mm ⋅d−1 (Fig. 3), which is two orders of magnitude lower than the radar-derived BMR (Fig. 1). Using an enthalpy-based formulation to include additional heat stemming from liquid water in basal ice supplied only another 0.29 W ⋅m−2 (Fig. 3, Eqs. and , ).

Fig. 3.

Sources and sinks of energy and BMR estimates. (A) Theoretical BMRs at S30 derived from estimates of heat lost by conduction into basal ice and sourced from geothermal heat flux, friction along ice base, enthalpy, and viscous heat dissipation when surface water is routed along the bed. (B) Magnification of A to illustrate small magnitude of contributions other than viscous heat dissipation. (C) Corresponding measurements of BMRs from ApRES.

Viscous Energy Dissipation in the Basal Drainage System.

Our study site is in an area of active subglacial drainage, with delivery of surface meltwater to the glacier bed through supraglacial lake drainage and hydrofracture resulting in moulins (36). We therefore explored the possibility that the high local BMR is driven by the heat generated by loss of gravitational potential energy as surface meltwater descends to, and flows at, the glacier bed (37). We derived the major drainage paths from hydrologic potential gradients established from surface and bed elevation datasets (38) and calculated the energy balance of surface meltwater traveling beneath the glacier (Eq. , ). In this model, all routed water flows down the hydrologic potential gradient (Eq. , ). As forcing, we used runoff from the RACMO2 regional climate model to prescribe daily inputs of surface water during the 2014 summer melt season (39). With a highly crevassed surface limiting the extent to which meltwater is transported supraglacially (Fig. 1), we made the simplifying assumptions that surface water reaches the bed in the grid cell in which it is produced and that all energy is subsequently dissipated as heat along basal drainage paths. At a spatial resolution of 500 m, we find a close overall agreement between modeled BMR in the central drainage path that passes our study site (S30 on Fig. 3) and the observed BMR (Fig. 3). Taking the rainfall event on 18 August, for example (Fig. 4), the model predicts 54 mm of basal melt (Fig. 3) in the basal drainage path near site S30 compared to the observed 57 mm on that day (Fig. 1). While model resolution does not change the routing of water according to hydrologic gradients, increasing (decreasing) the resolution will increase (decrease) the modeled BMR because water is routed to smaller (larger) grid cells. The close agreement between our observations and the model at 500-m resolution may reflect the approximate area over which the ice is in contact with flowing water. However, we note that the model is simple and does not feature all of the hydrological processes involved with subglacial drainage. We also cannot rule out that energy exchanges during the water’s descent to the bed might reduce the energy available for basal melting. If the latter is the case, our model would need a somewhat finer resolution in order to reproduce basal melting at the observed rate. We can nevertheless conclude that the area over which ice is in contact with flowing water at site S30 probably is on the order of some hundreds of meters.

Fig. 4.

Basal drainage and viscous heat dissipation. (A) Modeled accumulation of water in the basal drainage system of Store Glacier when hydrologic model transfers RACMO2 surface runoff on 18 August along the bed of the glacier catchment (blue colors). Spatial resolution of model is 500 m. Gray colors show the ice sheet’s surface elevation. (B) Basal melt from modeled viscous heat dissipation in the basal drainage system on 18 August. Red colors denote basal melting, and blue colors denote basal freezing, which occurs when energy from viscous dissipation alone cannot raise the temperature of water to the pressure-dependent phase transition. Circular Insets show water accumulation (A) and corresponding high basal melt rate from viscous heat dissipation (B) in a major subglacial drainage path (white dot) passing near site S30 (+), where ApRES/borehole records were obtained.

Discussion

The overall agreement between measured and modeled BMR after accounting for viscous heat dissipation shows that surface meltwater is a vast, yet overlooked, energy source. While enhanced basal melting from surface meltwater delivered at the bed has previously been inferred from ice-penetrating radar profiles in which internal reflectors dip toward the bed near the injection point (40), our ground-based measurements with ApRES and borehole thermistors provide direct measurements of the process’s magnitude and consequent melt rate. A recent large-scale study by Karlsson et al. (41) estimates GrIS basal ablation to be 5.2 ±1.6 Gt of ice per year, but notes that the spatial variability in this melting is high and still unconstrained, especially along subglacial drainage pathways where energy from the surface is likely to be focused. Our study resolves this uncertainty by showing that ice melting along basal drainage pathways beneath Store Glacier can reach levels similar to those recorded at the glacier’s surface in response to solar heating (Fig. 1 vs. Fig. 2). This discovery raises important questions concerning viscous dissipation and drainage efficiency. According to the classic theory of water flow in glaciers, flow of water at relatively high pressure in small cavities should be unstable and revert to relatively low-pressure flow in channels (42), with channel size reaching a steady state when wall melting by viscous heat dissipation exactly balances creep closure of the conduit (11). Central to this theory are two commonly used simplifying assumptions, which are 1) that the temperatures of the water and the ice wall are the same, fixed at the pressure-dependent phase transition temperature, and 2) that heat generated by viscous dissipation is used instantaneously either to melt the conduit walls or to keep the water temperature at the melting point (11). Although the assumed instantaneous heat transfer is practical and widely used (26), there is a physical inconsistency between assuming that the temperature of the water and the ice wall are equal and requiring that viscous heat dissipation in the flowing water leads to the instantaneous melting of those walls. Below, we develop further the implications of this contradiction.

Basal Heat Transfer.

In general, the rate at which heat is exchanged between a solid surface and a liquid flowing in contact with it is proportional to the temperature difference between the two (43). The first-order approximation is , where Q is the heat flux, h is the heat transfer coefficient, T is the bulk water temperature, and T is the pressure-dependent water–ice phase transition temperature of the conduit walls (44). Subglacial water flow can therefore either cause melting of conduit walls or have its bulk temperature equal to the temperature of the walls, but not both at the same time.

The heat transfer responsible for the high BMR we record at Store Glacier can be explained from colocated borehole temperature records at site S30 (30). The lowermost temperature sensor in this borehole (T1 in Fig. 2) may provide clear evidence of viscous heat dissipation in the basal drainage system, where temperatures ranged mostly between 0.2 ∘C and 0.8 ∘C over several months. The peak basal temperature of 0.88 ∘C occurred shortly after the late August rainfall event, which also resulted in a sharp rise in electrical conductivity (30), possibly as a result of ionic enrichment associated with an increased suspended sediment load in the turbulent flow. While the overlying temperature sensors (T2, T3, ) froze in and cooled, T1 remained warm and showed no sign of freezing (Fig. 2). The “warm” T1 record cannot be explained from measurements made below the ice base, because the geothermal heat flux is not sufficiently high, nor from mechanical friction alone (). Such warm conditions so close to the ice base contradict the simplifying assumption of instantaneous heat transfer, which dictates that temperatures at the base of ice sheets should be effectively bound by the pressure-dependent phase transition temperature, here –0.40 ∘C. Yet, previous studies have shown that water flowing through a glacier can sustain temperatures well above the freezing point (44–47). Indeed, the equilibrium water temperature reached when viscous heat dissipation in the water equals the heat flux into the surrounding colder ice can match our measured value of 0.88 ∘C in a conduit where the pressure-dependent melting point is –0.40 ∘C, for example, if the hydraulic radius is 2 m and the gradients in local elevation and hydraulic head both reach 5∘ (Eq. , ).

To estimate the heat transfer between the subglacial water and the basal ice at our study site, we assume that the heat flux from the highest BMR (57 mm ⋅d−1) is provided by the warmest water ( 0.88 ∘C), and vice versa (i.e., the lowest BMR of 10 mm ⋅d−1 for 0.19 ∘C). This gives a heat transfer coefficient of ∼60 W ⋅m−2 ⋅∘C−1 to 170 W ⋅m−2 ⋅∘C−1, which can be achieved for a large range of water depths of 0.01 m to 10 m, while the water velocities should be on the order of 1 cm ⋅s−1 to 10 cm ⋅s−1 (Eqs. and , ). Our observations also indicate that heat generated mechanically is advected downstream. To capture this effect, instantaneous heat transfer cannot be assumed, and hydrological glacier models will instead need to solve the energy transport using heat transfer coefficients that control the rate at which heat generated by mechanical energy dissipation is transferred to the walls (47). While the latter was included in original work on Jökulhlaup by Nye (23), Spring and Hutter (24), and Clarke (25), the energy transport equation has, so far, not been implemented widely apart from the special case of Icelandic outburst floods (23, 48).

Basal Drainage.

Our measured BMR indicates that the basal drainage system would require a dimension at least as large as the 25-m spatial footprint of the ApRES (). The high BMR could occur in an efficient system in which channels are much wider than high (49), or, alternatively, in a system of canals (12), which are nonarborescent but theoretically stable under ice sheets underlain by sedimentary beds (50). While channels can also form over sedimentary glacier beds (12), a system of canals is more consistent with the 45-m-thick layer of unconsolidated sediments reported at site S30 (31), as well as high basal water pressures observed close to (>90%) the ice overburden throughout the period of observation (Fig. 2). However, the high BMR may also occur if a thin film or water sheet grows larger than the laminar–turbulent transition and is stabilized by clasts protruding into the ice (51). The offset between minor diurnal peaks in basal water pressure and the daily maximum surface air temperature (Fig. 2) is consistent with the latter or canals forming at site S30, while channels may develop closer to the terminus (52).

Surface Driver of Basal Melting.

The 2014 melt season (1 June to 31 August) produced an average of 16 × 106 m3 of surface meltwater per day in the catchment (). Assuming that all water drained to the bed, the power for basal melting by viscous heat dissipation would range from 0.66 GW on 11 August when BMR was recorded at 4 mm ⋅d−1, to 8.6 GW on 18 August when BMR was 57 mm ⋅d−1 (Fig. 1). The recorded peak BMR corresponds with a peak in daily runoff of 80 × 106 m3 (). In comparison, turbines in the Three Gorges Dam in China produce 0.7 GW of power with a peak flow rate of . The average power generated when surface water is transferred to the bed at Store Glacier (∼3 GW) dwarfs all other basal heat sources, and is comparable with the power generated in the world’s largest man-made dams. Our observations, therefore, yield direct evidence for the sustained impact of the process proposed by Mankoff and Tulaczyk (37), who argued that present-day runoff on the Greenland Ice Sheet may deliver 66 GW to its base and that the power available for viscous heat dissipation could increase to 110 GW or 320 GW by 2100 under Intergovernmental Panel on Climate Change climate scenarios Representative Concentration Pathway 4.5 (RCP4.5) and RCP8.5, respectively.

Although we have assumed that runoff is directly transferred to the bed through the numerous crevasses and moulins on Store Glacier (36), analysis of echo strength and attenuation recorded in our ApRES measurements from site S30 indicates that some of the runoff is instead stored englacially (53). This effect is indirectly included in our study, because the modeled runoff used to quantify viscous dissipation excludes meltwater retention at the surface (39). This retention can be seen in the difference between our measurements of surface ablation and modeled runoff at site S30 (). We also note that the seasonally stored melt volume is a small fraction of the total runoff and that impounded water inferred previously at site S30 has a short (1 y) residence time (53). However, deep crevasses containing refrozen meltwater show retention of meltwater to greater depths than previously reported (54). Hence, we cannot rule out the possibility that some of the water we infer to reach the bed may be stored englacially, resulting in a continued delivery of water to the bed after the melt season has ended or cryohydrologic warming in the upper part of the glacier if the water refreezes. The former may explain why the observed BMR remains higher than modeled values when the melt season has ended (Fig. 1), and why oceanographic measurements show a sustained delivery of fresh basal meltwater from Store Glacier into Ikerasak Fjord even in winter (55). Cryohydrologic warming will, however, be important in places where crevasses are deep but cannot penetrate the full ice thickness (54). The extent to which crevasses will fully or partially penetrate the glacier is controlled by the mean state of stress in the ice which influences the ability of water to pond (56).

The high BMRs from ApRES are consistent with theoretical estimates of viscous dissipation in the basal drainage system. They are also supported by independent measurements of warm subglacial water. Our observations therefore call into question the assumption of the thermodynamic equilibrium of water flowing through and beneath glaciers. We have shown that the temperature of meltwater flowing in glacial conduits need not be at the pressure melting point and that heat is not transferred between meltwater and ice instantaneously as assumed by most theoretical studies (11). This means that hydrological glacier models based on Röthlisberger’s theory in general may be overestimating rates of conduit enlargement through melting, and, conversely, that viscous dissipation in distributed systems may not result in channel formation even when discharge grows large. This disruption may occur when ice slides rapidly over a rough bed, creating a setting in which protruding sediments or clasts stabilize a turbulent water sheet (51), or, alternatively, if fluvial erosion under the glacier produces a nonarborescent system of canals (50).

Although warming of subglacial water shows that only a portion of the available energy goes into melting, the transfer of water from surface to bed makes channels likely to form near supraglacial lakes, which drain rapidly (36) and form moulins (57) that often continue to deliver large fluxes of water from considerable heights to the glacier bed. While viscous dissipation in that water will promote the growth of channels at these locations in general, a delay between meltwater flux and channel growth may lead to a transfer of water and energy into the enveloping distributed system, creating the conditions we report herein. Numerical modeling of the hydrological system beneath Store Glacier shows that such nonequilibrium conditions are likely to be common (52).

Materials and Methods

The ApRES System.

In this study, we deployed an ApRES system with 16 cavity-backed bowtie antennas (eight transmitting, eight receiving) at site S30 on Store Glacier (58). The ApRES recorded the relative depths of internal reflectors and the ice base in a Lagrangian reference frame by transmission of a frequency-modulated continuous wave. The signal frequency increases linearly from 200 MHz to 400 MHz over 1 s, corresponding to a (coarse) range resolution of 0.43 m before processing steps (29). Combined with phase measurements embedded within each coarse-range bin, range detection with millimeter precision can be achieved given ideal (low signal-to-noise ratios) conditions (29). The ApRES was deployed to run autonomously, collecting radar reflection data at 4-hourly intervals from 26 July to 4 December 2014, after which the antennas were damaged during strong winds, and data collection ceased. The dataset was prepared for calculation of vertical deformation and BMRs following well-established phase processing procedures (29, 32, 34, 59), as summarized below.

Vertical Deformation Rates.

Every 4 h, the ApRES system transmitted a burst for each antenna pair, giving a total of 64 chirps per burst (58). The high precision of the system allows the range to englacial reflectors and the ice–bed interface to be resolved at every transmitted burst (). Vertical ice deformation rates and their respective errors were derived for burst pairs separated by 24 h as shown in and described in detail by Young et al. (32). For the ApRES data acquired from 26 July to 4 December 2014, we found internal layers’ vertical displacements to increase linearly through the ice column (), indicating that the vertical strain rate is depth independent. We note that only internal layers with a cross-correlation coherence threshold of between consecutive measurements were used to estimate the vertical deformation (32). With this threshold, vertical deformation rates were based on linear displacements of internal layers observed in this upper 80% of the ice column (). Although we cannot rule out the possibility of a different deformational regime occurring within the lowermost 20% of the ice column, we assume that the strain remained depth independent throughout the ice column, while noting that englacial deformation at depth will typically be either a continuation of the upper (overlying) strain regime or an enhanced expression thereof (32). If the latter were the case, basal melting would be higher than our radar-derived estimates, because the ice column is observed to thicken at site S30 (). For the data we present, a robust linear regression model () was the best fit. See Young et al. (32) for details.

BMRs.

The BMR was derived from ApRES data using the same approach adopted in studies of Antarctic ice shelves (33, 34). By differencing the total amount of vertical deformation, occurring over the ice column between consecutive measurements, with the concurrent change in ice thickness, we generated a 4 month long time series of daily BMRs beneath Store Glacier (). We estimated the total change in ice thickness at each time step by determining the coarse-range offset of the bin enveloping the basal reflector through identifying the amount of lag corresponding to the maximum amplitude of the cross-correlation (), and its respective fine-range offset through the phase of the complex cross-correlation function (). The SE of ice thickness measurements was derived from phase variations across all chirps in each burst, as described by Young et al. (32). The BMR error bounds were calculated as the square root of the sum of the squared errors tied to vertical deformation rate and ice thickness, thereby taking error propagation into account. The location of the basal reflector was defined as the first range bin in which all bursts coherently capture a single reflector from the dielectric contrast between ice and liquid water at the bed. The shape of the basal reflector within the study area is dictated by the topography of the underlying bed layer, and, therefore, on a fast-flowing glacier, the aperture footprint of the ApRES array would be incrementally offset due to the down-glacier movement of the ice surface. Range change was therefore conducted incrementally at a time step of 24 h to not only to minimize the change in basal reflector shape but also enable tracking the basal reflection peak through range and phase. As the radar system is resolution limited with a footprint of radius , where R is the total range from source to reflector and is the coarse-range resolution (0.43 m); this corresponds to a maximum daily offset in footprint radius and area of 1.8 m (8%) and 167 m2 (10%), respectively, given the observed maximum surface velocity of 672 m ⋅y−1 (). Because this offset is minor, we assume no influence on the ApRES measurements. However, we conservatively apply an additional 10% of the daily measured BMR values to their respective error bounds in a first attempt to account for these unknowns.

Basal Ice and Borehole Temperature Records.

At the ice–bed interface, the melting point temperature of ice, adjusted for pressure, varies according to the Clausius–Clapeyron gradient, which is 7.42 × 10−8 ∘C ⋅MPa−1 for pure ice and air-free water (). With an ice density of 917 kg ⋅m−3 and a pure ice column with no firn present, we estimate the pressure-adjusted melting point temperature beneath a nominal 604 m to 606 m of ice to be –0.40 ∘C. This melting point temperature is substantially lower than the T1 temperature sensor record used to infer viscous heat dissipation in the basal drainage system. While we cannot rule out the possibility of error in our measurements, Doyle et al. (30) present three lines of evidence that indicate that the T1 temperature sensor was calibrated and operational: 1) The thermistor ice bath calibration curve for T1 was consistent with those of all the other thermistors; 2) the temperature time series for T1 does not show the characteristic freezing curve observed for all the other thermistors, which suggests that the thermistor did not freeze in; and 3) damage to the thermistor cable caused by deformation or basal sliding would be likely to stretch the cables, which would increase its resistance and drive apparent temperature downward, not upward. The observation of T1 peak temperatures in unison with a spike in electrical conductivity shortly after the late August rainfall event, which brought warm air and precipitation over the ice sheet, also indicates that sensor T1 was working (30). During this event, surface ablation was measured at the seasonal peak rate of 56 mm ⋅d−1 (30), hence indicating that the coincident T1 peak was induced by viscous dissipation in an expanded basal drainage system carrying a large volume of surface meltwater. Hence, we infer the T1 record to capture the effect of water warmed by viscous heat dissipation in the basal drainage system. The T1 record is in good agreement with, and also independent of, the radar-derived BMR.

Basal Heat Budget.

To understand the high BMRs observed beneath Store Glacier, we quantified sources and sinks of heat at the base of the ice sheet,where G is geothermal heat flux; the second term is frictional heat calculated from basal shear stress τ and basal motion U; the third term is the conductive heat loss calculated from the basal ice temperature gradient θ and ice thermal conductivity K; the fourth term is latent heat of fusion calculated from the BMR (negative when ice base freezes); L is the coefficient of latent heat of fusion, and ρ is the density of ice. Thee four terms define the standard basal heat budget used in most previous work (60) and are here used in our initial heat budget calculation. The fifth term, Q, is added in order to also include energy released due to the viscous resistance in the water flow beneath the ice (37). Below, we describe how each term was quantified.

Geothermal Heat Flux.

In this study, we used a geothermal heat flux value of 60 × 10−3 W ⋅m−2 based on crustal magnetic field (61) and thickness (35). Although modeled geothermal heat flux over Greenland is highly variable, ranging from 40 × 10−3 W ⋅m−2 in the south to 140 × 10−3 W ⋅m−2 in the central north where the crust is thinner, this variability is not important in this study, as other basal heat sources are significantly higher.

Frictional Heat.

To quantify the frictional heat, we derived estimates of the basal shear stress, τ, and the rate of basal motion, U, from contemporaneous observations in colocated boreholes. Because the glacier is underlain by till (31), we used the Coulomb plastic failure criterion to describe the till’s shear strength, that is, , where is the effective pressure calculated as the difference between ice overburden pressure (, where h is the ice thickness measured by ApRES) and the basal water pressure (p, recorded by borehole pressure sensor). Due to precise measurements of both h and p (Fig. 2), we were able to quantify the effective pressure very accurately. The characteristic friction constant, , does not vary greatly across different glacial environments, and we assumed a value of 30∘ which is shared by most normally consolidated tills (62). Following Ryser et al. (63), the rate of basal motion was derived by subtracting tilt recorded in a borehole drilled to the bed from a contemporaneous GPS record of surface motion at the borehole site (Fig. 2). The tilt sensors were processed assuming the produced vertical gradients of horizontal velocity were all in the flow direction. The resulting time series was filtered with a two-pole, low-pass Butterworth filter with a 72-h cutoff period, and then binned into daily averages to match the time steps of other parameters. More detailed descriptions of these borehole records can be found in Doyle et al. (30).

Conductive Heat Loss.

The conductive heat loss of –60 mW ⋅m−2 was derived from the thermal conductivity of ice () and a basal ice temperature gradient of –0.0286 ∘C ⋅m−1 established from borehole temperature records shown in Fig. 2. The equilibrium temperatures were –0.86 ∘C for sensor T3 (installed at 596.5 m below surface), –0.76 ∘C for sensor T2 (600.5 m), and –0.64 ∘C for sensor M1 (603.3 m). Sensor T1 did not freeze in. Details of these records can be found in Doyle et al. (30).

Viscous Heat Dissipation.

When surface meltwater is injected to the bed, energy (Q in Eq. , ) is released due to the viscous resistance in the water flow. We partitioned this energy into gravitational and potential energy components, using the approach described by Mankoff and Tulaczyk (37).

Between the injection point and outflow, we assume all energy is dissipated as heat within the grid cell where the energy transfer occurs. Henceforth, as water flows down the hydraulic gradient, we calculate the energy released as heat based on the volume of water that is routed in each grid cell, including 1) the change in the hydraulic potential and 2) the change in the pressure-dependent phase transition temperature (37),where V is volume of water; is the hydraulic potential gradient, where the subscript p denotes the pressure component; C is the Clausius–Clapeyron gradient; c is the specific heat of water, and ρ is the density of water (see for parameter values).

Water Routing.

To quantity Q, the amount of energy available for viscous heat dissipation at site S30, we used a hydrological model in which water is routed subglacially in the catchment beneath Store Glacier. Specifically, the model tracks the flux of surface meltwater from source (i.e., surface runoff reaching the bed) to sink (i.e., subglacial discharge into fjord), in order to estimate the energy produced by pressure and elevation changes. The energy for viscous heat dissipation in the basal drainage system was estimated using daily values of surface runoff from the RACMO2.0 regional climate model (39) under the assumption that all surface water reaches the bed and that all energy is dissipated as heat (37). To route water, we used r.watershed tool in GRASS GIS as a directional routing algorithm in which cells with lower hydraulic potential receive a fraction of the outflow (37). The hydraulic potential was calculated as (64),where ρ is the density of water, g is gravitational acceleration, and z is the bed elevation, prescribed from BedMachine 3.0 topographic data (38); α is the flotation fraction, here set to 0.9 based on ice overburden pressure from measured ice thickness and basal water pressure (Fig. 2); ρ is the density of ice, and z is the surface elevation as prescribed by ArcticDEM. The resulting model output was gridded at a 500-m spatial resolution. Runoff was injected at the bed beneath the grid cell in which it was produced, and the water was assumed to be at the subglacial pressure-dependent phase transition temperature; that is, we ignore any warming at the surface from radiative sources while assuming that the water cools according to pressure change between the bed and the surface. The energy for viscous heat dissipation in our model occurs when there is a drop in either gravitational potential energy (first term on the right-hand side [RHS] in Eq. ) or pressure (second term). When water flows under thinning ice where the phase transition temperature increases, energy is used to warm the water, resulting in either less melting or a switch to basal freezing if the drop in gravitational potential energy cannot provide sufficient heat. Basal freezing may also occur if water flows uphill and the pressure drop cannot provide sufficient energy, whereas viscous heat dissipation will melt ice the fastest when there is a drop in gravitational potential and an increase in pressure.

Enthalpy of Basal Ice.

To supplement the thermomechanical model (Eq. ), we also calculated basal melting under the assumption that the basal ice is at the phase change temperature. In this case, we used a one-dimensional representation of the jump equation for enthalpy to derive a BMR (9),

The first and second terms on the RHS of Eq. are the frictional heat and geothermal heat flux (described above). The last term describes how changes in subglacial water pressure, , are related to energy fluxes when the ice is underlain by a subglacial water layer with enthalpy, , and thickness, η. The third term, , is the nonadvective heat flux into temperate basal ice expressed in terms of pressure (p) and enthalpy (H),where is the thermal conductivity of the temperate ice–water mixture, with k for pure ice and k for liquid water, ω is the water fraction, and K0 is temperate ice diffusivity. Fig. 3 shows the additional energy for basal melting, when melt rates from Eq. based on parameter values shown in are compared with those derived from Eq. .

Equilibrium Water Temperature.

Energy dissipation occurs inside glacial conduits due to viscous resistance in the flow. As the dissipating energy warms the water (see ), the heat loss into the conduit ice wall also grows, which gives an equilibrium condition when the two are equal. In a straight inclined conduit with stable water flow, the equilibrium temperature is (44),where R is hydraulic radii, s is the hydraulic slope based which combines the gradients of elevation and pressure head, and c is an empirical constant for turbulent flow at 0 ∘C (see for value). Eq. shows that will always be higher than Tm and that water temperature of 0.88 ∘C can be reached when T is –0.40 ∘C, for example, if the gradients of elevation and pressure head are 5∘ each and R is 2 m.

Heat Transfer.

We estimated the heat transfer to be 60 W ⋅m−2 ⋅∘C−1 to 170 W ⋅m−2 ⋅∘C−1 by assuming that the highest (lowest) observed BMRs of 57 mm ⋅d−1 (10 mm ⋅d−1) were driven by water temperatures measured at 0.88 ∘C (0.19 ∘C). To make a first order estimate of the associated water flow rate, v, we assumed a linear relationship with the heat transfer coefficient (44),where c is the constant for turbulent flow at 0 ∘C (). The observationally derived heat transfer of 60 W ⋅m−2 ⋅∘C−1 to 170 W ⋅m−2 ⋅∘C−1 can therefore be associated with theoretical flow velocities of 2.2 cm ⋅s−1 to 6.6 cm ⋅s−1.

To link the heat transfer with a first-order estimate of the water depth, D, we turned to an empirical relationship developed for heat transfer to a river ice cover (65),where B is an empirical constant (). This equation suggests that the heat transfer coefficient is relatively insensitive to the water depth and that the main control comes from flow velocity. A heat transfer coefficient of 60 W ⋅m−2 ⋅∘C−1 to 170 W ⋅m−2 ⋅∘C−1 can be achieved for a large range of water depths between 0.01 m and 10 m, while the water velocities should be on the order of 1 cm ⋅s−1 to 10 cm ⋅s−1.