Greenland and Canadian Arctic ice temperature profiles database

Here, we present a compilation of 95 ice temperature profiles from 85 boreholes from the Greenland ice sheet and peripheral ice caps, as well as local ice caps in the Canadian Arctic. Profiles from only 31 boreholes (36 %) were previously available in open-access data repositories. The remaining 54 borehole profiles (64 %) are being made digitally available here for the first time. These newly available profiles, which are associated with pre-2010 boreholes, have been submitted by community members or digitized from published graphics and/or data tables. All 95 profiles are now made available in both absolute (meters) and normalized (0 to 1 ice thickness) depth scales and are accompanied by extensive metadata. These metadata include a transparent description of data provenance. The ice temperature profiles span 70 years, with the earliest profile being from 1950 at Camp VI, West Greenland. To highlight the value of this database in evaluating ice flow simulations, we compare the ice temperature profiles from the Greenland ice sheet with an ice flow simulation by the Parallel Ice Sheet Model (PISM). We find a cold bias in modeled near-surface ice temperatures within the ablation area, a warm bias in modeled basal ice temperatures at inland cold-bedded sites, and an apparent underestimation of deformational heating in high-strain settings. These biases provide process level insight on simulated ice temperatures.

1 Introduction

With the tremendous social implications of sea level change, the past decade has seen a proliferation of simulations of the current and future geometry and dynamics of the Greenland ice sheet (Bindschadler et al., 2013; Goelzer et al., 2020; Aschwanden et al., 2021). These present-day complex ice flow models build upon the legacy of simpler past thermomechanical models (Letréguilly et al., 1991; Huybrechts et al., 1991; Funk et al., 1994; Calov and Hutter, 1996). Due to the high sensitivity of ice viscosity to ice temperature, the thermal state of the ice sheet is a critical element of these simulations (Colgan et al., 2015). At present, however, the englacial temperature fields of even cutting-edge ice sheet simulations remain largely unevaluated against observed temperatures (Aschwanden et al., 2019). While recent studies show potential for deriving internal ice temperatures from satellite or airborne data (Macelloni et al., 2019; Jezek et al., 2022), these techniques are not yet widely employed. There are consequently diverse opinions on Greenland's basal thermal state across the current generation of thermomechanical ice flow models (MacGregor et al., 2016, 2022).

Several different methods for measuring ice temperatures have been used on the Greenland ice sheet and in the Canadian Arctic. The methods include the following: borehole logging, where a temperature sensor is moved up or down the borehole, measuring either “continuously” as the probe moves or stopping to measure at every depth (known as “stop-and-go”; Johnsen et al., 1995; Clow, 2008); sensor strings (this term includes both digital and analog), where thermistors are frozen into the ice, and ice temperatures are recorded at various depths at the same time (Iken et al., 1993; Ryser et al., 2014); and fiber optic distributed temperature sensing (DTS), where the ice temperature is measured near-continuously along the full cable length (Law et al., 2021). We summarize the methods and their advantages and disadvantages in Table 1.

Ice temperature profiles collected in Greenland and the Canadian Arctic have not been systematically compiled into a coherent database. In particular, many pre-1990 ice temperature profiles languished in non-digitized reports or gray literature. This presented a clear motivation to assemble ice temperature measurements into a consistent and comprehensive community resource. Here, we describe our compilation of ice temperature profiles from Greenland and the Canadian Arctic into an open-access database, with well-documented and uniform metadata for each entry. We include Canadian Arctic ice caps in our predominantly Greenland database, as these regions reside within the domain of some Greenland ice flow models (e.g., Tarasov and Peltier, 2004; Gowan et al., 2021).

The earliest temperature profile in our database is from 1950, when a 125 m profile was measured at Camp VI in West Greenland (Heuberger, 1954). Earlier temperature profiles may exist, and incorporating these profiles into the database is an ongoing process. We restrict our database to ice temperature profiles extending well below the depth of the seasonal temperature cycle. At cold, dry sites, this is often approximated to be 10 to 15 m below the surface (Cuffey and Paterson, 2010). Similar to a recent effort to compile surface mass balance observations into a readily accessible common framework (Machguth et al., 2016), we aim to create a community resource that facilitates comparisons between simulated and observed ice temperatures. Here, we describe the format of the database and the sources of our ice temperature data and metadata. We also discuss best practices for comparing observed ice temperatures with simulated ice temperatures from a 3-D thermomechanical model. Finally, we make an appeal to community colleagues to continually update this dataset and provide instructions for them to do so.

Table 1 Advantages and disadvantages of the main methods for measuring deep-ice temperatures. This compilation ignores the absolute accuracy of the ice temperature sensor, which can vary greatly within any one method type.

2 Temperature data

