Research ArticleGEOPHYSICS

Thermodynamics of a fast-moving Greenlandic outlet glacier revealed by fiber-optic distributed temperature sensing

See allHide authors and affiliations

Science Advances  14 May 2021:
Vol. 7, no. 20, eabe7136
DOI: 10.1126/sciadv.abe7136


Measurements of ice temperature provide crucial constraints on ice viscosity and the thermodynamic processes occurring within a glacier. However, such measurements are presently limited by a small number of relatively coarse-spatial-resolution borehole records, especially for ice sheets. Here, we advance our understanding of glacier thermodynamics with an exceptionally high-vertical-resolution (~0.65 m), distributed-fiber-optic temperature-sensing profile from a 1043-m borehole drilled to the base of Sermeq Kujalleq (Store Glacier), Greenland. We report substantial but isolated strain heating within interglacial-phase ice at 208 to 242 m depth together with strongly heterogeneous ice deformation in glacial-phase ice below 889 m. We also observe a high-strain interface between glacial- and interglacial-phase ice and a 73-m-thick temperate basal layer, interpreted as locally formed and important for the glacier’s fast motion. These findings demonstrate notable spatial heterogeneity, both vertically and at the catchment scale, in the conditions facilitating the fast motion of marine-terminating glaciers in Greenland.


Mass loss from the Greenland Ice Sheet (GrIS) has increased sixfold since the 1980s and is now the single largest cryospheric contributor to global sea-level rise (1, 2). Around half of this mass loss is from surface meltwater runoff, while the other half is driven by discharge of ice into the ocean (3). Surface melting and runoff can be observed by satellites and in situ methods in a relatively direct manner, facilitating validation of surface process models and a more straightforward exchange between modeling and empirical observations [e.g., (4)]. In contrast, ice dynamics, which are largely controlled by englacial ice properties and basal conditions, present substantial observational challenges and are therefore less well constrained in models, making such processes one of the largest sources of uncertainty in global sea-level projections (5, 6).

A crucial but scarce source of information comes from boreholes, which provide unambiguous records of subsurface ice properties, such as temperature profiles, deformation rates, and the conditions at the bed that influence basal motion. These records are especially relevant for Greenland’s marine-terminating glaciers, which flow rapidly (≥200 m a−1) through englacial deformation and basal motion (7). Such glaciers drain most of the ice sheet (2, 8), yet they are rarely the subject of borehole investigations due to logistical challenges including pervasive crevassing and high strain rates that damage installed equipment [e.g. (7, 9, 10)]. As a consequence, GrIS models often rely on physical constraints from theoretical considerations or from studies conducted on smaller and more accessible glaciers elsewhere, which may not be representative of larger or faster ice masses [e.g. (11, 12)].

Here, we report results from the first use of fiber-optic distributed temperature sensing (DTS) on the GrIS, installed in a 1043-m-deep borehole drilled to the base of the fast-moving Sermeq Kujalleq (Store Glacier in the literature and used hereafter), which has one of the largest glacial catchments in Greenland [34,000 km2, (13)]. The installed fiber was sampled every 0.25 m, providing an englacial borehole temperature record with a vertical resolution (0.65 m; see Materials and Methods for further information) two orders of magnitude higher than in previous comparable borehole studies using discrete sensor setups. With this highly resolved temperature record, we report on the thermodynamics of Holocene ice [formed since 11.7 thousand years (ka) in the current interglacial period], pre-Holocene ice (formed before 11.7 ka in the most recent glacial period and prior), and the Last Glacial-Interglacial Transition that separates the two sections. We also report on the transition from cold ice (with no liquid phase) to basal temperate ice (in equilibrium with a liquid phase at the ice-water phase transition temperature), all in unprecedented detail. These results indicate a heterogeneity in deformation between and within different ice sections that emphasizes a need for a multilayered approach to ice rheology and represent an important improvement in our observation-based understanding of ice deformation in fast-moving sectors of the GrIS.


We installed a DTS system on 5 July 2019 in a 1043-m hot-water drilled borehole located at site R30 on Store Glacier (70.57°N, 50.09°W; Fig. 1), where surface ice motion is ~600 m a−1. Borehole BH19c was drilled at the center of Lake 028, which had drained rapidly 2 months before drilling (14), leaving a large fracture and active moulin with continued inflow by supraglacial melt water (Fig. 1). The DTS system, which consisted of a Silixa XT-DTS, armored cable housing fiber optics, and thermistors for calibration purposes, works by measuring the temperature-sensitive components of Raman backscatter (Stokes and anti-Stokes wavelengths) and travel time, when a laser pulse is transmitted through the optical fiber (Materials and Methods).

Fig. 1 Map showing Store Glacier in West Greenland.

(A) Map showing flow of the GrIS and glaciers flowing into Uummannaq Fjord. Black lines show direction of ice flow (originating 80 km inland at 2-km spacing). Surface velocity (color scale) is 2018 annual velocity from MEaSUREs data (69, 70). Off-ice surface elevation (gray scale) and 1000-m ice-sheet surface elevation contour (green line) are from ArcticDEM v3 (71). Yellow box in Greenland inset shows study location. (B) Aerial image of site R30 showing location of boreholes BH19a-d (including BH19c where the DTS system was installed) and BH19e-g (red dots). Moulins (green dots) are fed by supraglacial streams. Black dashed line traces fracture that caused supraglacial lake drainage and moulin formation before boreholes were drilled. Image acquired by drone on 21 July 2019 (14).

Our DTS measurements span from 5 July to eventual cable failure at 889 m depth at 21:30 UTC on 13 August (see Materials and Methods for sampling time information). Although our observations are limited to this 40-day period, the temperature profile for measurements averaged over the last 4 days before cable failure is within 0.2°C of equilibrium ice temperatures calculated from theoretical cooling curves at all points above 925 m (Materials and Methods). Below, we describe our observed vertical temperature profile (Fig. 2) in three distinct sections: first, cold ice with smooth temperature gradient variations interrupted by isolated perturbations; second, cold ice characterized by a shallower temperature gradient, more frequent temperature gradient variations, and substantial heterogeneous cable deformation; and third, temperate basal ice.

