The cooling signature of basal crevasses in a hard-bedded region of the Greenland Ice Sheet

. Temperature sensors installed in a grid of 9 full-depth boreholes drilled in the southwestern ablation zone of the Greenland Ice Sheet consistently record cooling over time within the lowest third of the ice column. Rates of temperature change outpace cooling expected from vertical conduction alone. Additionally, observed static temperature profiles deviate significantly from modeled purely diffusional thermal profiles, implying strong non- 10 conductive heat transfer processes within the lowest portion of the ice column. Although numerous heat sources exist to add energy and warm ice as it moves from the central divide towards the margin such as strain heat from internal deformation, latent heat from refreezing meltwater, and the conduction of geothermal heat across the ice-bedrock interface, identifying heat sinks proves more difficult. After eliminating possible mechanisms that could cause cooling, we find that the observed cooling is a manifestation of previous warming in basal ice. Thermal decay 15 after latent heat is released from freezing water in basal crevasses is the most likely mechanism resulting in the temporal evolution of temperature and the vertical thermal structure observed at our site. Basal crevasses are a viable englacial heat source in the basal ice of Greenland’s ablation zone and may have a controlling influence on the temperature structure of the near basal ice. continuous record of englacial temperatures identify of heat flow at our field location. of 15CA and 14SA). Additionally, borehole pairs 14SA, 14SB and 15CA, 15CB are located within 10 m of each other and show similar rates and extent of cooling over time, establishing repeatability and viability of our measurements and methods. our field site in 3-dimensions, would result in cooling that we have not captured in our vertical 1-D modeling. It is plausible that boreholes intersect these regions of ice that have been thermally altered by basal crevassing.


Introduction
Rates of ice deformation and the onset of basal sliding are regulated by the thermal regime of ice sheets (Hooke, 1981;Weertman, 1957), which makes ice temperature an important physical variable to be measured in the field and modeled to determine current and future ice motion. As ice moves from the central divide towards the margins, heat is conducted into the ice column from the surface and basal boundaries and is introduced through sources such as 25 energy released by internal deformation and refreezing of meltwater (Cuffey and Paterson, 2010).
Mechanisms of heat transfer and heat sources in Greenland's ablation zone have been investigated by comparing sparsely collected temperature profiles to output from thermomechanical models (e.g. Iken et al., 1993;Funk et al., 1994;Lüthi et al., 2002;Harrington et al., 2015;Lüthi et al., 2015;Poinar et al., 2017). When attempting to replicate the vertical temperature structure, models consistently produced ice temperatures that are lower than 30 collected temperatures. Previous modeling and observational studies highlight the appreciable influence of refreezing meltwater as a source of latent heat in Greenland's ablation zone (e.g. Phillips et al., 2013;Harrington et al., 2015;Lüthi et al., 2015;Poinar et al., 2017;Seguinot et al., 2020). Outside of Greenland, observations of englacial warming from latent heat sources have been demonstrated as early as the 1970s (Jarvis and Clarke, 1974).
Here we investigate heat transfer processes using an extensive dataset collected from over 300 temperature 35 sensors emplaced in a grid of 9 full-depth boreholes in southwestern Greenland. Our data provided sub-daily continuous records of ice temperature over a 2-3 year period, longer than many similar temperature time-series datasets. The high temporal and spatial resolution allowed us to examine how the vertical thermal structure evolves in both time and space during the study period. Our continuous record of englacial temperatures enabled us to identify dominant drivers of heat flow at our field location. 40 2 Methods