The temperature database currently comprises a total of 95 temperature profiles. In two instances, the same borehole was logged twice at different times (Thomsen et al., 1991). Other boreholes which have been logged more than once also include GISP2-D, GISP2-G, NGRIP-2, NGRIP-S2, and NEEM-D; however, these are currently not included in the database. The database also contains partial temperature profiles, which are considered to make up one profile. This is the case for the two sites of FOXX6 and GULL1 (Ryser et al., 2014). Furthermore, five partial profiles at one site (Lüthi et al., 2002), all situated less than 25 m from one another, were also regarded as one borehole. This means that the 95 temperature profiles in the database are considered to originate from 85 unique boreholes.

See Fig. 1 for an overview of the drill site locations and the database for the complete list of borehole IDs. Borehole IDs denote the general location of the borehole, project-specific borehole IDs if multiple boreholes exist at that location, and the year of the measurement. Borehole temperature records were entered into the database in two different formats, namely (a) tabulated original measurements, provided either by community members or obtained from historical reports, or (b) published graphics of measured temperature profiles digitized using the publicly available digitization tool WebPlotDigitizer (https://automeris.io/WebPlotDigitizer, last access: 3 March 2022; Rohatgi, 2021).

Figure 1 (a) Overview of the drill site locations of temperature profiles contained in the database (plot created using QGreenland software; Moon et al., 2021). Note that for locations where several drill sites exist close to each other, only the common name is shown, e.g., Agassiz. The green and purple boxes indicate regions where numerous profiles have been drilled. (b) Enlarged view of the region near Jakobshavn Glacier. (c) Enlarged view of the region near Isunnguata Sermia glacier. Panels (b) and (c) show positions plotted on top of Sentinel-2 multispectral satellite imagery, with 10 m resolution, from 2019 (MacGregor et al., 2020). All maps are displayed in EPSG:3413 polar stereographic projection, with the easting and northing coordinates given in meters.

Overall, 71 of the 95 profiles were available as tabulated temperature measurements (Heuberger, 1954; Classen, 1977; Stauffer and Oeschger, 1979; Clarke et al., 1987; Thomsen et al., 1991; Iken et al., 1993; Hansson, 1994; Cuffey et al., 1995; Thomsen et al., 1996; Cuffey and Clow, 1997; Dahl-Jensen et al., 1998; Fischer et al., 1998; Lüthi et al., 2002; Kinnard et al., 2006; Buchardt and Dahl-Jensen, 2007; Kinnard et al., 2008; Lemark and Dahl-Jensen, 2010; Rasmussen et al., 2013; Ryser et al., 2014; Harrington et al., 2015; Hills et al., 2017; Zekollari et al., 2017; Doyle et al., 2018b; Seguinot et al., 2020; Hubbard et al., 2021a; Harper and Meierbachtol, 2021; Law et al., 2021). The remaining 24 profiles were digitized from figures (Hansen and Landauer, 1958; Davis, 1967; Paterson, 1968; Classen, 1977; Paterson et al., 1977; Colbeck and Gow, 1979; Gundestrup and Hansen, 1984; Blatter and Kappenberger, 1988; Gundestrup et al., 1993).

Here, we describe both the initial Dataverse database release (Mankoff et al., 2022) and the GitHub living repository (Mankoff, 2022). The Dataverse database (Mankoff et al., 2022) contains four curated files, namely two comma-separated value (CSV) files with temperature and depth-normalized temperature at each site, a Keyhole Markup Language (KML) file of borehole locations, and a metadata file. The GitHub living repository (Mankoff, 2022) is where raw data are collected, curated, and documented in order to create the final data. The repository is termed “living” because it is also the entry point for new data in future database versions.

The GitHub living repository (Mankoff, 2022) consists of high-level, post-processed CSV files ready for use and low-level folders with original sources and notes for each temperature profile. The CSV files include two data files combining all profiles, where one file has depth in units of meters with a step size of 1 m, and another file has normalized depth as 0–1 with a step size of 0.01. This standardization of the profiles from coarse discrete measurements to uniformly finely spaced values was done by interpolating between points, using cubic spline interpolation with no overshoot. A third CSV file combines all metadata for each profile, and finally, a geographic information system (GIS) file (format KML) provides borehole locations. Folders also contain a file named “data.csv” with the original data at the original depth resolution, a README.org file with any relevant notes, including “cleaned” email correspondence when personal communications were involved in recovering the data, or other files detailing the data acquisition.

Where ice temperature profiles are derived from a graphic (Barnes_(B4,D4)_1974, Barnes_T0(975_1976, 91_1977, 81_1978, 61_1977, and 20_1979), Camp Century, Devon72, Devon73, DYE-3, Isua_10-14, Laika_75a-75e, Meighen67, SiteII, and Tuto_D-11), the graphic is included in the folder. These folders also include a ^*.tar or ^*.JSON file which can be opened in WebPlotDigitizer; the file contains the digitized profile and shows the digitization method. Graphical temperature data were digitized in the following way: the figure was loaded into WebPlotDigitizer as an xy plot. The axis of the figure was then calibrated by defining two points on each axis, with these definition points preferably being as far away from each other as possible. The program has two options for digitizing, namely “manual” or “automatic” extraction, and typically a combination of both methods was used. Using the automatic mode, profiles were defined by creating a mask, and an automatic extraction was used to produce an initial set of points. The manual mode was then used to adjust these points where necessary to fine-tune a visual match to the graphic. For figures that contained both in situ temperature points and a fitted temperature line, only the points were detected using WebPlotDigitizer's “blob detector” automatic algorithm.

Digitizing data from graphics can introduce non-trivial digitization uncertainty associated with the potentially improper alignment of the defining axes (systematic uncertainty) and/or the improper alignment of the data points (random uncertainty). The magnitude of this digitization uncertainty is proportional to the sizes of the graphic and data range, which vary from profile to profile. To highlight the potential importance of the digitization uncertainty, we propagate the uncertainty associated with a ±2 pixel error in each of the axes, thus defining points for both temperature and depth in the digitized profiles with the thinnest ice (Tuto_D-11; 48 m) and the thickest ice (DYE-3; 2038 m). The Tuto_D-11 graphic (Davis, 1967; Fig. 2a) is ∼1000 pixels tall and spans 47.5 m, which yields a ±2 pixel depth uncertainty of ±0.15 m. The Tuto_D-11 graphic is ∼600 pixels wide and spans 32 ^∘C, which yields a ±2 pixel temperature uncertainty of ±0.1 ^∘C. In comparison, the DYE-3 graphic (Gundestrup and Hansen, 1984; Fig. 2b) is ∼1200 pixels tall, spans 2200 m and yields a ±2 pixel depth uncertainty of ±3.6 m. The DYE-3 graphic, ∼900 pixels wide, spans 8 ^∘C and yields a ±2 pixel temperature uncertainty of ±0.02 ^∘C. These end-member scenarios highlight how digitization uncertainty varies from graphic to graphic.

Figure 2 (a) Graphic from which the Tuto_D-11 profile is digitized (Davis, 1967). This figure is reproduced from a U.S. government publication. (b) Graphic from which the DYE-3 profile is digitized (Gundestrup and Hansen, 1984). This figure is reproduced with permission from the International Glaciology Society.

All ice temperature profiles – from both digitized and tabulated sources – were standardized by interpolating between points, using cubic spline interpolation with no overshoot, in order to resolve the two common depth scales, namely the absolute depth scale at 1 m vertical resolution (Fig. 3a, c, e) and the normalized depth scale at non-dimensional 0.01 vertical resolution (Fig. 3b, d, f). The latter allows easy comparison between sites of different ice thicknesses and is useful for overcoming slight differences in ice thickness when comparing observed and modeled temperature profiles. The normalized depth-scale temperature file is named “temperature_dnorm.csv”. The absolute depth-scale temperature file, which is expressed as temperature with depth from the surface, is named “temperature.csv”. In each file, “NaN” means that no temperature data are available at the given depth, and in the absolute depth file, “−999” refers to an elevation below bedrock. Ice thicknesses vary significantly across the database profiles, from 3085 m at NGRIP to 48 m at Tuto_D-11, so the number of −999 below-bed null values varies between sites. This is because the temperature vector length is constant across the database entries, but the ice thickness varies across database entries.

Figure 3 Overview of all ice temperature profiles in the database expressed in both absolute depth (a, c, e) and normalized depth (b, c, f). For visibility, the profiles are divided into local ice caps (a, b), marginal ice sheet sites (c, d), and inland ice sheet sites (e, f). The local ice thickness coverage of each individual profile is given as a percentage.

3 Metadata

The database includes additional supplementary information. For each borehole, this information is stored in the file “meta.csv”, which contains a total of 18 metadata fields (Table 2). Every measurement entry in the database is labeled with a unique borehole ID, and a non-unique, alternative, and more descriptive place name. This alternative place name was deemed useful as, in some cases, ice temperature profiles have been measured over different campaigns by different people who have used different nomenclature for the same sites. Furthermore, borehole IDs are generally defined by either the group carrying out the measurements or by the site of the measurement. These borehole IDs are somewhat inconsistent, however, based on the differing conventions used by previously published studies.

Whenever specific locations are available, then that of each ice temperature borehole is provided in degrees east longitude and degrees north latitude. In some cases, locations have been shown in figures in the accompanying literature from which approximate locations were georeferenced. In addition to the position of the borehole, three further sets of information are provided relating to geographic setting. The site type describes whether the profile samples the Greenland ice sheet or a local ice cap. The geographic location provides a more descriptive general location of borehole sites using recognized regional geographic names. Last, the source of the location information is given.

The local ice thickness at the borehole location is provided in meters, along with the ice thickness source. When the original literature does not report the local ice thickness, an estimate is provided using BedMachine v3 (Morlighem et al., 2017). In these instances, the ice thickness is useful for determining the fraction of the ice depth, from the surface to the bed, covered by the given temperature profile and given as the percentage coverage in another field. Also listed are the depth in meters from the surface of the ice sheet to the uppermost temperature measurement of a given profile (depth of top) and the lowermost temperature measurement (depth of bottom).

For each profile, the data source describes the source of the temperature data, including links to open-access data repositories used by authors to share the original data, and the science source lists any original scientific publication that describes or mentions the temperature data. The related DOIs for both the data and science source categories are also reported when available. The date describes the common era year during which the temperature profile was measured. In the best cases, this includes the month and day that a given profile was measured. Finally, a field for additional notes is included, presenting any further information (relating to, for example, temperature, thickness, or location) that does not fit into the above categories.

Despite the best efforts of our community, there are some missing metadata fields in the database that are indicated with a question mark. The database is being continually updated. Any additional information, such as missing, updated, or new temperature data, may be added to the database by contacting the author team here or opening a GitHub issue at https://github.com/GEUS-Glaciology-and-Climate/greenland_ice_borehole_temperature_profiles/issues (last access: 1 March 2022). We use GitHub issues to resolve and identify and clarify problems in the ice temperature datasets, and we also welcome any new or missing data into the database via GitHub issues.

Table 2 Description of the database metadata structure.

4 Data summary

Of the 95 deep-ice temperature profiles making up the database, 70 originate from Greenland and 25 from the Canadian Arctic (Fig. 1). Of the 70 profiles from Greenland, 66 are from the ice sheet, and 4 are from peripheral glaciers and ice caps. All 29 peripheral glaciers and ice cap temperature profiles are plotted together in Fig. 3a and b. The peripheral glacier and ice cap profiles are relatively shallow; the deepest one reaches 420 m at a site where the ice thickness is estimated to be 520 m. Profiles from the ice sheet are divided into 50 marginal profiles (Fig. 3c and d) and 16 inland profiles (Fig. 3e and f), according to proximity to the main flow divide. The profile at GISP2, reaching bedrock at depth 3053 m and with a local ice thickness coverage of 98 %, is the deepest measurement in the database. Overall, 9 of the 95 profiles reach depths >1000 m.

The inland profiles also exhibit the coldest temperatures, with temperatures reaching a minimum of around −32 ^∘C halfway down the profile at normalized depth 0.48. In total, 34 of the database's 95 profiles are from locations characterized as “cold bedded”, where the ice appears to be frozen to the bedrock. This characterization is assigned based on the observed bottom values of the profile. At some sites, we cannot speculate on the basal thermal state based on temperature observation alone. A total of eight of these cold-bedded profiles extend to >95 % of the local ice thickness (i.e., database coverage >95 %). This contrasts with 56 “warm-bedded” sites where the ice at the bed is at pressure melting point. In total, 15 of these warm-bedded profiles extend to >95 % of local ice thickness. The basal thermal state of the three sites, CampIII_78, CampVI_50, and SiteII, cannot be characterized, as these boreholes covered less than 25 % of the local ice thickness. Overall, 39 of the 95 profiles have a depth coverage >90 %.

The temporal distribution of the profiles shows how the number of available ice temperature measurements has increased over time (Fig. 4). The temperature profiles span a period of 70 years. This implies that the shallow measurement points cannot necessarily be expected to reflect current ice temperatures. Arctic air temperatures have warmed markedly during the past 7 decades. The response time for deep-ice temperature variations, however, is usually much longer than 70 years (Cuffey and Paterson, 2010). For example, at a low-accumulation site in the ice sheet interior, a surface air temperature signal from ca. 50 years ago can likely influence present-day ice temperatures at >50 m below the ice sheet surface (Muto et al., 2011; Karlsson et al., 2020).

Figure 4 Histogram showing the temporal distribution of the temperature profile measurements, which are color coded based on the format in which the original data were available. The profiles entered into the database were available as a digitized graphic (orange), a digitized table (blue), or already published and digitally available (green).

The majority of the ice temperature profiles presented here have not been previously available in digital format. Profiles from only 31 boreholes (ca. 36 % of the database) have been previously published in open-access data repositories (Ryser et al., 2014; Harper, 2017; Doyle et al., 2018a; Seguinot et al., 2020; Hubbard et al., 2021b; Law et al., 2021). The public availability of these profiles highlights the laudable open-science trend in the natural sciences. While we present derived products of these profiles, we provide links to the open-access repositories in which the data were originally deposited and encourage database users to recognize and cite these previous data repository sources where applicable. Profiles from 54 boreholes (ca. 64 % of the database) are made digitally available here for the first time. This includes profiles from 24 boreholes that were provided as tabulated data by community members, profiles from 24 boreholes that were digitized from previously published graphics, and 6 that were previously available as published tables.

5 Sources of uncertainty

Borehole ice temperature measurements are associated with at least three main sources of uncertainty, including measurement uncertainty, depth uncertainty, and disturbance uncertainty. There are other sources of uncertainty that we do not consider here, such as self-heating of the sensor, sensor calibration, and leakage path, to name a few (e.g., Clow et al., 1996; Clow, 2008). Measurement uncertainty reflects the accuracy with which a temperature sensor can measure ice temperature. This uncertainty is influenced by both the design and calibration of the temperature sensor. It is reasonable to suggest that the temperature sensor accuracy has improved by at least an order of magnitude from 1950 to today (from perhaps ±0.1 ^∘C in the 1950s to an effective limit of better than ±0.01 ^∘C today). This trend primarily reflects improvements in the down-borehole technology deployed in glaciological applications and not the absolute accuracy of temperature sensors available for public purchase.

The issue of depth uncertainty is greatest in ice temperature profiles measured by manually deployed loggers. Manual loggers can have a characteristic depth uncertainty of ±1 %. By contrast, automated winch-deployed loggers typically use rotary encoders to measure logger depth with an uncertainty better than ±0.01 %. For example, the GISP2 borehole has a drill depth of 3053 m but a logger depth of only 3050 m, which represents a discrepancy of ±0.001 % (Cuffey et al., 1995). At deep interior sites, the temperature gradient in the lower part of the ice column is typically of the order of 0.02 K m⁻¹ (see Fig. 3e and f). At fast-moving sites, the temperature gradient may be even larger (see Fig. 3c and d), thus increasing the temperature uncertainty associated with depth uncertainty. At shallow sites, temperature gradients are typically small (see Fig. 3a and b), making the depth uncertainty negligible. For example, at Camp Century, the temperature profile was likely measured with a depth accuracy higher than ±0.1 %. With a measured basal temperature gradient of 17.5 ^∘C km⁻¹, this is equivalent to an ice temperature uncertainty of ±0.025 ^∘C at the maximum borehole depth of 1388 m. There are, however, some site settings – especially when using manually deployed loggers, logging non-vertical boreholes, or measuring extreme basal temperature gradients – in which depth uncertainty can actually play a larger role in ice temperature uncertainty than temperature sensor accuracy.

Disturbance uncertainty, which is associated with ice temperatures equilibrating after drilling, is a third source of measurement uncertainty. In boreholes filled with drilling fluid, the disturbance uncertainty can include the development of convection cells within the equilibrating drilling fluid. Convection cells are expected to occur when the temperature gradient in the ice exceeds a critical value dependent on the sum of a lapse rate term and the critical potential-temperature gradient (Clow, 2014). The size and strength of convection cells is a poorly defined function of drilling fluid density, borehole diameter, and ice temperature gradient. In the Agassiz77 borehole, slight thermal instabilities in the drilling fluid produced temperature offsets of ±0.05 ^∘C around the mean ice temperature gradient. In the NEEM 2011 profile, however, the convection cells between 1500 and 2000 m depth reach ca. 50 m vertical length scale and produce ±0.5 ^∘C offsets around the mean ice temperature gradient (Fig. 5). The GISP-D temperature log (not presented here) also shows pronounced laminar convection cells between 1600 and 2000 m depth (Clow et al., 1996).

Figure 5 Offsets from the mean ice temperature gradient highlight the influence of large (ca. 50 m) convection cells forming in the borehole fluid in a deep section of the NEEM ice temperature profile measured in 2011. These borehole fluid convection cells contribute to the disturbance uncertainty associated with an equilibrating borehole temperature profile when the borehole-logging technique is used. This uncertainty is avoided when sensors or cables are frozen into the ice in other measurement techniques.

Given the ultra-low density of air, air-filled boreholes cannot support analogous convection cells. The disturbance uncertainty associated with convection cells is accordingly negligible in air-filled boreholes. However, temperature disturbances due to changes in atmospheric pressure have been observed in air-filled bore holes in Antarctica. This effect is largest in the top 20 m of the borehole, just below the ice surface, with magnitudes on the order of up to ±0.003 ^∘C (Clow et al., 1996). In water-filled boreholes, such as those created by hot-water drilling, the disturbance uncertainty is associated with the sensible heat released from hot drilling water and the latent heat released by borehole refreezing. This thermal disturbance of hot-water drilling often requires many months to fully dissipate, as the latent energy released by the refreezing borehole warms the surrounding ice. Empirical cooling curves of temperature measurements collected multiple times post-drilling can usually constrain this thermal disturbance effect to better than ±0.1 ^∘C (Humphrey and Echelmeyer, 1990; Ryser et al., 2014; McDowell et al., 2021). Mechanically drilled ice core holes also experience thermal disturbance, resulting primarily from the sensible heat carried into the hole by the drill fluid. The associated disturbance uncertainty is, however, much smaller than for hot-water drilled holes (Clow, 2015). The temperature disturbance caused by mechanical drilling with fluid-filled boreholes dissipates to the level of precision within 5 d (Zagorodnov et al., 2012).

The total measurement uncertainty varies substantially between profiles in this database. Generally, these sources of uncertainty can be viewed as being independent from each other and therefore summed to estimate the total measurement uncertainty for a given ice temperature profile. Assessing this total measurement uncertainty is most critical for studies focused on absolute temperature, for example, when remeasuring present-day ice temperatures to compare with past ice temperatures at a specific site. We provide guidance for assessing the total measurement uncertainty by contrasting the uncertainty budgets of an older (Tuto_D-11; 1962) and younger (Store_R30; 2019) ice temperature profile (Table 3). The Tuto_D-11 ice temperature profile is from an air-filled borehole in which the ice temperature was measured with ca. 60-year-old sensor technology using a manually deployed logger. The Tuto_D-11 data also include an additional graphical digitization uncertainty (described in Sect. 2). The Store_R30 ice temperature profile was measured using a fiber optic string frozen into a hot-water drilled borehole that was continuously monitored to establish empirical cooling curves. The contrasting uncertainty budgets of these two boreholes (see Table 3) highlights the value of consulting the original science source publication for site-specific uncertainty assessments.

Table 3 Contrasting total measurement uncertainty budgets for the Tuto_D-11 and Store_R30 ice temperature profiles. The conversion of depth uncertainty into a temperature uncertainty is done, assuming a characteristic temperature geothermal gradient of 20 ^∘C km⁻¹. Here, the word “Negligible” means that the uncertainty is <0.01 ^∘C.

6 Comparison with simulated ice temperatures

Numerical ice flow models are critical tools for projecting the future sea level rise contribution from the Greenland ice sheet. While the simulated initial states are typically evaluated against observed ice sheet form (i.e., thickness) and flow (i.e., velocity), the relative paucity of ice temperature data means that they are typically not evaluated against the observed ice sheet thermal state (Goelzer et al., 2020). Given the sensitivity of ice viscosity to ice temperature, biases in the thermal state of ice sheet simulations likely contribute to biases between simulated and observed recent Greenland ice loss (Goelzer et al., 2017; Aschwanden et al., 2019; Law et al., 2021). Identifying and addressing these biases in thermal state is critical for improving projections of Greenland ice loss (Aschwanden et al., 2021). To help shape a best practice for using the database that we present, we describe an approach for evaluating a contemporary ice sheet simulation against all observed Greenland deep-ice temperature observations.

Several process level studies of potential heat sources have been performed, which compare individual observed temperature profiles from local areas with temperature profiles modeled by a thermal, or thermomechanical, ice flow model (Iken et al., 1993; Lüthi et al., 2002; Harrington et al., 2015; Lüthi et al., 2015; Meierbachtol et al., 2015; McDowell et al., 2021; Law et al., 2021; Maguire et al., 2021). Although these studies featured different local areas, the comparisons generally showed that models tend to underestimate englacial temperatures and thus need to incorporate additional heat sources in order to reproduce the observed ice temperature profiles. Suggested additional heat sources include cryo-hydrological warming – which transfers latent heat when surface meltwater flows through englacial pathways and refreezes – deformational heating, and basal water heat flux (Funk et al., 1994; Wohlleben et al., 2009; Phillips et al., 2013; Lüthi et al., 2015; Zekollari et al., 2017; Karlsson et al., 2020). Modeled ice temperatures are likely also influenced by the choice of geothermal heat flow map, however, there is a diversity of opinion regarding the magnitude and spatial distribution of geothermal heat flow beneath the ice sheet (Rezvanbehbahani et al., 2017; Colgan et al., 2022).

The contemporary thermomechanical ice sheet simulation we adopt is the Parallel Ice Sheet Model (PISM; Bueler and Brown, 2009) simulation of Aschwanden et al. (2016). This simulation represents the form, flow, and thermal state of the Greenland ice sheet in ca. 1990, following a 125 kyr paleoclimatic spin-up and followed by another 2 kyr of transient equilibrium, with mass flux adjustment forcing to minimize the misfit against the observed ice sheet thickness and extent. This ice sheet simulation has a horizontal resolution of 900 m and a vertical resolution of 20 m. The geothermal heat flow map by Shapiro and Ritzwoller (2004), variable in space but not time, was used as a basal thermal boundary condition. PISM uses an enthalpy scheme for the conservation of energy calculation, in order to accommodate heat transfer in both freezing and temperate ice (Aschwanden et al., 2012). While the spin-up is meant to approximate the Greenland ice sheet at a specific time slice (ca. 1990), we are comparing this simulated thermal state with the temperature profiles observed over a 70-year time span. Shallow profiles located close to the margin can experience significant changes in ice thickness and temperature on this timescale.

At locations of observed temperature profiles, modeled vertical temperature profiles were extracted from the PISM simulation based on the nearest-neighbor grid point. Only temperature profiles from the Greenland ice sheet were included in this analysis. The PISM vertical temperature profiles were transformed from the height above bed to the normalized depth below surface and linearly interpolated to the vertical resolution of the normalized depth field of the database. Furthermore, at a few drill site locations, the modeled ice was too thin and did not have enough vertical grid points to interpolate the modeled temperatures to a proper normalized depth axis. It was therefore necessary to exclude profiles from the analysis where the simulated ice thickness was fewer than five vertical grid points (50 m). Additionally, in cases of temporally repeated observed profiles, only the most recent measurement was included in this analysis (see Table 4).

Out of the 95 temperature profiles from the temperature database, 66 were ultimately suitable for inclusion in this comparative analysis of observed and modeled ice temperatures. The observed temperature profiles were divided based on three characteristic regimes. First, whether the profiles were located in the accumulation area (16 profiles) or ablation area (38 profiles). Second, whether the basal thermal state of the profiles was considered warm bedded (44 profiles) or cold bedded (9 profiles). Finally, whether profiles were located in high-strain regions (10 profiles) or low-strain regions (46 profiles). The surface mass balance regime is determined by whether the borehole is located below the snow line, in the ablation area, or above the snow line, in the accumulation area, in contemporary satellite imagery (Fig. 1). The basal thermal state regime is based on whether the ice–bed interface is measured to be below the pressure melting point temperature (i.e., cold bedded) or not (i.e.,  warm bedded). In instances where the borehole does not reach the bed, we extrapolate the basal thermal state where reasonable (i.e., FladeIsblink06 is likely cold bedded), or we list the basal thermal state as being “unknown” where the extrapolation distance seems unreasonable (i.e., CampVI_50 is unknown). Finally, the ice dynamic regime is classified as being high strain when the ice flow is channelized and low strain when sites are located in sheet flow or divide flow. Figure 6 shows the observed temperature profiles for each of the three regimes.

Table 4 Overview of the number of profiles in the three regimes before and after excluding profiles not usable for the model comparison analysis.

Figure 6 Observed temperature profiles on a normalized depth scale for three characteristic regimes. (a) Surface mass balance regime with profiles located in the ablation area (pink) or the accumulation area (gray). (b) Basal thermal state regime with profiles either cold bedded (purple) or warm bedded (green). (c) Ice dynamic regime with profiles located either in high-strain-rate regions (orange) or low-strain-rate regions (blue).

The difference between the modeled and observed profiles was calculated at 0.01 normalized vertical resolution and then aggregated into coarser vertical layers at 0.1 normalized vertical resolution (Fig. 7). This method of examining profiles in vertical sections provides better insight to the depth variations in model bias than reducing every profile into a single metric of best fit. This vertical box plot analysis can also provide some insight with respect to which terms in the energy balance – advection, diffusion, heat production, or latent heat transfer – may be responsible for the differences between modeled and the observed profiles, as the relative importance of these terms varies with depth. Where our comparison box plots show negative temperature differences, modeled temperature values are smaller than the observed temperatures, implying that this PISM run is underestimating temperatures. The reverse occurs where box plot temperature differences are larger than 0 ^∘C.

Figure 7 Each plot shows the difference between the modeled minus observed normalized depth profiles aggregated into coarser vertical layers of 0.1 resolution. (a, d) Surface mass balance regime. (b, e) Basal thermal state regime. (c, f) Ice dynamic strain regime. The N value shows the number of profiles that the box plot is based on. Points plotted on top of each box plot show individual temperature difference values for all of the profiles at that depth within that category.

The surface mass balance regime comparison generally highlights a better PISM fit to observations in the accumulation area than the ablation zone (Fig. 7). This is generally expected, since the accumulation area more closely represents the theoretical steady-state temperature profiles upon which classic ice flow heat equations are predicated. The accumulation box plot shows that the mean difference value shifts from being positive in the upper part of the profile to negative, with a greater spread in differences, in the bottom 20 % of the profile. In comparison, both the ablation and temperate box plots, which share strong overlap in comprising profiles, are clearly less well captured by PISM. These warmer ablation area profiles show a larger spread in difference values throughout the entire vertical profile. Importantly, both box plots show a systematic negative mean difference at all depths, meaning that the simulation consistently underestimates ice temperatures. This may suggest that one or more heat sources are being poorly represented in warmer ablation areas in this particular PISM simulation. This is a general problem for solving the heat equation in all ice flow models. Previous studies similarly showed poor representation of the ablation area and the need to include additional heating sources to reproduce the higher observed englacial temperatures (McDowell et al., 2021).

After excluding profiles for which the modeled thickness is less than 50 m, the cold-bedded regime box plot represents the smallest sample size of just five profiles (Camp Century, Dye-3, GISP2, GRIP, and Station Central). These profiles are all drilled at inland locations where the ice is moving relatively slowly. This generally limits the role of horizontal heat advection at these sites. In contrast to the warm-bedded and ablation box plots, the cold-bedded box plot shows a systematic positive mean difference over all vertical layers, meaning that the modeled temperatures generally overestimate the observed temperatures over the entire column of ice. The largest differences between modeled and observed temperatures (+4 ^∘C on average), however, are found near the ice sheet bed. Here, where vertical ice advection becomes small, the diffusion of geothermal flow should be the most influential term of the heat equation. As heat diffusion is generally well-represented and well-constrained in all ice flow models, this may suggest that this particular PISM simulation employs too high a geothermal heat flow into the ice sheet base. We note that the geothermal flow is poorly constrained under the Greenland ice sheet (Rezvanbehbahani et al., 2017; Colgan et al., 2022) and over long spin-up periods, so even small deviations in the geothermal heat flow can result in large variations in ice temperature.

The shape of temperature profiles in fast-flowing, high-strain regions is expected to have a C shape associated with strong horizontal advection of upstream ice temperatures (Cuffey and Paterson, 2010; Iken et al., 1993). While this advection of relatively cold upstream ice cools the middle depths of the ice temperature profile, in comparison to the surface boundary condition, the deformational heating warms the lower depths (Cuffey and Paterson, 2010). The high-strain box plot shows that the PISM simulation appears to systematically underestimate ice temperatures in the lowermost 50 % of the ice sheet. This bias reaches up to −5 ^∘C, on average, between 80 % and 90 % depth. This underestimation of ice temperatures near the bed may result from this particular PISM simulation underestimating the strain heating associated with enhanced ice deformation in softer Last Glacial Period ice (e.g., Lüthi et al., 2002; Karlsson et al., 2013). Generally, however, it is challenging to simulate temperature profiles in high-strain ice flow regimes, as this requires correct parameterization of both local and upstream heat budgets, as well as the horizontal advection profile linking these heat budgets, in order to correctly represent steep gradients in the ice temperature profile (Cuffey and Paterson, 2010).

The low-strain box plot shows a systematic negative mean difference value across all ice depths. This suggests that, in slow-flowing regions, the model generally underestimates ice temperatures. This might be explained by this particular PISM simulation employing a spin-up surface boundary condition that had either an air temperature that is too cold or an accumulation rate that is too high. As vertical ice advection is dependent on the accumulation rate, overestimated accumulation rates can result in overestimated vertical advection rates, which enhance the advection of relatively cold near-surface ice deeper into the ice sheet. The low-strain regime profiles are diverse, however, and comprise both ablation and accumulation area profiles, with both cold and warm beds. The low-strain regime generally has the largest spread of difference values. This makes it difficult, in comparison to the other box plots presented here, to interpret a process level mechanism for the systematic cold bias across these diverse sites.

In this modeled run, the key parameters influencing the modeled temperatures are the enthalpy, which requires understanding of the liquid water fraction in the ice column and its evolution in time; the geothermal heat flow, which also has the potential to affect the presence of water in a given location; and the surface boundary conditions, i.e., accumulation rate, which control vertical temperature advection, and ice surface temperatures. These parameters are all very difficult to constrain; however, this is especially the case for the enthalpy and the geothermal heat flow parameter.

7 Outlook

We have compiled 95 ice temperature profiles from 85 individual boreholes from Greenland and the Canadian Arctic into a comprehensive and consistent community resource. This constitutes an increase of +54 ice temperature profiles over what was previously available in open-access data repositories (64 % of the database). The scientific community is encouraged to contribute missing or new data or metadata to the database either by contacting the author team or using the GitHub issue pathway. While we expect the primary use of this database to be evaluating simulated 3-D ice temperature fields from thermomechanical ice flow models, there may soon be a sufficient number of ice temperature profiles from the Greenland ice sheet and peripheral ice caps to use machine learning approaches to spatially interpolate ice temperatures between boreholes on the basis of independent geophysical variables. This database may possibly serve as a modest first step towards addressing this emerging challenge, and other challenges, associated with understanding the thermal state of ice in Greenland and the Canadian Arctic.

We provide an initial comparison between observed ice temperature profiles and those simulated by the PISM ice flow model. This comparison highlights the inherent challenge of comparing 1-D observational datasets with an uneven spatial distribution to a simulated 3-D temperature field. These challenges include the uneven ice thickness between modeled and observed geometry, especially at marginal sites; resolving the difference between the systematic biases and random error in the modeled temperature profiles; and confronting the fundamental uncertainty in the basal thermal state due to complex processes. In this comparison, the poor representation of cryo-hydrological warming, or the transfer of latent energy, appears to be contributing to a model cold bias at peripheral ice sheet sites. Similarly, the poor representation of deformational heating and vertical ice stratigraphy appears to be contributing to a model cold bias at high-strain-rate sites. Clearly, thermomechanical models provide the most reliable solution to the heat equation in low-strain, cold-bedded inland areas.