Fig. 2 Vertical ice temperature profiles.

(A) Full temperature record from DTS measurements averaged over the 96 hours prior to the last recorded measurement before cable failure at 21:30 UTC on 13 August. The combined sampling time was 3.8 hours, with 10 min of active recording every 4 hours over this period. Solid black line is recorded temperature, with 95% confidence interval shown in light gray shading. Red line is the equilibrium temperature estimated from observations and theoretical freezing curves, with dark gray shading showing root mean square error. Horizontal gray dashed line (at 889 m depth) is the point of cable failure on 13 August and the inferred location of the Last Glacial-Interglacial Transition. (B) Close-up of Anomaly-208 with same axis units. (C) Temperature gradient, with orange circle highlighting a temperature gradient anomaly at 100 to 111 m. (D) Close-up showing temperatures in the bottom part of the borehole below 880 m. Orange line is the pressure-dependent melting point assuming a linear Clausius-Clapeyron slope of 9.14 × 10−8 K Pa−1. Black circles are thermistor data. The highest thermistor at 1033 m is interpreted to be an outlier. (E) Inset shows temperatures in the lowermost 100 m.

First, the temperature decreases progressively from −1.2°C at the surface to a minimum of −21.1°C at 580 m depth, where the temperature gradient switches from negative to positive downward, reaching a maximum sustained (varying smoothly over >20 m) value of 0.13°C m−1 at 865 m. This concave profile is a characteristic feature of fast-moving GrIS outlet glaciers (15), where advecting ice from the cold, high-elevation interior has not yet been warmed by heat diffusion from the warmer surface air temperature at lower elevations, basal geothermal heat flux, and strain heating concentrated toward the bed (16). This section terminates at the location of cable failure (889 m), where the temperature is −5.3°C following a slight decrease in temperature gradient to 0.11°C m−1. Herein, we observe (i) a local minimum in temperature of −6.9°C at 8 m depth, followed by a local maximum of −6.7°C at 16 m depth; (ii) two localized gradient increases and reversals at 100 to 111 m and 208 to 242 m (hereafter Anomaly-208); (iii) a sharp (4 m) and temporary 60% dip in temperature gradient of 0.07°C m−1 at 837 m, below which is a minor (10%) increase in temperature gradient of 0.1°C m−1. Anomaly-208 is the largest of these anomalies (0.25°C) and is clearly distinguishable throughout the entire time series (Fig. 3), indicating that it is a stable feature over the length of observation.

Fig. 3 Ice temperature variations during the measurement period.

(A) Full temperature record from 5 July to 13 August. The red dashed line is the boundary of close-ups shown in (B) and (C). White dashed lines are temperature contours with 5°C spacing. (B) Close-up showing the temperate zone below 970 m. (C) Difference between temperature measurements shown in (B) and theoretical ice-water phase transition temperatures calculated with the Clausius-Clapeyron slope shown in Fig. 2 (9.14 × 10−8 K Pa−1). The bottom of the cold-temperate transition zone at the end of the measurement period (982 m) is shown as a dashed black line throughout in (B) and (C). The time series was created by averaging DTS measurements over 8-hour periods regardless of active DTS sampling time, giving 8 hours of active sampling time between 5 and 21 July, 96 min between 21 and 23 July, and 20 min beyond 23 July (with 1 in 24 outputs from a 10-min sampling period; see Materials and Methods for further information).

The top of the second section is defined by a step change in temperature gradient from 0.13°C m−1 to 0.07°C m−1 (53%) at 889 m (Fig. 2, A and B), below which the temperature continues to increase gradually downward, although with more fluctuations compared to the section above. We observe a notable temperature anomaly with magnitude 0.25°C at 912 m, associated with a temperature gradient fluctuation of magnitude 0.15°C m−1 over 5 m, as well as numerous temperature fluctuations between 889 and 935 m with a maximum amplitude of 0.1°C m−1 over distances as short as 0.6 m. Because this length scale approaches the DTS spatial resolution, we cannot rule out the possibility that the latter, smaller temperature perturbations, which occur after 8 August, stem from notable localized (<1 m) cable deformation (see Materials and Methods and figs. S1 and S2 for further details). Below 935 to 959 m, the temperature gradient oscillations reduce to the same magnitude observed in the uppermost section and the temperature gradient gradually settles below 959 m to a consistently low negative value (−0.002°C m−1) at 982 m. We use this gradual transition in gradient over 23 m to define a cold-temperate transition zone (959 to 982 m) as opposed to the cold-temperate transition surface used in previous theoretical studies [e.g., (17)].

Third, while we define the cold-temperate transition zone as a zone extending to 982 m depth, we observe temperatures within the boundaries of the expected pressure and solute-dependent ice-water phase transition temperature from 970 m downward (fig. S3). Hence, we define 970 m as the top of the basal temperate zone giving it a thickness of 73 m, noting that the upper 12 m also forms part of the cold-temperate transition zone. Within this lowermost section, we also observe spatiotemporal temperature changes (Figs. 2E and 3C), indicating that processes influencing the ice-water phase transition temperature are operating at a time scale of days to weeks. These subtle variations are especially notable within the lowermost 20 m of the ice column. The temperature at 1040 m drops 0.04°C to −0.85°C from installation on 5 to 14 July, increasing gradually to −0.83°C over the remainder of the record. While our recorded basal temperate zone thickness decreases significantly as the borehole refreezes in the days after installation (Fig. 3), this adjustment becomes gradually slower and ceases altogether on 10 August when our measurements stabilize. We take the stabilized record after 10 August (Fig. 2) to represent the true basal temperate zone thickness of 73 m.