Field site
Nine boreholes were drilled in the southwestern portion of the Greenland Ice Sheet in a sector that discharges in the land-terminating outlet glacier, Issunguata Sermia. The field site was located approximately 33 km from the margin and had a surface elevation of 980 m above sea level with ice depth ranging from 640 m to 680 m. The dimensions 45 of the borehole grid were approximately 250 m N-S and 500 m E-W. Figure 1 shows the layout of the instrumented field site.
Borehole experiments conducted immediately after drilling indicated the ice rests on bedrock with perhaps a veneer of sand and gravel approximately a decimeter thick . The subglacial drainage system was indicative of one overlying bedrock and governed by hard-bedded physics. Mean annual air temperature 50 measured at the site was -10.5°C and ablation rates were 2-3 m yr -1 (Hills et al, 2018). Ice moved to the west (266°) at ~100 m yr -1 and inclinometry-derived deformation rates showed that basal sliding comprised ~96% of ice motion, with internal deformation contributing the remaining 4% (Maier et al., 2019). Airborne radar measurements indicated that the bed elevation increases slightly from east to west across the site such that the ice flowed up a minor reverse bed slope (Paden et al., 2014).

Data collection
A hot water drill was used to melt boreholes ~10 cm in diameter through the full depth of the ice column. Cat-5 cables instrumented with TMP102 temperature sensors with 12 bit analog-to-digital converters were installed in the boreholes after drilling and logged with a custom data logger at the surface. Each temperature sensor was controlled by a separate downhole microprocessor that ensured the sensor was powered only as minimally required to obtain a 60 reading. The paired microprocessor handled all digital communications, power control, and data processing to minimize errors caused by self-heating. Sensors were spaced every 10 m in the lowest 150 m of the ice column and every 20 m for the remainder of the depth. The uncalibrated accuracy of the sensors was 0.5 °C and the resolution was 0.0625 °C. We performed a freezing-point calibration prior to complete borehole freeze-up using a Clausius -Clapeyron gradient of -7 x 10 -8 °C Pa -1 (Cuffey and Paterson, 2010). However, due to uncertainty in the Clausius -65 Clapeyron slope given dissolved air and impurities within the ice and the ambient pressure at the time of calibration, we are confident in the sensors' accuracy to approximately 0.5 °C. Data loggers at the ice surface recorded temperature readings every 4 hours between the date the borehole was drilled in either July 2014 or July 2015 and the end of the data collection period in July 2017. Each temperature measurement that the sensor records was an average of 64 temperature readings, which was then registered at the nearest resolution step of the sensor.