Fig. 4 Integrated differential attenuation along the borehole averaged over two time periods.

(A and D) Averaged integrated differential attenuation profiles recorded between 21 and 24 July (10.5 hours total sampling time) and the 96 hours before cable failure at 21:30 UTC on 13 August (3.8 hours total sampling time), respectively. Close-ups are shown in (B), (C), (E), and (F). The orange circle in (D) shows a slight perturbation that matches the location of the temperature anomaly also identified by an orange circle in Fig. 2C.

The raw DTS data also give a qualitative indication of cable deformation. Integrated differential attenuation (Fig. 4) quantifies the difference in cumulative signal attenuation between the Stokes and anti-Stokes wavelengths along the fiber and is used in calibration to correct for optical perturbations that may arise from fiber deformation through bending or extension [(18), Materials and Methods]. Laboratory strain gauge testing, using the same cable as deployed in the field, showed that an integrated differential attenuation signal is first clearly visible at a strain of 1.00% and increases in magnitude up to a maximum tested strain value of 1.85% (figs. S4 to S6). We see a growing integrated differential attenuation signature of deformation throughout the measurement period in the exact location of Anomaly-208 (Fig. 4, B and E). We also see substantial variations between 889 and 935 m that vary on a meter scale (Fig. 4, C and F), suggesting more heterogeneous deformation than over Anomaly-208 where integrated differential attenuation varies smoothly. Below the cold-temperate transition zone, integrated differential attenuation remains as smooth as in the central section of the profile, suggesting that the cable has undergone less deformation there than in the region directly above.


The Store Glacier temperature profile at R30 (Fig. 2A) is broadly consistent with theoretical estimates of ice temperature governed by advection—which in glaciers such as Store brings cold ice from high elevations to warmer settings at lower elevations—and vertically directed diffusion (15). The advective component explains why ice temperature at 580 m depth is at −21.1°C, even though near-surface ice is much warmer given a mean annual temperature of −9°C (7). Variations in ice temperatures within ~10 to 15 m of the surface are tied to seasonal air temperature cycles (16). Temperatures higher than the mean annual air temperature below this depth can be explained by latent heat of fusion released when water freezes in firn pores and within crevasses in the ablation area (19), which can extend to depths of several hundred meters (20). Below 580 m, the temperature increases toward the basal zone, a result of geothermal heat flux and strain dissipation heating concentrated near the bed [e.g., (7)]. Where the bed is wet across large regions of the GrIS’s margins (21), and basal ice is therefore pinned at the phase transition temperature, frictional heat from basal sliding and sediment deformation and viscous dissipation of subglacial water also supply energy (22), with freeze-on of temperate ice occurring if the pressure-dependent melting point drops rapidly (23). Under such conditions, a basal temperate zone consisting of ice at the phase transition temperature may ultimately develop if strain heating exceeds the conductive heat flux away from the cold-temperate transition zone [e.g., (17, 24)] or if liquid water in veins or fractures freezes and releases latent heat within the ice (24). While our observations are consistent with previous records of these processes at other sites using discrete sensors, typically spaced tens of meters apart or more [e.g., (2527)], our continuous DTS measurements enable us to resolve both previously identified and presently undocumented thermodynamic processes in unsurpassed detail.

We separate the ice profile at R30 into three major and distinct sections, inferring (i) cold Holocene ice (0 to 889 m) characterized by smooth temperature changes and featuring the steepest sustained gradient, indicative of low and uniform deformation with isolated strain bands; (ii) cold pre-Holocene ice (889 to 970 m) characterized by a transition to a shallower, nonuniform temperature gradient and substantial variations in integrated differential attenuation, indicative of strongly heterogeneous deformation throughout; and (iii) temperate but still pre-Holocene ice (970 to 1043 m) with subtle departures from a linear Clausius-Clapeyron slope and with relatively low cable deformation. These sections and their implications are discussed in turn below.

Localized temperature anomalies in the upper Holocene ice

Although Anomaly-208 is unique in resolution, local temperature anomalies that could feasibly be fitted with a shape similar to Anomaly-208 have been detected in previous borehole-based studies, albeit at a lower resolution (7, 19, 26, 27). This suggests that anomalies such as the one described here may be widespread. We consider it unlikely that Anomaly-208 is a result of the thermal influence from drilling given (i) drilling rates were fast (~1.2 m min−1) and designed to produce a borehole with a uniform diameter of approximately 12 cm using a specialized model (fig. S7) through solid ice with no firn present, (ii) borehole drilling rates for proximal boreholes also varied smoothly over this depth range (fig. S7), and (iii) additional heat input in this depth range during drilling would increase the borehole diameter and lengthen the time until borehole freezing, which was not observed (Fig. 3).

Anomaly-208 may then hypothetically be a result of a latent heat source, historic deposition of snow at elevated temperatures, or enhanced deformation. A latent heat source from freezing of a nearby water-filled crevasse would, however, lead to a thermal influence over a much greater vertical extent than we observed, as recorded at the terminus of Bowdoin Glacier further north (10), and would not lead to the sustained cable deformation we report (Fig. 4, B and E). In addition, water freezing in crevasses or conduits cannot explain the heat sink above the anomaly, which is required to maintain its shape on multiyear time scale (Fig. 5). Deposition of snow at elevated temperatures in a past climatic event can also be ruled out as this would produce a far more diffuse profile (28). We therefore hypothesize that strain heating can explain the temperature anomaly, with additional evidence provided by the integrated differential attenuation recorded in the optical fiber at exactly the same depth (Fig. 4, B and E). The integrated differential attenuation record shows that the cable froze in with no measurable deformation, but increasing strain was placed on the cable at the exact location of Anomaly-208 over the duration of the measurement period.

Fig. 5 Temperature anomaly reproduced in a simple strain band model.

(A) Ice temperature distribution of Anomaly-208, with dotted black line showing the final observed temperature profile. Solid black line is a theoretical temperature profile after 400 days of vertical diffusion and no horizontal advection. Dashed yellow line is temperature profile obtained with a strain band model with square-wave heat source and sink as shown in (B). (B) Energy flux input and sink used in the model. (C) Schematic representation of strain band (darkening blue band) with (a) possible deep fracture, (b) advecting cold ice, and (c) concentrated deformation.

In a simple vertical one-dimensional (1D) heat diffusion model, Anomaly-208 rapidly (~1 year) loses much of its shape and magnitude if there are no heat inputs or sinks (Fig. 5) but is maintained when a square-wave heat input of 30.8 mW m−2 is added between 222 and 227.5 m together with a heat sink of 10.8 mW m−2 between 208 and 216 m (Materials and Methods). This energy input and sink can be used to estimate the associated deformation rate. If ice is assumed to deform in simple shear, we obtain a deformation velocity integrated over the shear band thickness of 14 m a−1 (see Materials and Methods for calculation and creep parameter uncertainty). When the heat sink is assumed to be the result of an advective flux with constant velocity through depth, we obtain a deformation velocity of 8 m a−1 (Materials and Methods). Although, once installed, the borehole cable moves with the ice (a Lagrangian reference frame), the influence of ice advection and deformation that occurred prior to installation (a Eulerian reference frame) is recorded by the installed cable, as the magnitude of subsequent changes relative to the inherited signal is small. The advective flux described here is assumed to be composed of incoming cold ice that has not yet been thermally affected by localized strain heating, similar to the interpretation of Hills et al. (29), who argued that a negative change in ice temperature in a setting where vertically directed heat diffusion and strain heating would be expected to increase temperature must result from advection. This requirement for an advective heat sink implies that heat production and Anomaly-208 are spatially localized features. A deformation rate of 8 to 14 m a−1 accounts for roughly 2% of surface velocity (600 m a−1 at R30), or 7% of internal ice deformation if internal deformation and basal motion are partitioned similarly to previous borehole observations at site S30, located 8.9 km away (Fig. 1) (7).

Enhanced deformation could result from stress perturbations, local rheology changes, or both. Volcanic ash would reduce viscosity, in the same manner that wind-blown dust softens glacial ice (30), and volcanic events are well recorded in GrIS ice cores (31) over the radio stratigraphy age estimate for Anomaly-208 [0 to 3 ka, (32)]. However, the relative smoothness of Anomaly-208 over a depth range of 30 m, which contrasts with temperature and integrated differential attenuation perturbations within pre-Holocene ice (discussed below), suggests that it is not caused by one or a series of discrete thin horizons. As crevasses prevent non-orthogonal stress transfer (33) and water-filled crevasses were observed to a maximum depth of 265 m in a borehole drilled 1 km away at site R29 (Fig. 1) (20), we instead hypothesize that Anomaly-208 is a result of increased deviatoric stress below the depth of maximum crevasse propagation or a result of longitudinal and lateral stress transfer from spatially varying bed conditions [e.g., (34)]. With similar integrated differential attenuation signatures tied to other temperature anomalies at 100 to 110 and 837 m depth, we further hypothesize that strain heating can occur within Holocene ice as well as within pre-Holocene ice in which variable deformation is more commonly observed [e.g., (7, 25, 26)].

Enhanced deformation in pre-Holocene ice

Pre-Holocene ice is typically identified geophysically by lower reflectivity in radio-echo sounding data (35) or an anisotropic reflection in active seismic data (36). Pre-Holocene ice has a higher dust and impurity content than Holocene ice, impeding grain growth such that crystals remain smaller in the former. It also has a strongly anisotropic rheology (37), with a near-vertical c-axis orientation that develops from bed-parallel shear (38), meaning the strain rate of this ice is sensitive to the orientation of the applied stress field. Consequently, deformation rates are typically 2.5 to 5 times greater in pre-Holocene ice than in Holocene ice under similar temperature and stress conditions (30, 38, 39). The inferred transition from Holocene to pre-Holocene ice (i.e., the Last Glacial-Interglacial Transition) at 889 m depth (85% of the ice thickness) matches the transition from isotropic to anisotropic ice recorded at ~880 m depth in the same borehole using distributed acoustic sensing with 10-m spatial resolution (40). The inferred Last Glacial-Interglacial Transition also agrees well with the depth inferred previously from seismic studies of this transition at site S30 on Store Glacier [84%, (36)], and it falls within the estimated range for the Last Glacial-Interglacial Transition in a wider regional study [82 to 85%, (35)]. While pre-Holocene ice properties are well studied, we are the first to report the thermal and physical manifestation of these properties, and those of the Last Glacial-Interglacial Transition, on a glacier in a marine-terminating drainage basin.

The drop in temperature gradient (53%) across the Last Glacial-Interglacial Transition is consistent with increased strain heating through the pre-Holocene ice, while the larger-scale variations (>5 m) in temperature gradient (e.g., 915 m depth) are consistent with enhanced strain heating similar to that associated with Anomaly-208, discussed above. If these variations are a result of deformational folds (41), these folds would need to occur on a small scale (~5 m). The highly heterogeneous integrated differential attenuation response (Fig. 4) further suggests that enhanced strain may occur at a submeter scale, consistent with deformation trends determined from ice core records at the North Greenland Eemian Ice Drilling (NEEM) site for pre-Holocene ice at premelting temperatures [>−10°C, (38)]. The location of cable failure falls exactly at this step change in temperature gradient, suggesting that the Last Glacial-Interglacial Transition is a high-strain interface. This is expected as the less viscous Holocene ice would be more capable of transferring stress, as hypothesized at the ice stream margin of Jakobshavn Isbræ (42), in slower flowing ice ~30 km further north of Jakobshavn Isbræ due to changes in the subglacial hydrological system at the end of the melt season (34) and in much smaller valley glaciers [e.g., (43)]. Our results indicate that this behavior also occurs at Store Glacier.