70
Our data allowed us to examine both the vertical temperature field within the ice column as well as the temporal changes in temperature. The temperature of a hot water drilled borehole decays towards the ambient ice temperature over time. After 1 year < 0.1% difference remains between the borehole temperature and ambient ice temperatures. For ambient ice temperatures below -1 °C, this is well below sensor resolution and does not affect our observational data. Additionally, regressing equation 24 in Humphrey and Echelmeyer (1990) between 1 and 3 years 75 after borehole closure results in a rate of residual cooling on the order of 10 -11 °C yr -1 . To avoid any potential contamination of both the static and temporal trend of the sensor temperature data from the thermal disturbance of installation, we disregarded temperatures recorded in the first year after drilling in our analysis. To display a snapshot of the ambient vertical temperature field with minimal effects of the thermal disturbance from hot water drilling, we report temperatures collected one year after boreholes are drilled (Humphrey and Echlemeyer, 1990; 80 (Supplementary Materials).
We also neglected the first year of measured temperatures in the temporal records of temperature change to eliminate the drilling disturbance. After this time, the measurements of temperature in the environment isolated from electrical noise resulted in a time series of temperature changes that is mostly noise-free and at the absolute limit of the sensor's resolution. The data was cleaned of the few readings recording more than ± 2 resolution steps from the 85 preceding and succeeding measurements, since extreme jumps were likely caused by digital transmission errors between the sensor and data logger at the surface. Although englacial temperatures likely vary smoothly over time, the digital representation from sensor data shows temperatures changing over time in a stair-stepping pattern. Since the temporal temperature variations are small in relation to the sensor resolution, the dataset is highly characterized by this stepped behavior (Fig. S1).

90
The rate of temperature change recorded by a sensor was determined by a least squares regression of the data recorded over time. Often the first and last digital steps of this record had fewer entries recorded due to the restrictions of the sampling window. Linearly regressing our stepped data to determine the rate of change of temperature required careful treatment to prevent a biased estimate from truncated observation periods. Assuming that temperature change is close to linear with time, we expected that with a longer record, the sensor would record 95 approximately the same number of readings at each resolution step. Short data lengths occurring at either the beginning or end of the records were likely truncated by the restricted time period of our study. To correct for our restricted sampling time, we padded the first and last steps in our data to match the length of a fully-recorded resolution step. This was performed by copying the required length of data from a full resolution step, and pre-or post-appending it to the truncated data, while equating the temperatures. This procedure is illustrated in Fig. S2, 100 which also shows the difference in derived rates of temperature change over time (Supplementary Materials).

Results
The vertical temperature field collected from the nine boreholes one year after temperature sensors were installed is shown in Fig. 2, with two boreholes only being instrumented with temperature sensors in only the lowest half of the ice column. Temperatures are lowest in approximately the middle of the ice column, ranging from ~ -12 °C to -10 105 °C and increase to ~ -11 °C to -4.5 °C at the surface. Temperatures also increase towards the bed, reaching temperatures that range from approximately -0.9 °C to -0.5 °C. To establish a depth reference for sensor comparison, and since we are most interested in near-basal ice temperatures, we have plotted temperature profiles with respect to height above the bed. signature is seen in the lowest third of the ice column in most boreholes, with some boreholes warming or cooling in locations outside of this region. Most boreholes cool at a rate of a few hundredths of a degree Celsius per year in the lowest portion of the ice column. Only five full-depth profiles are shown, as temperature sensors in boreholes 15N and 14W failed for most of the study period. Boreholes 15CB and 14SB were only instrumented in the lowest portion of the ice column.

115
While the observed temporal temperature changes are small, we believe the trends are a real signal and not an artifact of sensor drift or random noise. Analog-to-digital systems can exhibit both noise and drift over time; however, the expected magnitude of drift or aging is approximately an order of magnitude less than our observed trends. In addition, the coherence and lack of obvious bias in the data gives us confidence. We observe coherent sections of ice that cool (the bottom of 14SA, 14N, and 15CA) while other portions warm (top of 15CA and 14SA).

120
Additionally, borehole pairs 14SA, 14SB and 15CA, 15CB are located within 10 m of each other and show similar rates and extent of cooling over time, establishing repeatability and viability of our measurements and methods.

Discussion
We focus on observed basal cooling, since cooling of ice near the bed is enigmatic. Multiple heat sources are available to warm ice as it flows from the interior towards the margins. Heat sinks, however, are not obvious, other 125 than vertical thermal conduction from the relatively warm ice at the bed into the cold central core. We proceed by investigating the potential causes of cooling in the lower part of the ice column.

Diffusive cooling driven by the observed static temperature field
We compare the magnitude of observed cooling over time to the temporal temperature changes that would be produced from purely vertical thermal diffusion given the shapes of our collected temperature profiles. We use the 130 1-D thermal diffusion equation in the vertical dimension: where T is temperature, and t is time; κ, , and C are respectively, the thermal diffusivity, density, and specific heat capacity of ice. The constants used in this study are reported in Table S1. We set the source or sink term, Q, to zero, assuming only thermal diffusion occurs for the following calculations. The spatial second derivative is derived from our observed temperature profile. The temperature data has small amounts of noise as well as calibration offsets 135 between sensors. Derivatives of real data, with any noise or error, tends to amplify the noise. To reduce noise levels and focus on the large-scale characteristics, we smooth the profiles using a 5 th degree polynomial. We choose a 5 th degree polynomial because it smooths over sharp kinks within the profile from calibration errors while maintaining the original curvature, particularly the negative concavity of the temperature profiles observed near the bed. The full-depth observed rates of temperature change in Fig. 3 are also smoothed using a 5 th degree polynomial to 140 facilitate comparison between the changes in temperature expected from collected profiles.
Comparison of our data and expected temperature changes from the analytical diffusion equation (Eq. 1), indicates the curvatures in our temperature profiles cannot produce the observed temporal changes of temperature over time through diffusive vertical heat flow alone. To illustrate the discrepancies between calculated temperature changes caused by vertical thermal diffusion and the observed temporal temperature change, we calculate the 145 additional sources or sinks of heat required throughout the ice column to sustain the observed temperature change through time. Figure 4 presents a vertical distribution of the required Q from the thermal diffusion calculations to produce the recorded temperature changes. Although there is considerable variability, all boreholes require additional heat sinks within portions of the lowest half of the ice column, demonstrating that cooling near the bed cannot result from vertical thermal conduction alone.

150
In addition, we also highlight regions of the temperature profiles that are not only warming or cooling anomalously, but that may be warmer or cooler than would be expected given the average vertical temperature field at our site. The shape of the site-averaged temperature profile results from ice integrating changes in boundary conditions and all upstream thermomechanical processes. Any deviation from this averaged profile within each collected temperature profile will highlight thermal regions that reflect changes either to the local boundary 155 conditions or heat sources or sinks. These changes must have occurred recently for sensors within the closely spaced boreholes to capture temperature fields that do not match the average vertical temperature field at the site.
Discrepancies between each collected profile and the site-average temperature profile are shown in Fig is colder throughout the full-ice depth than the mean thermal profile, but this borehole also is where sensors record the most warming throughout the profile. Spatially, the southeast portion of our field site does not have a cooling 165 signature in relatively warm ice near the bed, which is evident throughout the other boreholes that have time-series of temperature measurements.

Possible heat sinks
Pervasive cooling in the lowest portion of boreholes where ice is warmer than the site-average profile in some, but not all, of our borehole data is enigmatic, since multiple heat sources exist near the bed, but identifying heat sinks is 170 problematic. In addition, the non-uniform nature of cooling near the bed requires explanation. We discuss several possible heat sinks below.

Melting in basal temperate ice
Ice near the bed at our field site is at the calculated pressure-melting temperature (Fig. 2). Some melting can occur in this temperate ice as the pressure field varies due to ice flow, and transient increases in pressure decrease the melting reaches the pressure melting temperature in the lowest ~5-10 meters at our site, with an even smaller temperate layer upstream. This process cannot explain cooling outside of this small area near the bed.

Basal refreezing from a thinning ice column
The ice sheet thins slightly along flow at our field site; the bed elevation increases while the surface elevation 180 remains relatively constant. With a lower ice overburden pressure, the pressure-melting temperature increases (Cuffey and Paterson, 2010). Any interstitial water present in the ice would then be colder than the new melting temperature and would refreeze. Lüthi et al. (2002) estimated the water content of the temperate layer of ice near the bed of Jakobshavn Isbrae has a moisture content of ~ 1%. However, refreezing releases latent heat into the surrounding ice, resulting in warming rather than cooling. Because the temperate layer is confined adjacent to the 185 bed, this process would also be localized.

Downward advection of the cold central core
All full-depth temperature profiles display a characteristic cold central core (Fig. 2

Long-term changes in boundary conditions
Variations in the boundary temperatures will propagate into the interior of the ice column. Here we address several of the most likely boundary temperature changes to assess how they could affect the temporal evolution of the ice at our site. The The length-scale for changes in surface temperatures during the Little Ice Age 500 years ago is ~130 m, which suggests hardly any thermal signature from the surface boundary exists below ~ 400 m depth.
Equation (2) also provides an estimate of the timescale of diffusion, illustrating changes in ice temperature from the end of the last interglacial period would have diffused away through entire ice column after ~12,000 years.
The most significant thermal change to the basal boundary is the transition from frozen to temperate, which 210 occurs upstream of our site. This change introduces a warmer basal boundary relatively recently in terms of the ice flow; however, this would create a region of colder ice near the bed that is being warmed by conduction from the warm bed, not cooled as we observe.

Changes in ice motion
The interior ice sheet, which is frozen to the bed, moves largely by internal deformation and accumulates dissipated 215 heat near the bed where strain rates are highest. At our site ice moves primarily by basal sliding (Maier et al., 2019) and this strain heat source is negligible, although frictional heating still exists at the base. Downstream from the region where this transition from deformation to sliding occurs, the ice near the bed would lose the strain-derived thermal energy and would cool as the heat diffused upward into the cold core. Although the location of the upstream transition from a frozen to thawed bed conducive to basal sliding is not precisely known, it is on the order of 200 km 220 up-glacier from our field site (MacGregor et al. 2016), which translates to several thousand years for ice flow. Any heat anomalies would have decayed using Eq. (2). Additionally, we have shown that the shape of the temperature profiles that would result from previously accumulated strain heat near the base will not produce the magnitude of observed cooling from vertical thermal conduction alone. An additional heat sink must be required.
Although the frozen to sliding transition is probably the largest change to the strain heating term, analysis

Thermal decay from basal crevasses
With the lack of obvious heat sinks, we investigate the possibility that the observed cooling in the lower depths of the ice occurs as thermal decay after latent heat is released from freezing water. Temporally and spatially uniform 235 heating by refreezing of disseminated water would not lead to our observations, since horizontal thermal diffusion would not occur in laterally extensive sections of isothermal ice. However, widely spatially separated transient heating events, such as dispersed basal crevasses or moulins, only warm local ice, leaving significant regions of colder, unheated ice to facilitate horizontal diffusion, allowing the heated ice to cool back to ambient temperature.
Since heat would flow horizontally, not vertically, it would not be revealed by our vertical thermal diffusion 240 modeling.
A similar pattern of cooling in ice warmer than the 1-D thermal diffusion model produces is found in the upper portion of 14N (Fig. 4). Hills et al. (2017) discussed the thermal anomaly exhibited in borehole 14N; the authors recognized the relatively high temperatures of ice in the upper portion of the borehole when compared to others in the grid (Fig. 2) and noticed that temperature sensors recorded decreasing temperatures over time. Hills et 245 al. (2017) attributed this localized thermal anomaly to a refrozen crevasse that was observed on the ice surface nearby. The lowest 100 m of ice in boreholes 14SA, 14N, and 15CA all exhibit nearly identical thermal characteristics to the uppermost portion of 14N: ice is warmer than the site-average profile and it cools through time.
Basal crevasses can have a thermal influence on discrete but large regions of ice. Because these fractures are distinct features, the fact that we do not see a clear cooling signature in relatively warm ice throughout the base 250 of the ice column at our site is not surprising. Slight warming in cold ice in borehole 15E and no clear signal in borehole 15S provides evidence that the ice is being affected by a thermal process that is not distributed in nature.
Multiple fractures can open similar to crevasses at the surface, and their growth is likely vertically restricted above the bed due to the cold, stiff central core of the ice column which reduces fracture propagation.
Therefore, these features only affect the lower portion of the ice column. Additionally, basal crevasses likely are not 255 simple planar features; strong deformational gradients near the bed can alter their shape and cause them to intersect more recently formed crevasses. Although basal crevasses may be more complex than an en échelon arrangement of vertical fractures, for simple modeling we focus on the thermal disturbance by basal crevasses that allow the influx of a plane of water from the bed up into the cold ice. The refreezing of this injected water releases heat into the surrounding ice. Accurately modeling the thermal disturbance created by a basal crevasse is hampered by our lack of 260 knowledge of these features, since they have only been directly observed under mountain glaciers (Harper et al., 2010). Even one-dimensional thermal modeling in ice surrounding a crevasse requires knowledge of the width of the crevasse and the presence of other nearby crevasses, since closely spaced water filled crevasses greatly reduce the refreezing rate (Jarvis and Clarke, 1974). In addition, even the refreezing problem is complicated by the Stefan boundary condition as the thickness of the crevasse is reduced by refreezing (Carslaw and Jaeger, 1959).

265
We do not attempt detailed modeling of the thermal disturbance of basal crevasses at our field site, but we can approximate the thermal disturbance created by a single or by multiple basal crevasses to compare with our data.
We are interested in the thermal behavior of ice in the region of a basal crevasse at times much longer than the crevasse remained open, so we model the crevasse as an instantaneous planar heat source, considering only the total latent heat released during the refreezing. We examine the thermal disturbance at a long timescale in the half-space 270 adjacent to the crevasse using an analytic solution presented by Carslaw and Jaeger (1959): Here T is temperature, l is the linear distance away from the heat source, t is the time that has passed since the heat was released, Q is the amount of heat released per unit area of the plane, ρ, C, and κ is the respectively the density, specific heat capacity, and thermal diffusivity of ice.
We assume a crevasse thickness of 1 m, matching the upper limit of previously observed basal crevasse 275 widths (Harper et al., 2010). The release of heat initiates a wave of thermal energy propagating away from the source (Fig. 6A). The diffusion of the warming front causes ice temperature to increase before cooling asymptotically back to the ambient ice temperature (Fig. 6B). After the crevasse refreezes, the remaining period of warming is short compared to the timescale of cooling, thus the likelihood of temperature sensors installed in a previously crevassed area recording warming is low. To compare this theoretical result with our data, we average the 280 rates of cooling observed in the lowest 100m of the boreholes and show the range of these averages in Fig. 6C. The theoretical curves depend strongly on the original size of the crevasse, and on the distance away from the crevasse; however, we have chosen values that show that this process can produce rates of cooling of the same magnitude observed in the borehole array.
It is tempting to estimate the time since water-filled basal crevasses were formed using our extensive 285 temperature data. However, such an attempt is under-constrained since similar thermal signatures could be produced by many crevasse geometries and arrangements in time and space. Instead we examine the same crevassing scenario above as a baseline to investigate the amount of time elapsed since refreezing to produce observed cooling rates given various distances between boreholes and hypothetical crevasses. Equation (3)  assuming an ice velocity of 100 m yr -1 , the ice currently at our field site would have been 1-3 km up-glacier, which is a location of current surface crevasses and a significant change in slope.

295
Basal water pressures at or above ice overburden pressure promote the growth of basal crevasses 10s of meters above the bed (van der Veen, 1998). Longitudinal extension further facilitates crevassing at distances over 100 m above the bed (van der Veen, 1998). Subglacial water pressures in Greenland's ablation zone have consistently been measured at or near ice overburden pressure (e.g. Iken et al., 1993;Meierbachtol et al., 2013;Andrews et al., 2014;, van de Wal et al., 2015Doyle et al., 2018), and have peaked at ~ 1.15 300 times the ice overburden pressure at our field location (Wright et al., 2016). Longitudinal extension, while not temporally or spatially prevalent, may result from ice flow over uneven bed topography or spatial changes in basal friction. Additionally, observed stick-slip behavior and hydraulic jacking from rapid delivery of surface meltwater through passive seismic surveys could provide another mechanism for basal crevasse formation (Moore et al., 2013).
It is unlikely that a single basal crevasse is responsible for the cooling that we observe at the field site.

305
More plausibly, the ice at our field site was fractured by multiple crevasses up-glacier. We cannot precisely determine the number or spacing of crevasses; however, a 1-D modeling analysis shows that three 1 m thick refrozen crevasses would lead to dispersed cooling over a large area of the bed after 20-30 years (Supplementary Materials). The 1-D temperature field evolution illustrates the warming that occurs due to the release of latent heat immediately as the crevasse completely refreezes, followed by dispersed cooling surrounding the crevasses over 310 time (Fig. 7). Given negative heat fluxes are required in the lowest portion of all boreholes (Fig. 4) to create observed cooling, thermal energy likely diffuses away from our field site in 3-dimensions, which would result in cooling that we have not captured in our vertical 1-D diffusion calculations. It is plausible that boreholes intersect these regions of ice that have been thermally altered by basal crevassing.
Moulins are unlikely to create the observed thermal observations. While the internal architecture of moulins 315 is likely complex, as a simple model, we conceptualize moulins as linear pipes that release thermal energy from the surface to the bed because of how these features are manifested at the surface. Cooling in relatively warm ice is not observed throughout the entire ice column in each borehole, which we would expect to see from a linear heat source extending over the full ice depth. Moulins have been shown to migrate horizontally with depth (Holmund, 1988), which could possibly avoid this objection; however, numerous moulins would be needed to generate the cooling we 320 observe in multiple boreholes. In the ablation zone surrounding Issunguata Sermia, active moulins have a mean density of 0.1 per square kilometer (Smith et al., 2015), and Catania and Neumann (2010) found relict moulins to be spaced ~ 1 km apart, which suggests that it is unlikely that multiple moulins have thermally altered the instrumented block of ice with a surface area of ~ 0.125 km 2 .

325
The existence of basal crevasses has not been widely considered outside of alpine settings; however, these features likely impact the bulk water storage, basal water flow, and ice rheology in Greenland's ablation zone. Harrington et al. (2015) examined full-depth ice temperature data collected in boreholes at a site down-glacier, but on the same flowline as this study, and found temperate layer growth could not be reproduced by thermal modeling incorporating strain heat and longitudinal compression. The authors suggest basal crevasses allow latent heat to be released well 330 above the bed and create the thick temperate layer.
The lack of a thick basal temperate layer indicates the transient thermal influence of basal crevasses. Ice viscosity is highly temperature dependent, thus basal crevasses could create regions of enhanced deformation, although previous observations in Greenland have shown that deformation rates in temperate ice are often lower than theory predicts (e.g. Doyle et al., 2018;Maier et al., 2019). In addition, water content in 335 basal temperate ice also influences rates of glacial ice deformation. A previous study at the down-glacier site, found a higher water content in the basal temperate ice layer revealed by ground-based radar, than is used in modeling studies (Brown et al., 2017). While unable to detect macro-spaces for bulk water storage in the radargrams, the authors speculated that basal crevasses provide macro-spaces for water storage based on the thermal evidence from Harrington et al. (2015), which results in the high water content observed.

340
Although basal crevasses have only been directly observed in the alpine setting, their existence under the Greenland Ice Sheet has strong implications for the basal hydrology. Basal crevasses have been observed creating englacial water flow above the bed (Harper et al., 2010) and even allowing complex three-dimensional crossing water flow paths (Huzurbazar and Humphrey, 2008). These features have the potential to modulate subglacial discharge and efficiently transfer water through the subglacial drainage system.

6 Conclusion
Temperature measurements that began one year after sensor installation in 9 boreholes and span 2-3 years show that ice cools in the lowest portion of the ice column, with all boreholes displaying cooling within the lowest third of the ice column. The main thermal sink is conduction of heat from this region up into the cold central core of the ice.
Vertical diffusion calculations shows this cooling to be insufficient to explain our data. In addition, the lowest 350 section of ice in some boreholes is warmer than the mean thermal profile at our field site, which suggests relatively recent and localized changes to the boundary conditions or near-basal heat sources. The observed cooling regions of the boreholes are often the same regions of ice that are warmer than the site-average profile.
After eliminating various potential heat sinks, we conclude the observed thermal anomalies are created by the cooling of basal ice that has been previously warmed by basal crevassing occurring upstream of our observation 355 site. Latent heat released by refrozen basal crevasses locally warm ice near the bed, but the effects are distributed as thermal energy propagates outwards resulting in dispersed cooling at rates of the same magnitude as our observations. Our results provide additional thermal evidence for the existence of basal crevasses in Greenland's ablation zone, which has strong relevance for the temperature dependent rheology of the basal ice and intriguing implications for the complexity of the basal hydrologic system. Author contributions. NH, JH, and TM designed the field study and installed temperature sensors. IM performed the analysis and wrote the manuscript with contributions from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.