Heterogeneous ice deformation and fabric in pre-Holocene ice are extensively studied at ice core sites at interior locations near ice divides, where ice flow is slow [e.g., (38, 39)]. Our results, specifically the integrated differential attenuation, indicate that meter- or submeter-scale deformation heterogeneity in pre-Holocene ice may also be important for the fast motion of marine-terminating glaciers such as Store. This is important as this variability is not well represented in the widely applied Nye-Glen isotropic flow law (44, 45) and as pre-Holocene ice is generally interpreted to be responsible for most of the ice deformation occurring in fast-moving outlet glaciers [e.g., (7, 25)]. The reason that such heterogeneous ice deformation has not previously been observed on glaciers such as Store may stem from strain being estimated by discrete sensors in cabled borehole studies, which may not have sufficiently high spatial resolution to fully capture the spatial heterogeneity [e.g., (7, 26)], with the exception of Ryser et al. (25) who captured this heterogeneous behavior within the limits of 40-m sensor spacing. The more limited spatial resolution in discrete inclinometer studies compared to DTS may mean that strain-banding location has been missed or overlooked in previous studies, possibly leading to an overestimation of basal motion.

The cold-temperate transition zone itself is often described as a discrete, isothermal surface in ice-sheet models [e.g., (17)] and inferred as such between discrete temperature sensors [e.g., (7, 26)]. However, we observe its presence as a change in temperature gradient over 23 m, which supports Lliboutry (46, 47) who proposed a “fuzzier” or diffuse transition due to changes in pressure in veins and lenses and a dependence of melting temperature on the scale of ice crystals. While 23 m (959 to 982 m) falls below the resolution of most models, a “fuzzy” cold-temperate transition zone, as observed here, may have important implications for hypothesized processes that shift the cold-temperate transition zone position, such as water movement through the ice matrix (24) or basal crevassing (48).

Temperate ice beneath the cold-temperate transition zone

The basal temperate ice recorded at R30 is substantially thicker than that at site S30, located only 8.9 km away [Fig. 1; (7)], where thermistors close to the bed constrain temperate ice thickness to at most a few meters. Furthermore, boreholes drilled only 1 km away at site R29 in 2018, where ice thickness is 950 m, limit the basal temperate zone to a maximum of 14 m there (20). This shows that the vertical extent of temperate ice can vary greatly over distances as short as one ice thickness. Investigations of Jakobshavn Isbræ (9, 26, 49, 50) also reveal a markedly varying temperate zone thickness, with as much as 300 m or ~10% of the total ice thickness inferred as basal temperate ice near the center of the ice stream, although the borehole at that location did not reach the bed. However, Funk et al. (49) interpret this to arise from convergent flow into a bedrock trough that is 1.5 km below sea level, whereas our study region exhibits lower convergence (Fig. 1) and bedrock topography amplitude is <500 m (fig. S8). This distinction between Jakobshavn Isbræ and Store Glacier implies that basal temperate zone thickness may be attributable not only to convergent flow into a bedrock trough but also to spatially varying basal conditions.

Because diffusive heat flux is directed away from the cold-temperate transition zone in both vertical directions (precluding conductive transfer of basal heat inputs into cold ice above the cold-temperate transition zone), a continuous heat source is needed to maintain it. Potential heat sources are advective heat flux, strain heating, and latent heat of freezing from movement of water through grain boundaries or fractures (24). Basal melting from geothermal heat flux, high viscous heat dissipation from turbulent melt water (22), and sliding across high-traction sticky spots may also reduce the thickness of the temperate zone. The variability in the thickness of the temperate zone ice we record in this sector of Store Glacier (~5 to 73 m) therefore indicates that the conditions responsible for temperate zone thickness are similarly variable over distances of ~1 to 10 ice thicknesses. We hypothesize that the greater temperate ice thickness at R30 compared to R29 may be a result of local variations in strain heating as the ice passes through a bedrock saddle (fig. S8) or encounters variations in subglacial sediment properties and drainage conditions. Booth et al. (40) show that borehole R30 at Store Glacier is likely to be underlain by consolidated sediment, the shear strength of which may vary depending on porosity and water supplied in subglacial drainage pathways (51), and Hofstede et al. (36) report notable differences in the availability of water at the ice-sediment interface beneath Store Glacier across distances of only a few kilometers. Alternatively, latent heat may be supplied to the cold-temperate transition zone via basal crevassing by surface water from the drainage of the overlying lake and subsequent moulin development (14). As temperate ice exhibits strain rates 5 to 10 times greater than cold ice under equivalent stress (24, 52), this may have substantial ramifications for ice flow modeling in terms of both deformation profiles and basal conditions. Our findings and previous studies at Jakobshavn Isbræ (9, 26, 49) demonstrate a clear need to advance our understanding of temperate zone processes in fast-moving glaciers.

The temperate basal zone does not have a constant temperature gradient predicted by a spatiotemporally constant Clausius-Clapeyron value, water content, mean stress, and solute concentration (Fig. 2E). Drawing on original studies on alpine glaciers [e.g., (53, 54)], we conclude that the temperature of temperate ice in Store Glacier varies over time scales of days and weeks. However, the temperate ice we observe is of “Canadian” type (55), where water content is derived from the melting of basal ice and possible inclusion of subglacial water, rather than primarily the metamorphosis of wet firn into ice observed at valley glaciers such as Blue Glacier, Alaska (54). Calculations of the influence of Clausius-Clapeyron value, shear stress, and solute concentration on temperate ice temperature (fig. S3) (16) show that the first two, and to a lesser degree, solute concentration, can explain the variation we observe. The variation may also result from incomplete return to the phase transition temperature following rapid decompression (56). While it is challenging to effectively isolate the primary cause of these variations, our high-resolution temperature record helps to constrain the parameter space of water, heat, and solute transport mechanisms within basal temperate ice (fig. S3), which may influence the cold-temperate transition zone and exert a direct control on ice deformation rates.

Last, we do not record a substantial integrated differential attenuation response within the basal temperate zone, which is consistent with low deformation rates within basal temperate ice in previous borehole studies (25, 26, 57). This reduced deformation may be a result of the larger crystal fabric associated with temperate ice inhibiting grain boundary sliding and hence limiting deformation (58) or alternatively be a result of concentrated high deformation rates over the cold-temperate transition zone in the same manner as over the Last Glacial-Interglacial Transition. The integrated differential attenuation response at the top of the cold-temperate transition zone (fig. S1) supports the latter interpretation, which is further supported by previous work showing a marked reduction in ice viscosity at temperatures close to, but below, the pressure-dependent melting point due to premelting effects [e.g., (38, 59)].

New insights from DTS

Our results demonstrate the transformative potential of high-resolution DTS in the field of glacier thermodynamics. We have shown that substantial heat can be generated by previously unrecorded localized strain within Holocene ice, which is otherwise characterized by a smooth variation in temperature gradient. We record fundamental differences in the temperature profiles of Holocene and pre-Holocene ice, demarcated by a step change in temperature gradient and a potential high strain interface at the Last Glacial-Interglacial Transition. The underlying pre-Holocene ice is characterized by temperature anomalies and substantial integrated differential attenuation variations, suggesting strongly heterogeneous ice deformation. We identify a diffuse 23-m transition from cold pre-Holocene ice to a basal temperate zone, with less cable deformation within the temperate ice than in overlying cold ice, suggesting that deformation may be concentrated in the cold-temperate transition zone as well as the Last Glacial-Interglacial Transition. Measurements at our study site show a basal temperate zone thickness of 73 m, considerably thicker than previous measurements of <14 m at a site just 1 km distant (20). As temperate ice has a viscosity 5 to 10 times lower than cold ice (24, 52), ice deformation may therefore vary greatly over distances as short as one ice thickness. These findings demonstrate that ice deformation and temperature, both vertically in borehole profiles and spatially within catchments, exhibit greater heterogeneity than previously considered.


Drilling and installation

Fiber-optic cables were installed in borehole BH19c, at site R30 within the drained lake L028, drilled on 4 to 5 July 2019 using the system of Doyle et al. (7) with an additional heater pressure unit to deliver a total of 60 liters min−1 at a surface temperature of roughly 80°C and a 1350-m hose. To reduce the thermal disturbance to background ice temperature and thus the time for the refrozen borehole to reach equilibrium and the overall drilling time, Greenler et al.’s (60) model was used to calculate optimum drilling speeds for a borehole diameter of 12 cm at a time 4 hours after installation using the ice temperature profile for S30 (Fig. 1) (7). An average drill speed of 1.2 m min−1, modulated to account for colder ice in the center of the profile, was used from 0 to 1020 m, decreasing to 0.5 m min−1 for the final 10 m due to loading and unloading behavior in the drilling system. This contrasts with previous boreholes on Store Glacier (7), which were drilled with a more constant speed, resulting in a variable borehole diameter that was in some places larger than necessary, and which therefore extended the time for borehole freezing. While faster speeds were theoretically possible with the supplied heat input, the narrow aperture of the drill nozzle meant that this was not mechanically feasible. Drilling ceased when the borehole rapidly drained. Figure S7 shows the modeled borehole diameter for each of the three boreholes drilled to the bed in 2019.

The fiber-optic cable was attached to a thermistor-piezometer cable weighted with a ~2-kg chain 280 mm above the piezometer diaphragm; this assembly was lowered to 1043 m below the ice surface, whereupon pressure readings stopped increasing, and then raised 50 mm above the inferred bed location. We used a Solifos BRUsens fiber-optic sensing cable containing two single-mode fibers [OS2, for separate distributed acoustic sensing covered by Booth et al. (40)] alongside four multimode (OM3) bend-insensitive fibers for DTS in a duplex arrangement. A stainless steel basal turnaround with an outside diameter of 15 mm and internal fiber loop of 5 mm diameter was used. The cable diameter was 4.8 mm, with fibers enclosed in a gel-filled stainless steel capillary tube further protected by stainless steel wires and a polyamide outer sheath. The hydrostatic pressure resistance of the cable was 30 MPa, with a specified installed tensile stress resistance of 3 kN. Our own tensile tests show that cable failure at 21°C actually occurs in excess of 5 kN (fig. S5).

DTS system and calibration

Silixa XT-DTS with a 0.25-m sampling resolution and ~0.65-m manufacturer-stated spatial resolution was configured to take measurements from the two multimode (OM3) fiber loops in sequence, connected by a basal turnaround assemblage. The manufacturer-stated spatial resolution refers to the distance over which 80% of a ~30°C step change in temperature can be observed along the fiber [(61), see (62) for further details]. The DTS unit was powered via a low-voltage disconnect by two 116-Ah deep-cycle lead acid batteries connected to three 45-W solar panels and contained in a large Pelican case. The batteries were topped up by a generator during field operations. To remove the need for fusion splicing in the field, E2000/APC pigtail connectors were preattached using IP68 conduits leading from a small Pelican case containing the fiber-optic splices. The measurement averaging time was originally set at 2 min with near continuous operation from 5 to 21 July whereupon it was set to 10 min with a rest time of 40 min, with the rest time increased to 4 hours on 23 July to reduce power consumption for unattended operation.

We used the open source Python “dtscalibration” package (63) to calculate a temperature time series, T (K), with uncertainty bounds at time t (s) and depth z (m), which builds on the standard DTS equationT(z,t)=ζln(PS(z,t)PaS(z,t))+C(t)+0zΔα(z)dz(1)PS and PaS are the Stokes and anti-Stokes power and ζ is the sensitivity of (anti-)Stokes scattering to temperature which is dependent on fiber material and does not vary in time. C is a dimensionless function of the alignment and sensitivity of the optical system, connectors, and splices that can vary with interrogator temperature and movement. The integral term is the integrated differential attenuation, where ∆α(z′) (m−1) is equal to the attenuation of the Stokes wavelength minus the attenuation of the anti-Stokes wavelength and is nonzero as attenuation is related to wavelength. Strain, sharp bends, fiber microbending, and splices are known to affect differential attenuation (18, 62, 64). To calculate the integrated value of this term at each point, a double-ended calibration must be used, where the DTS shoots laser pulses first in a forward direction and then in a reverse with the reasonable assumption that temperatures are equal in the forward and reverse directions (the basal turnaround assemblage marks the halfway point). The correction from integrated differential attenuation may not be perfect for highly localized, high-magnitude, cable deformation occurring on the scale of the DTS spatial resolution (0.65 m). This may explain the very fine scale changes of <0.1°C between 889 and 935 m that occur after 8 August (Fig. 2D and figs. S1 and S2); these could result from heterogeneous and high-magnitude cable deformation leading to minor over- or undercompensation in the calibration procedure.

We further used the dtscalibration package to align and correct for asymmetric losses (greater attenuation in one direction than the other) in the turnaround splice, thought to be caused by differential modal spreads in (anti-)Stokes wavelength following a splice within 50 mm of a 5-mm bend radius turnaround. The dtscalibration package uses a Monte Carlo simulation to approximate the probability density function and confidence intervals of the estimated temperature and its variation with depth and time. A full description of the dtscalibration package can be found within des Tombe et al. (63) and accompanying GitHub documentation. The spatial temperature change within the absolute 95% confidence intervals produced with this method (Fig. 2, A to E) has an instrument stated precision of <0.01°C under our calibration procedure (61); that is to say, we have very high confidence in the observed temperature changes with depth with the possible exclusion of the <0.1°C changes discussed above but that the entire profile can shift within the absolute confidence intervals. The DTS calculates distance using the refractive index of the cable and gives a total depth of 1062 m; all depth values were therefore scaled by a factor of 0.9821.

If ζ is assumed constant as in van de Giesen et al. (64), then only one reference section is required for calibration purposes. We use a ζ value of 476.53 K, obtained in a separate laboratory-based calibration procedure using the same DTS unit and cabling as used in the field and two precalibrated (class A, ± 0.15 + 0.002×T °C) PT100s placed into well-mixed polystyrene-insulated water baths with 25-m fiber sections at 0°C and ambient room temperature (18°C, left to settle overnight). Two hours of measurements were taken, but only the 20 min showing the lowest Stokes wavelength variance were used in the calibration. We additionally maintained an ice bath during field operations measured by the same thermistors described above. However, insufficient water mixing and a high resultant temperature variation (±0.4°C) led us to use a section of the temperate zone between 100 and 1010 m in place of the ice bath for field calibration of C and integrated differential attenuation, as the temperate zone exhibited far lower temperature variability (fig. S9). The ice bath was ultimately dismantled for unattended operation due to concern about high winds. The temperate zone temperature was calculated asTm=T0γ(PiP0)(2)where Tm (K) is the depth-dependent melting point of ice, T0 (K) is the triple point temperature of ice, γ is the Clausius-Clapeyron constant (K Pa−1), and P0 (Pa) is the triple point pressure of ice. The local (hydrostatic) pressure, Pi (Pa), for a parallel-sided slab of ice is calculated asPi=ρgHcos(θ)(3)where ρ is the ice density taken as 910 kg m−3, g is the gravitational acceleration (m s−2), H is the height of overlying ice (m), and θ is the ice surface and bed slope (taken as 0.96° using a linear regression along flow line using ArcticDEM surface elevation data over 10 ice thicknesses). A Clausius-Clapeyron constant of 9.14 × 10−8 K Pa−1 was used to provide the best fit possible for the thermistors over the last 4 days of the record and assumed constant for the full measurement period. This method has the limitation of needing to assume a constant temperature within the temperate zone but is used here for simplicity and in the absence of a more reliable method of calculating the exact temperature. A separate cable with digital temperature sensors was installed in the borehole; however, this deployment failed. Our recommendation for future studies is to use the cold central section and a temperate zone, if present, for calibration. The section between 1000 and 1010 m, rather than the basal section next to the thermistors, is used as (i) different heat capacities and thermal conductivities mean that the two sensing systems equilibrate at different rates, which has a large impact especially for the first few days, and (ii) no fit could be obtained incorporating all thermistors. As manufacturer-stated thermistor uncertainty is very low, this probably arises from a small thermistor calibration error. (iii) Initial calibration (not incorporated into the final record) using the ice bath as a reference section suggests that this section of the temperate zone exhibits the least temporal temperature variation, and (iv) this approach allows more noticeable temperature changes at the top and base of the temperate zone to be measured. A comparison of thermistor and DTS recorded temperature in temperate ice can be found in the Supplementary Materials (fig. S9).

Equilibrium ice temperature

The undisturbed ice temperature (Teq in K) and root mean square error were estimated following (34, 65) asTeq=T(t)(Q4πk(ts))(4)where Q (W m−2) is the heat energy released per unit length of the borehole during drilling, k= 2.1 W m−1 K−1 is the thermal conductivity of ice, s is the delay in seconds until the onset of asymptotic cooling, and T(t) is the temperature in Kelvin at time t in seconds. Following Ryser et al. (25), the parameters Q, s, and T0 were determined by fitting Eq. 4 to temperature at the given depth during asymptotic cooling using the SciPy optimize.curve_fit nonlinear least squares fit function. This method produced errors at lower temperature change rates where the final observed temperature is closer to 0°C and was not used beyond 930 m depth. Borehole diameter modeling following Greenler et al. (60) for a target diameter of 12 cm ensured that latent heat energy input was as low as possible. The final observed temperature (13 August, 6 weeks following installation) is within 0.2°C of the predicted equilibrium ice temperature at all points above 930 m (Fig. 2, A and B). A further 40 days would be required for the fitted exponential cooling curve to fall within 0.05°C of the predicted equilibrium temperature.

Strain from differential attenuation

Integrated differential attenuation was also used separately to the dtscalibration package as a measure of strain. This was calculated following (64) aszz+ΔzΔα(z)dz=ln(PS(z+Δz)PaS(z+Δz))ln(PS(z)PaS(z))+ln(PS(z)PaS(z))ln(PS(z+Δz)PaS(z+Δz))2(5)where ⇒ and ⇐ indicate forward and reverse measurement directions, respectively, and ∆z is the sampling resolution. As Raman scattering has not previously been used as a measure of strain, we conducted laboratory testing using a Testometric X350-10 stress and strain gauge in a temperature-controlled (21 ± 0.5°C) room. A 2-m section of cable of the same specification as used in the field was strained by 1.85% (ductile deformation began at 0.65%) at a maximum force of 5.8 kN. A 10-min measurement averaging time was used at each strain increment. This produced a similar disturbance to that seen in Fig. 2 (fig. S4).

Anomaly-208 heat transfer and deformation

We used a simple 1D heat diffusion model adapted from (66) to investigate temperature evolution in the region of Anomaly-208 with a vertical resolution of 0.5 m under (i) no heat inputs or sinks and (ii) heat inputs and sinks prescribed manually as square waves to maintain the anomaly at a steady state. The energy input and sink required to sustain the anomaly can be used to estimate the deformation rate. Strain heating, Qs, is expressed as Qs=tr(τϵ·), where τ is the deviatoric stress tensor and ϵ· is the strain rate tensor (10, 16). Stress and strain are related via ϵ·=Aτen1τ (44, 45), where τe is the effective stress and τe2=12tr(τ2), and A is the creep parameter, which we take as 7.7 × 10−25 s−1 Pa−3 following aggregate values in Cuffey and Paterson [(16), table 3.3] with the range 6.7 × 10−25 to 8.7 × 10−25 s−1 Pa−3, used to test parameter sensitivity and the flow exponent n ≈ 3. Following Seguinot et al. [(10), equation 6], effective strain rate (dimensionless), ϵ·e2=12tr(ϵ·2), is related directly to energy dissipation asϵ·e=(12QsA1n)11+1n(6)

We follow the assumption of Seguinot et al. (10) that deformation occurs within a 2D cross section with vertical, z, and horizontal along-flow, x, dimensions (with y perpendicular to flow), yielding effective strain in Cartesian components, ϵ·e2=12(ϵ·xx2+ϵ·zz2)+ϵ·zx2, and giving stress units of N m−2 and strain heating units of W m−2. Assuming incompressibility, ·ϵ·=ϵ·xx+ϵ·zz=0, gives effective strain as a function of longitudinal and shear components, ϵ·e2=ϵ·xx2+ϵ·zx2. We neglect the longitudinal component as we are not in a zone of rapid ice thinning, but note that the effect of a longitudinal component would be to reduce the simple shear component and relate effective strain to deformation velocity, u (m s−1), as (67) (equation 2.6a)ϵ·e12uz(7)

Finally, assuming a constant deformation velocity with depth and combining Eqs. 6 and 7 givesu2z(12uz)(8)where ∆z= 5.5 m is the height over which deformation is occurring. This gives a deformation velocity of 14.2 m a−1 with limits 13.7 to 14.7 m a−1 for creep parameters 6.7 × 10−25 and 8.7 s−1 Pa−3, respectively.

We further constrain deformation velocity by assuming that the heat sink is a result of an advective flux prior to cable installation with constant velocity through depth within a z x plane to give [e.g., (68)]u=QaρCT(9)where Qa is our advective energy sink (W m−2). Here, the depth over which advection occurs does not need to be considered. We take the specific heat capacity of ice at −17°C as 1.972 kJ kg−1 K−1. Taking a value of ∆T of 0.2 K, the difference between the peak and trough of the anomaly gives a velocity of 8.2 m a−1.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We thank A. Andreasen and Uummannaq Polar Institute for hospitality. We are grateful to S. Tyler and CTEMPs Instrument Facility for technical guidance; L. Greenler for providing code for modeling borehole diameter; M. MacKie, S. Peters, M. Prior-Jones, and E. Dawson for assistance in the field; and N. de Battista for help with strain gauge testing. We thank A. Aschwanden and two other anonymous reviewers for comments. Funding: This research was funded by the European Research Council as part of the RESPONDER project under the European Union’s Horizon 2020 research and innovation program (grant 683043). R.L. and T.R.C. were supported by Natural Environment Research Council Doctoral Training Partnership studentships (grant NE/L002507/1). B.H. was supported by a HEFCW/Aberystwyth University Capital Equipment Grant. Author contributions: R.L. and P.C. designed the DTS experiment. P.C. led the field operations. B.H. and S.H.D. led the drilling work. R.L., P.C., B.H., S.H.D., T.R.C., C.M.S., and A.B. collected data in the field. R.L. analyzed the DTS data with support from B.d.T. and B.S. C.K. advised in laboratory experiments for integrated differential attenuation. R.L. wrote the manuscript with help from all coauthors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper, the Supplementary Materials, and the Apollo repository (available at

Stay Connected to Science Advances

Navigate This Article