Potential to recover a record of Holocene climate and sea ice from Müller Ice Cap, Canada

Abstract

Müller Ice Cap sits on Umingmat Nunaat (Axel Heiberg Island), Nunavut, Canada, ~ 80°N. Its high latitude and elevation suggest it experiences relatively little melt and preserves an undisturbed paleoclimate record. Here, we present a suite of field measurements, complemented by remote sensing, that constrain the ice thickness, accumulation rate, temperature, ice-flow velocity, and surface-elevation change of Müller Ice Cap. These measurements show that some areas near the top of the ice cap are more than 600 m thick, have nearly stable surface elevation, and flow slowly, making them good candidates for an ice core. The current mean annual surface temperature is −19.6 °C, which combined with modeling of the temperature profile indicates that the ice is frozen to the bed. Modeling of the depth-age scale indicates that Pleistocene ice is likely to exist with measurable resolution (300–1000 yr m⁻¹) 20–90 m from the bed, assuming that Müller Ice Cap survived the Holocene Climatic Optimum with substantial ice thickness (~400 m or more). These conditions suggest that an undisturbed Holocene climate record could likely be recovered from Müller Ice Cap. We suggest 91.795°W, 79.874°N as the most promising drill site.

Introduction

A high-resolution paleoclimate record from as far north as possible could provide valuable information about past meridional variations in climate, particularly polar amplification (Bekryaev and others, 2010). Similarly, a proxy record of sea-ice from the far north would fill a crucial gap in knowledge of past conditions since direct records span only decades for photo records (e.g. Thomson and others, 2011) to hundreds of years for written and oral records (Curley and others, 2021). Ice cores can provide longer climate records, since glacier ice holds proxies for temperature (e.g. δ¹⁸O or δD in the ice itself) while also trapping air bubbles that provide samples of past atmospheric composition (e.g. Fisher and others, 1995; Kinnard and others, 2006). Moreover, proxies have recently been developed that are sensitive to proximal or regional sea-ice conditions (Vallelonga and others, 2021), so an ice core near the Arctic Ocean could provide information about regional sea-ice and climate conditions in the past.

Preservation of a usable climate record in glacier ice only happens under favorable conditions. Cold surface temperatures are a prerequisite to preserve past conditions in stratigraphic order; melt moves material between layers, obscuring ages and affecting gas retention (Koerner and Fisher, 1990; Orsi and others, 2015). Slow ice flow is important to ensure that the paleoclimatic record of the core reflects changes in climate rather than in the source location and elevation of the ice (e.g. Fudge and others, 2019; Gerber and others, 2021).

While cold conditions and slow flow are effectively prerequisites for a good paleoclimate record, other competing effects must be balanced based on project-specific objectives. For example, there is typically a tradeoff between time span and age resolution of a core (Fischer and others, 2013; Van Liefferinge and Pattyn, 2013; Obase and others, 2023). Even in extremely cold areas in Antarctica, basal melt is possible at moderate geothermal heat fluxes if the ice is sufficiently thick and accumulation is sufficiently low so that the cold is not advected or diffused to the base too quickly (e.g. EPICA community members, 2004). Thick ice can then lead to a shorter time span in the core compared to thinner but otherwise similar areas where no basal melt takes place (Lilien and others, 2021). Conversely, large accumulation rates and basal melt can lead to excellent temporal resolution over a limited time span (Buizert and others, 2015).

A handful of ice masses in Canada, Greenland, Svalbard, and Russia extend north of 80°N (Pfeffer and others, 2014), and some of these may be suitable for obtaining a reliable paleoclimatic record. In Canada, on Umingmak Nuna (Ellesmere Island) the Northern Ellesmere Ice Fields reach 82°N and Agassiz Ice Cap reaches 81°N, both with surface elevation up to 2000 m a.s.l., while Müller Ice Cap on Umingmat Nunaat reaches 80.5°N and with surface elevation up to 1900 m a.s.l. To the west, Meighen Ice Cap reaches 80°N but its surface elevations reach only 270 m a.s.l. There are many small peripheral ice caps in Greenland that extend beyond 80°N, the largest of which are Hans Tausen Ice Cap (up to 83°N, 1400 m a.s.l.) and Flade Isblink (up to 81.8°N, 750 m a.s.l.). In Svalbard, Vestfonna and Austfonna on Nordaustlandet reach 80.2°N, with surface elevations up to 600 and 800 m a.s.l., respectively. In Russia, the Franz Joseph Land archipelago extends to 81.8°N with several ice caps, but none exceed 600 m a.s.l., while the Severnya Zemlya archipelago reaches 81.2°N, with several ice caps with surface elevations up to 800 m a.s.l. All ice-cap extents and elevations herein are from the Randolph Glacier Inventory (Pfeffer and others, 2014).

Most cores drilled on arctic ice caps have found short records or low-resolution climate information. For example, a 724-m core drilled on Akademii Nauk (80.5°N, 800 m a.s.l.), Severnaya Zemlya, Russia, found only a ~2500-year paleoclimatic record (Fritzsche and others, 2005). High accumulation rates allowed precise dating of only 1100 years of that core (Opel and others, 2013). Moreover, the low surface elevation of the site led to substantial refrozen surface meltwater in the core despite the high latitude. Other ice caps in Svalbard and the Russian Arctic have similar or lower surface elevations (Porter and others, 2018), so any cores from such areas are also likely to be plagued by substantial melt (Koerner and Fisher, 2002). Other cores have been recovered at higher elevations from northerly ice caps, such as a 345 m core drilled at 82.5°N, 1270 m a.s.l. on Hans Tausen Ice Cap, Greenland (Hammer and others, 2001), a 425 m core from Flade Isblink (81.3°N, 630 m a.s.l.), Greenland (Lemark and Dahl-Jensen, 2010), and 100–300 m cores from Agassiz Ice Cap (80.7°N, up to 1730 m a.s.l.), Canada (Fisher and others, 1995). These cores have obtained records with less melt but with lower temporal resolution or shorter duration. The oldest ice at Hans Tausen Ice Cap dated to about 11 ka, suggesting that the ice cap may have nearly disappeared during the Holocene Climatic Optimum (HCO, 8−4 ka; Madsen and Thorsteinsson, 2001), and the oldest ice at Flade Isblink was modeled to reach only 4 ka (Lemark and Dahl-Jensen, 2010). The Agassiz cores recovered Pleistocene ice (more than 11.7 ka; Vinther and others, 2008). Existing cores thus hint that northern Canadian ice caps have been more stable than the peripheral ice caps in Greenland and experience less melt than those in the Russian Arctic or Svalbard. This difference may result from isostatic rebound relating to changes in the Greenland Ice Sheet, residual ice remaining from the Innuitian and Laurentide Ice Sheets in northern Canada, or differences in melt or accumulation dominating surface mass balance (SMB; Koerner and Fisher, 2002). While the cause of this geographic variability is not firmly established, available records suggest that northern Canadian ice caps are the most promising target for recovering an extended, far-northerly paleoclimate record.

Study site

Müller Ice Cap on Umingmat Nunaat, Nunavut, Canada, extends beyond 80°N with a maximum surface elevation of about 1900 m a.s.l. (Porter and others, 2018). This northerly latitude and elevation are almost unique globally, as only the smaller Northern Ellesmere Ice Fields reach similar elevations at comparable latitude (Pfeffer and others, 2014). Müller Ice Cap covers an area of approximately 4200 km² including its outlets (Pfeffer and others, 2014), some of which terminate on land to the east while others reach the ocean or terminate in low-elevation fjords after flowing through the Prince Margaret Range to the southwest (Fig. 1). The ice cap is thought to have existed for more than 100 kyr and to have been part of the much larger Innuitian Ice Sheet at the Last Glacial Maximum (LGM; Dalton and others, 2022). For an extended period ~ 25 ka, the present-day location of Müller Ice Cap was likely a local ice dome with flow toward ice streams in present-day Massey, Eureka, and Nansen sounds (England and others, 2006; Dalton and others, 2020). Radiocarbon ages and geomorphology in the area suggest that retreat was underway by 10.5 ka (summarized in Bednarski, 1998; Atkinson, 2003), and there is overall agreement that the ice margin retreated from the edge of Umingmat Nunaat to near its present position from 10 to 8 ka (Dalton and others, 2020). Since evidence of a more retreated margin would have been destroyed by re-advance, these geological data do not exclude the possibility of further retreat or even absence of Müller Ice Cap during the HCO.

Figure 1. Overview of study area. (a) Location of (b) within Canada (red box). (b) Sentinel-2 composite of Umingmat Nunaat, courtesy of the European Space Agency. Yellow box shows area plotted in (c) and subsequent figures. (c) Surface elevation of Müller Ice Cap and surroundings from ArcticDEM (Porter and others, 2018). Thick black lines show ice-covered area. Thin black lines show contours of surface elevation. Black star indicates the most promising ice-core site. Red box indicates location of (d). (d) Zoom-in of area ~ the most promising ice-core site (black star). Markers show locations where measurements of a snow pit, firn cores, and weather were taken. Background as in (b).

Müller Ice Cap underwent substantial investigation in the 1950s and 1960s in a program led by McGill University (Müller-Battle, 1960). By digging a shaft that extended through 41 years of firn, they found that SMB was 370 kg m⁻² yr⁻¹ at the top of the ice cap. This can be compared to the meager 58 kg m⁻² yr⁻¹ accumulation at Eureka to the east (Cogley and others, 1996), though imprecise dating of the firn shaft suggests that the ice-cap SMB has considerable uncertainty and Cogley and others, 1996, considered 370 kg m⁻² yr⁻¹ anomalous. Shallow coring atop the ice cap in 1974, measuring 12-year average SMB down to a clear melt layer from summer 1962, indicates that SMB is much lower, approximately 160 kg m⁻² yr⁻¹ (Koerner, 1979). From melt layers within the firn, that work found that melt during the 12-year period ending with 1974 averaged 5.4 kg m⁻² yr⁻¹, or ~3% of the accumulation rate.

Since 1959, a mass balance record has been maintained just south of Müller Ice Cap on White Glacier (Fig. 1) (Cogley and others, 2011; Thomson and others, 2017). While the overall SMB of White Glacier has been negative (~ −200 kg m⁻² yr⁻¹ over the whole glacier), the glacier saw thickening above 1400 m between 1960 and 2014 (Thomson and others, 2017), suggesting that the higher-elevation Müller ice cap may have experienced positive SMB as well. However, overall mass balance at White Glacier has become increasingly negative, going from −0.1 m w.e. yr⁻¹ during 1960–1970 to −0.3 m w.e. yr⁻¹ during 2012–2022 (Thomson and others, 2017; WGMS, 2024). The equilibrium line altitude rose from 1000 to 1270 m during a similar time period (Thomson and others, 2017). However, glacier extents around Umingmat Nunaat decreased by <1% from 1959 to 2010, with advance or minor retreat on Müller Ice Cap's outlets (Thomson and others, 2011), suggesting that this White Glacier record does not translate directly to conditions on the ice cap.

Remote-sensing datasets provide additional context for glaciological conditions on Müller Ice Cap. Several velocity products are available for the area, all based on feature/speckle tracking of optical or radar imagery (Van Wychen and others, 2014, 2020; Gardner and others, 2019; Millan and others, 2022; Zinck and Grinsted, 2022).

While accurate for determining flow speeds of the outlets of the ice cap, these datasets have low relative precision in slow flowing areas and thus do not provide useful information about the flow of the ice cap's interior. Ice-penetrating radar surveys have been flown by NASA's Operation IceBridge and by the British Antarctic Survey (BAS) using CReSIS's Multichannel Coherent Radar Depth Sounder (MCoRDS) and the UTIG-JPL High-Capability Radar Sounder (HiCaRS), respectively (Paden and others, 2018; Benham and others, 2020). These surveys provide ice thickness on lines approximately 5–10 km apart, but in all freely available echograms any returns from internal layers are overshadowed by the surface multiple and other clutter. Zinck and Grinsted (2022) attempted to infer ice thickness in more detail and estimated that some areas have more than 750 m of ice, though the inferred distribution varied substantially depending on the method used; while that work provides a general picture of where thick ice might exist, it likely misses important details given the steeply varying topography that has been observed in airborne radar data. The global dataset of Millan and others (2022) provides another estimate of ice thickness with similar limitations. Available data do not establish whether the ice cap has well-layered ice or how ice thickness varies on the kilometer scale.

Here, we complement the existing measurements of Müller Ice Cap with ice-penetrating radar surveys, measurements of accumulation and temperature, and remote sensing of velocity and surface-elevation change. We use these data to identify the site on the ice cap that is most likely to retain old ice. Simple models of temperature and vertical ice flow are then used to estimate the basal temperature and age at that location to determine whether a useful paleoclimate record is likely to be preserved.

Field and remote-sensing data

Ice-flow speed

We calculated ice flow speeds using interferometric processing of synthetic aperture radar (InSAR) data acquired by Sentinel 1. As opposed to existing datasets in the region (Short and Gray, 2005; Van Wychen and others, 2016; Gardner and others, 2019; Zinck and Grinsted, 2022), which are well suited to study fast-flowing regions and their variability (Lauzon and others, 2023, 2024), the velocities calculated using InSAR have relatively high precision in the slow-flowing areas which might be suitable for an ice core.

Sentinel 1a and 1b data over Mueller's Ice Cap were interferometrically processed using the European Space Agency's SNAP toolbox. We processed ascending and descending single-look complex images acquired in the sensors’ interferometric wide swath mode during 2016–2022 to produce 6- or 12-day interferograms. Processing followed standard InSAR methodology, involving orbit correction, coregistration, debursting, topographic phase removal, and Goldstein phase filtering. All data that extended across the ice cap to bedrock on both sides were processed, resulting in a set of 28 ascending and 40 descending interferograms on four different ground tracks, all acquired in winters 2016–2022. The interferograms were unwrapped using the Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping (Chen and Zebker, 2002). Unwrapped results were then converted to radar line-of-sight (LOS) displacement and geocoded using the SNAP toolbox, and then converted to LOS speed per year.

LOS speed in multiple directions was corrected for phase ambiguity and then converted to ice-flow velocity. The oblique LOS speed was first projected to the horizontal plane and then the average speed over all bedrock (determined using a land cover database (Commission for Environmental Cooperation and others, 2020)) was removed to account for phase ambiguity and the average bias of the measurement. After this correction, RMS bedrock speeds for each scene vary from 0.08 to 0.72 m yr⁻¹. All measurements for each ground track were combined using a weighted average to produce mean horizontal LOS speeds in the four radar sight directions. The weighted mean products have RMS bedrock speeds of 0.17 to 0.40 m yr⁻¹. The LOS speeds were combined following Anderson and others (2020) to determine the two components of the flow vector in an orthogonal coordinate system (Statistics Canada Lambert; EPSG:3348).

InSAR processing shows that the outlet glaciers on the northeast side of Müller Ice Cap have a complex network of tributaries that extends deep into the interior of the ice cap (Fig. 2a). Nearly the whole ice cap drains to the northeast rather than through the Princess Margaret Range and down the large outlets to the southwest. There are several regions of slow flow (~2 m yr⁻¹) between the tributaries of the outlet glaciers and a few small regions of slow flow along the divide just northeast of the Princess Margaret Range.

Figure 2. Ice-flow speeds (a) from InSAR processing of Sentinel-1 data from 2015 to 2022 and (b) ITS_LIVE (Gardner and others, 2019) data from optical feature tracking. Black star indicates most promising drill site. Background images from Sentinel-1.

Of the previously available velocity measurements, ITS_LIVE (Gardner and others, 2019) best captures this complexity (Fig. 2b), yet in regions moving less than 5 m yr⁻¹ it misses the fine structure and misplaces the ice divide. This difference does not imply that previous datasets were low quality, simply that they chose methods suitable for the fast-flowing regions in which they were interested.

Surface-elevation change

We calculated rates of surface-elevation change from ICESat-2 and from ArcticDEM strip products, derived from optical satellites. These two datasets are complementary; though ICESat-2 has greater precision, it has lower spatial resolution, while the ArcticDEM strips have lower precision but greater spatial resolution and span a longer time window.

We used all available ArcticDEM data over the area to calculate the surface-elevation trend (Porter and others, 2018). These data consist of strips, formed from individual stereo pairs of WordView/GeoEye imagery, and span 2012–2021. The tiled ArcticDEM 2-m product was used as a reference DEM. Each strip was individually co-registered to the reference DEM over bedrock areas following Nuth and Kääb (2011). Coregistration used the demcoreg software package (Shean and others, 2016) and was necessary to remove small biases in the strips that can prevent meaningful inference of a trend. We excluded any strips with fewer than 2 × 10⁶ valid pixels, or where the mean absolute difference compared to the reference DEM was greater than 10 m. Within the remaining strips, bad pixels remained (due to clouds, poor correlation between the stereo pairs, or other errors in DEM creation); these were filtered out by removing all pixels where the difference compared to the reference DEM was more than three standard deviations, more than 100 m, or where the surface slope exceeded 40°. Differences at seams between neighboring tiles in the tiled ArcticDEM product led us to then create a new reference DEM from the median of these filtered strips.

The above steps were then repeated, using this new seamless DEM as the reference product. The strips were then divided into six-month bins; for each bin, we created a median mosaic with 16-m posting. Finally, we calculated the surface-elevation trend over 2012–2021 using a Theil-Sen estimator (Theil, 1950). The final product was then median sampled to 256-m resolution. Uncertainty on the ArcticDEM-derived surface-elevation change was estimated using the RMS elevation trend over bedrock. These results show relatively stable surface elevations near the center of the ice cap and thinning near its margins (Fig. 3a).

Figure 3. Rate of surface elevation change from (a) ArcticDEM and (b) ICESat-2. ArcticDEM changes were calculated from strips spanning 2012–2021 (Porter and others, 2018) while ICESat-2 changes were calculated from ATL11 from 2018–2022 (Smith and others, 2020). Black star indicates most promising drill site.

To determine thinning rates from ICESat-2, we simply calculated the average change from the ATL15 gridded product and used the published uncertainty estimate (Smith and others, 2020, 2022). ICESat-2 data for 2019–2022 also indicate that the outlets of the ice cap have changed significantly while the center of the ice cap has remained stable (Fig. 3b). Thinning rates on the land-terminating eastern outlets exceed 1 m yr⁻¹ in many places and are generally higher than those derived from ArcticDEM near the margins (Fig. 3). Since errors are 0.06 m yr⁻¹ for ICESat-2 and 0.08 m yr⁻¹ for the ArcticDEM-derived thinning, this 0.5 m yr⁻¹ increase in thinning represents a significant change. The pattern of thinning and thickening over the western outlets is like that from ArcticDEM, though significant exposed rock surrounding narrow glacier tongues make interpretation of the limited-resolution ICESat-2 data in this area problematic. Near the top of the ice cap, the ICESat-2 data show −0.009 ± 0.04 m yr⁻¹ of elevation change over 2019–2022, in good agreement with ArcticDEM and indistinguishable from zero.

Ice thickness and stratigraphy

We measured ice thickness and internal stratigraphy using ground-based and airborne ice-penetrating radar. Two surveys were flown by the Alfred Wegener Institute (AWI), Helmholtz Centre for Polar and Marine Research's Polar 5 aircraft (Alfred-Wegener-Institut, 2016) on May 8^th and 9^th, 2023. The flights followed nearly the same route (Fig. 4a) but used different bandwidths. The first operated at 150–520 MHz and the second used 180–210 MHz. The wider bandwidth provided higher resolution (~0.35 m vs ~4.3 m), while the narrower bandwidth allowed rapid processing to ensure that reliable data had been acquired. General system characteristics can be found in Hale and others (2016) and Franke and others (2022). The system used 8 channels, each co-polarized with the direction of travel, though during portions of the narrow-band flight only 7 channels were operational. All data were pulse compressed and SAR processed using the CRESIS toolbox (CReSIS Toolbox https://github.com/CReSIS/cresis-toolbox). While SAR processing corrects for off-nadir reflections along track, across-track reflections remain after processing, complicating identification of the bed reflection. To trace the bed through the clutter introduced by the complex topography, we had three individuals independently trace the bed. All layer tracing was done in ImpDAR (Lilien and others, 2020). Example profiles from the AWI surveys are shown in Supplementary Figure 1.

Figure 4. Ice-penetrating radar survey overview. (a) Survey map. Lines indicate tracks, with color indicating the survey. (b) Ice thickness from all surveys. Colors show ice thickness. Red box shows area of Figure 5. Black contours show 5 m yr⁻¹ velocity from Sentinel-1 interferometry.

The ground-based survey was conducted on May 22nd and 23rd, 2023. The radar operated at 180–480 MHz, sending out 2 μs chirped pulses at 5 kHz pulse-repetition frequency. The transmit and receive were both polarized, and the system used switches to acquire 4 different polarizations in sequence. The pulse order was HH, HV, VV, VV, VH, where H is across track and V is along track (e.g. HV means transmit across track, receive along track) and the VV mode was repeated to increase signal-to-noise while mapping layers. The system used 5 transmit and 6 receive channels. The radar electronics were placed in a PolyPod SnowCamper (enclosed sled) and towed behind a pair of roped-together snowmobiles, with the antennas on a balloon towed behind the PolyPod. Survey length was limited by rough surface conditions resulting from unusually high precipitation and winds during the preceding week. The soft, rough snow also prevented driving in straight lines, since the snowmobiles were unable to pull the radar directly uphill. In total, about 11 line-km of ground-based data were collected, focusing on the flattest area of the ice cap with surface elevation above 1850 m. Data processing consisted of coherent integration, pulse compression, interpolation to constant trace spacing, SAR focusing using the range-doppler algorithm, and median filtering. Tracing and interpretation again used ImpDAR (Lilien and others, 2020); due to fewer ambiguities, only one picker traced these data.

In the airborne data, the independent tracers all picked a bed reflector on 68% of traces, finding the remaining 32% ambiguous. The tracers agreed on the ice thickness to within 30 m on 77% of the data they traced, 52% of the total survey. The combined radar surveys provide a dense map of the ice thickness (Fig. 4b). Ice thickness varies between 180 and 750 m over spatial scales of a few kilometers. Near the top of the ice cap, surface topography and ice flow do not appear to correlate with this variation in ice thickness, but rather form a typical dome shape. Toward the eastward edge of the ice cap, deep troughs are evident beneath several of the outlet glaciers. Overall, there is a weak correlation (r² = 0.22) between thickness and ice-flow velocity. Previously available gridded products based on inverse modeling (Millan and others, 2022; Zinck and Grinsted, 2022) differ from the measured ice thickness by at least 135 m (RMS). The Millan and others (2022) product and Zinck and Grinsted (2022) ‘combined’ product are relatively bias free (−4.7 m and 20 m mean difference compared to all radar measurements, respectively), suggesting that they capture the ice thickness on a broad scale but miss kilometer-scale variations. The remaining ice-thickness maps from Zinck and Grinsted (2022) show significant bias compared to measurements (−57 to −156 m mean difference).

Internal layers are not visible in any airborne survey. In the ground-based data, layers are clearly visible through most of the ice column in much of the surveyed area (Fig. 5). From kilometer 2 to 3 in profile a–a’ (Fig. 5a) and from kilometer 0.5 to 1.5 in profile b-b’ (Fig. 5b), internal stratigraphy is relatively flat and is visible from about 200 m (when the switches to the amplifiers on the radar opened) to about 500 m depth. At the site that we find most promising, layers are visible down to 550 m, at which depth they become obscured by off-nadir reflections.

Figure 5. Ground-based radar survey. (a) Profile trending northeast. Yellow line shows bed. Blue arrow shows candidate drill site. (b) Profile trending southeast. (c) Location of profiles in other panels with colors indicating ice thickness. Background lines show survey from which the measurements originate, as in Figure 4a. Profile S1 to S1’ is shown in Supplementary Figure 1.

Temperature

Temperature was assessed from a weather station on the ice cap and from measurement at the bottom of a shallow borehole drilled at the candidate deep core site. The weather station was installed in May 2021 at 91.645°W 79.862°N, 1804 m a.s.l., and has measured temperature, humidity, wind speed, wind direction, and the distance to the snow surface through May 31st, 2024. Data were recorded hourly and telemetered every 4 h. Wind speed and direction are likely biased due to rime during the winter. It is possible that the temperature measurements have been affected by rime as well, though this error cannot be easily quantified (Huwald and others, 2009). However, temperature measurements are generally less affected by rime than wind since they involve no moving parts, so similar to other work relying on remote weather station data we treat the temperature data as reliable (Ettema and others, 2010). We supplement this record with data from the weather stations at Eureka (WMO identifier 71917), which have been in operation since January 1^st, 1953. During this time, three different stations have carried this WMO identifier, and in 2016 the station was repositioned from its old site at 10 m a.s.l. near the main buildings to near the runway, 2.5 km away at 83 m a.s.l. (Environment Canada, 2024). However, since 2013, a second weather station has run in the area (WMO identifier 71613, located between the old and new stations at 10 m a.s.l.). Using the difference between each station carrying the 71 917 designation and the 71 612 station during their periods of overlap, we stitched together a consistent record from January 1, 1959 (when hourly records began) to May 31^st, 2024, matched to conditions at the original station site (at 10 m a.s.l.). This resulted in a downward shift of the recent temperature measurements by about 0.82 °C to match the lower elevation, which was likely needed due to strong inversions in the area in winter (Serreze and Barry, 2014). We assessed the correlation between temperatures at Eureka and those on Müller Ice Cap during May 2021–May 2024 and used the trends in the Eureka data to infer how temperatures on the ice cap are likely to have changed. As an additional check on the temperature, we measured the temperature at the bottom of the 6.62 m borehole using an EasyLog EL-USB-1-PRO temperature logger. The logger had 0.2 C precision and was left to equilibrate overnight.

To extend the observational record on top of the ice cap, we used Sentinel-1 imagery to determine the melt season. As in previous work (Wang and others, 2005; Sharp and others, 2011; Zinck, 2021), we used reduction in radar backscatter intensity as an indication of surface darkening due to the presence of meltwater. We used an already-orthorectified and radiometrically calibrated library of all available Sentinel 1 scenes (Mullissa and others, 2021), out of which we selected all scenes collected in the extra-wide swath mode, HH polarization. We binned the surface topography at the summit of the ice cap into three bands: 1600–1700 m, 1700–1800 m, and >1800 m, and averaged radar backscatter across all pixels by elevation band to reduce noise. Following previous work (Zinck, 2021), we used −5 dB as cutoff for melt days, though due to the abrupt nature of the observed changes in backscatter, the results are insensitive to cutoff values between −4 and −10 dB. Calculations were performed in Google Earth Engine (Gorelick and others, 2017).

The mean annual surface temperature from May 3rd, 2021 to May 31st, 2024, was −19.6 °C (calculated as the average of all 3-year periods within this window). Temperatures varied between about −35 and 0 °C seasonally, with minimum −47 °C and maximum 21 °C over the entire record (Fig. 6b). The weather station likely captures the start of four melt seasons, with the day when temperatures first reached above zero ranging from May 5^th to May 26^th. The temperature measured at the bottom of the 6.62 m borehole was −21.10 °C. The difference between this value and that from the weather station can be attributed to the residual seasonal cycle at the borehole bottom, which is explored with a model below.

Figure 6. Temperature and melt observations. (a) Temperature at Eureka and on Müller Ice Cap. Colored lines show trends pre- and post-1973 for the Eureka data. (b) Data from the ice-cap weather station. Labels indicate first day of each year above freezing for both the high and low. (c) Radar backscatter from Sentinel-1 averaged over different elevation bands. Labels indicate the day of each year that the backscatter of the >1800 m band first dipped below −5 dB (a commonly held threshold for melt).

The temperature at the weather station on top of the ice cap is highly correlated with that at Eureka (r² = 0.94), with a mean offset of −2.55 °C to the harmonized Eureka record. This offset varies seasonally, with the ice cap 7.3 °C colder than Eureka in summer and 2.8 °C warmer than Eureka in winter. This variation is consistent with observations at Eureka showing decreasing temperature with altitude most of the year and strong inversions during winter (Lesins and others, 2010). While winter temperatures at the weather station on top of Müller Ice Cap were warmer than Eureka, they were not as much warmer as would be expected in the free atmosphere, which warms by about 10 °C in the first 1200 m before cooling about 1 °C km⁻¹ above (Walden and others, 1996). This difference is likely due to radiative cooling of the ice cap surface causing an inversion there as well, similar to over other large ice masses and consistent with the pervasive wintertime inversions in the Canadian Arctic (e.g. Schwerdtfeger, 1970; Lesins and others, 2010; Tikhomirov and others, 2021). In summer, the offset implies a vertical lapse rate of 4.0 °C km⁻¹, lower than canonical values of 6.0 to 6.5 °C km⁻¹ but similar to lapse rates found over mountainous regions (Minder and others, 2010; Shen and others, 2016). In spring and fall, the offset between the ice cap and Eureka was intermediate to those in summer and winter (the ice cap was 1.4 °C cooler than Eureka in spring and 4.3 °C cooler in fall), consistent with some inversion-dominated and some typical lapse rate conditions. The data from Eureka indicate a temperature increase of 0.071 °C yr⁻¹ for the 50 years ending in 2022 (Fig. 6a). During the preceding 13 years, temperatures were decreasing at 0.069 °C yr⁻¹. This result is consistent with past analysis of temperature trends at Eureka, which found 0.088 ± 0.017 °C yr⁻¹ between 1972 and 2007 (Lesins and others, 2010). Overall, these trends suggest that the modern surface temperature at Eureka, and thus on top of Müller Ice Cap, is about 2.6 °C warmer than in 1959, when records began. However, satellite remote sensing indicates more rapid warming (0.32 ± 0.01 °C yr⁻¹) above 1400 m on Umingmat Nunaat during 2000–2015 (Mortimer and others, 2016), suggesting that this may underestimate warming at the top of the ice cap.

Radar backscatter from Sentinel 1 indicates that the surface of the ice cap above 1800 m darkened during 6 of the 8 summers from 2016 through 2023 (Fig. 6c). During these summers, this high-elevation area had reduced backscatter for an average of 18.2 days, suggesting melt took place in that period (compared to 21.2 and 26.0 days for 1700–1800 and 1600–1700 m, respectively). However, for the >1800 m area, the date when backscatter first decreased in 2022 and 2023 was over a month later than when the temperature first exceeded 0 °C as measured by our weather station. In 2021, the temperature exceeded 0 °C almost daily from May 16th to August 24th, though no reduction in backscatter was observed. Rather than corresponding to the date when the temperature exceeded freezing, reduced backscatter loosely corresponds (within 5 days) when the daily low was first above freezing. However, the final day with reduced backscatter did not appear related to the final day when the daily low was above freezing. Thus, while the radar backscatter may indicate surface water, it does not appear to be sensitive to small amounts of melt and is not tightly coupled to the temperatures exceeding freezing. We hypothesize that the reduction in radar backscatter corresponds to a quantity of meltwater that is unable to percolate into the firn, and thus remains near the surface causing a reduction in backscatter; further examination is left to future work.

Surface mass balance

We determined the SMB from shallow ice cores, a snow pit, and the weather station. The pit was dug to 67.5 cm on May 19^th, 2023, at 91.687°W 79.870°N, 1812 m a.s.l. (Fig. 1d), at which point a firm layer was impenetrable for shovels. The pit was continued with a hand saw for another 30 cm below this layer.

Two firn cores were drilled with a Kovacs hand drill over the course of May 21^st to May 25^th, 2023, at 91.687°W 79.870°N, 1812 m a.s.l. and 91.795°W 79.874°N,1863 m a.s.l., respectively (Fig. 1d). The first core was stopped at 1.60 m depth since many solid ice layers prevented further progress. The second core found easier drilling and was stopped at 6.62 m due to time limitations. Visual stratigraphy was recorded for both cores. Samples from each were taken for δ¹⁸O, and density was measured on the 1.6 m core; problems with the scale and inability to bring back frozen samples prevented density measurements on the second shallow core. δ¹⁸O was measured in 2.5-cm increments on the 1.6-m core and 5-cm increments on the 6.62-m core using a cavity ringdown spectrometer (Picarro model L2130-i Isotopic H₂0 couple with vaporization module A0211). Two microliters of water were injected into a vaporization chamber, and the vapor was then passed into an infrared absorbance cavity. Nine injections were taken from each sample with the first three excluded to remove any residual carryover from the previous sample. Vapor content of δ¹⁸O values were calculated relative to three standards, Vienna Standard Ocean Water 2 (VSMOW2), Greenland Ice Sheet Precipitation (GISP), and Standard Light Antarctic Precipitation 2 (SLAP2), which were measured at intervals of ten samples. The measured δ¹⁸O is a proxy for temperature, with lower values indicating colder temperature at the time that the ice accumulated (e.g. Dansgaard and others, 1969) and can thus be used to infer the time span of the firn core using the seasonal temperature cycle.

To obtain SMB from the 6.62-m record, we identified likely winter temperature minima in δ¹⁸O, which we used to estimate the time span of the core. To better constrain the length of the firn-core record, we also cross-correlated the δ¹⁸O profile with the temperature profile at Eureka in two ways. First, we did a direct cross-correlation of δ¹⁸O vs ice-equivalent depth to the temperature vs time, effectively assuming the δ¹⁸O represents the temporally averaged temperature. Second, we compared δ¹⁸O vs ice-equivalent depth to temperature vs accumulation, effectively assuming the δ¹⁸O represents the average temperature during a fixed amount of accumulation. For each of these comparisons, we tested various assumptions about the length of the firn-core record, ranging from 10 to 25 years in daily increments and allowing up to two months of offset between the Eureka record and that in the firn core. To obtain the net accumulation during the period spanned by the firn core, we assumed that the density profile followed the other core to 1.6 m, and then had constant density below that. Since the density from 1.6 to 6.62 m is highly uncertain, we assumed it follows a truncated normal distribution with mean 700 kg m⁻³, std dev. 100 kg m⁻³, and limits 350 and 910 kg m⁻³. This distribution was chosen to span the density of surface snow to that of ice, incorporating the large uncertainty in density into the inferred SMB.

The weather station has measured distance from a sonic sensor to snow every hour from May 2021, to May 2024, though some measurements were affected by rime. The sensor was raised annually from 2022 to 2024 to prevent burial and change in height was recorded to maintain a continuous record. Rime caused the range to snow to be recorded as 0 m, so these samples were easily discarded before calculating net SMB. The sensor has ± 1 cm accuracy. From this record, we calculated annual SMB over a two-year period from May 2021 to May 2023, assuming the snow had a constant density of 550 kg m⁻³, based on the density at the top of our snow pit.

Visual inspection of the firn core revealed 11 isolated 2–10 cm melt layers surrounded by dry firn (Fig. 7a). Dry firn between ice layers ranged from ~0.1 to ~2 m in length. We did not use these layers to estimate the time span of the core, since there is not necessarily a one-to-one correspondence between melt layers and summers; some summers may experience too little melt to form a layer or melt from multiple summers may percolate down and form a single layer. The irregular spacing of the melt layers further suggests that the melt layers are not annual.

Figure 7. δ¹⁸O in firn cores. (a) δ¹⁸O records from 6.62-m core (black) and 1.6-m core (maroon). Blue region shows depth of coarse winter snow. Red dashed line shows expected depth of time when weather station was installed, May 4^th, 2021. Blue dashed lines show solid ice layers in the 6.62-m core. Markers indicate likely winter minima, with color and shape indicating confidence that the minimum represents winter. (b) Comparison between δ¹⁸O (black) and temperature at Eureka (full-resolution, dark gray; averaged for comparison, light gray). 12 years of temperature data plotted, leading to maximum correlation with δ¹⁸O (r ² = 0.31). Y-axis for δ¹⁸O has been converted to ice-equivalent depth. (c) As in b, but Y-axis for temperature is accumulation at Eureka rather than time, leading to maximum correlation with a 22-year record (r ² = 0.38).

From the 67.5-cm snow pit, SMB was found to be 260 kg m⁻² yr⁻¹, representing only one year's worth of data. Average SMB from the weather station was found to be 270 kg m⁻² yr⁻¹, derived from just over two years of data. The δ¹⁸O data from the 6.62 m core have 11 high-confidence winter minima in (Fig. 7a), defined as local minima that are the lowest in a neighborhood with 45 cm full width. There are an additional 5 medium confidence winter minima that are the most negative in a 35-cm neighborhood and 5 low-confidence winter minima that are the most negative in 15-cm neighborhoods. We assign 95% confidence that the high confidence minima represent winters, 66% medium confidence, and 25% for the low confidence. The probabilities that the layers represent winter minima combine to form a Poisson binomial distribution for the total length of the record, indicating that there are 14.8 ± 2.3 years of data in the record. Combined with the density and its uncertainty, this suggests that the ~15-year average accumulation on top of the ice cap is 309 ± 46 kg m⁻² yr⁻¹.

Correlation with the temperature record at Eureka was maximized with a 4416-day (12.09-year) record when δ¹⁸O is assumed to represent a temporal mean, and with a 7983-day (21.86-year) record when it is assumed to represent the mean over a fixed amount of accumulation. The correlation is approximately bimodal for each comparison, peaking around 12 years and 22 years (Supplementary Figure 2). These would suggest accumulation rates of 373 ± 36 and 182 ± 18 kg m⁻² yr⁻¹, respectively. However, correlations are weak, (r² = 0.31 and 0.38, respectively), and there is little visible resemblance between the temperature and δ¹⁸O records (Fig. 7b–c), so we place greater confidence in the peak-based estimate of the record length. Visual inspection of the snowpack showed coarse snow (indicative of cold temperatures, thus interpreted as winter accumulation) from 47.5 to 60 cm, which aligns with the first high-confidence minimum in the δ¹⁸O curve (Fig. 7a) but does not match with the 2022–2023 winter in the Eureka temperature record at 12- or 22-year scaling. Based upon the total surface-height change at the ice-cap weather station, the May 4^th, 2021 surface should be at 1.06-m depth. 1.06 m lies between the second and third high-confidence minima, consistent with those picks, but also does not match with the Eureka temperature record at 12- or 22-year scaling, further supporting peak-counting to determine the SMB.

Potential ice-core site

We now consider what location on Müller Ice Cap contains the longest climatic record. With the limited thickness of the ice cap (mostly <750 m) and reasonably cold surface temperatures (mean annual temperature of −19.6 °C), we anticipate no basal melt, so thicker ice should lead to a longer record and better temporal resolution. We set a minimum of 500-m thick ice to be suitable for a core, since otherwise age resolution near the bed will likely be very low. If we require that ice-flow speeds be less than 5 m yr⁻¹, we can exclude the top portion of the deep troughs that underly the onset of the eastward-flowing outlets, and instead consider deep pockets of ice near the cap's summit. Given the number of days above freezing each year, high surface elevations are important to avoid surface melt, further indicating that the best site to drill would be near the top of the ice cap; we thus set a minimum surface elevation of 1800 m a.s.l. for a potential ice-core site. We also desire radar-detectable internal stratigraphy, since it can exclude the possibility of flow disturbances or significant degradation of the record due to surface melt; while such stratigraphy might be found anywhere in the ice cap, it can only be confirmed in the region we surveyed with ground-based radar. The highest point with thickness measurements is underlain by a bedrock peak and has relatively thin ice (91.859°W, 79.888°N, ice thickness ~435 m, Fig. 1d), so the best site is likely somewhere near this peak surface elevation but over a trough or valley that permits thicker ice.

The most promising site we can identify on Müller Ice Cap is at 91.928°W, 79.884°N (Fig. 1d). That location has ice about 606 m thick and surface elevation 1863 m a.s.l., thus satisfying considerations for finding thick ice with high surface elevation. The site lies slightly east of the ice divide but has slow ice flow (2.1 ± 0.4 m yr⁻¹). To within error, the ice has not changed thickness during the last 11 years. In addition, this location has clearly visible radar stratigraphy down to about 550 m (Fig. 5), as deep as anywhere we surveyed. This location is near top of the ice cap, about 20 m below the highest elevation that we were able to survey and 50 m below the highest point overall, suggesting that surface melt is likely to be close the cap-wide minimum.

Modeling

Having identified a promising core site, we turn to modeling temperature and age to determine whether a long, high-quality paleoclimate record is likely to exist there.

Modeling methods

Ice temperature

We modeled the ice temperature to ensure that the basal temperature at the potential site is below the pressure melting point. Temperature was simulated using the standard 1-D advection-diffusion equation in the vertical (e.g. Zwinger and others, 2007)

(1)$$\displaystyle{{\partial T} \over {\partial t}} = w\displaystyle{{\partial T} \over {\partial z}}-\kappa \displaystyle{{\partial ^2T} \over {\partial z^2}}.$$

where w is the vertical velocity (positive downward) and κ is the thermal diffusivity of ice (assumed to be a constant 32.19 m⁻² yr⁻¹). The vertical velocity was determined using a Lliboutry model (Lliboutry, 1979)

(2)$$w = \dot{b}\left({1-\displaystyle{{\,p + 2} \over {\,p + 1}}\left({1-\displaystyle{z \over H}} \right) + \displaystyle{1 \over {\,p + 1}}{\left({1-\displaystyle{z \over H}} \right)}^{\,p + 2}} \right), \;$$

Eqn (2)where $\dot{b}$ is the surface mass balance, H is the ice thickness, and p is a shape factor, here taken to be 3 following recent work where greater information about vertical strain was available (Chung and others, 2023). Eqn 2 assumes there is no basal melt (so that w goes to zero at the bed) and that the ice thickness is unchanging. Some of our simulations use variable ice thickness, so using Eqn 2 to describe the velocity amounts to assuming that the vertical velocity instantaneously adjusts to match the new ice-cap geometry. We solved Eqn 1 using discontinuous Galerkin finite elements with upwinding, implemented in the open-source Python package Firedrake (Rathgeber and others, 2016). The domain was discretized into 100 intervals (see Supplementary Figure 3 for a demonstration that this resolution is sufficient for numerical accuracy). The model was run for 40 000 years using a fourth order Runge-Kutta scheme with a 12-daytime step (Shu and Osher, 1988). Only the most recent 20 000 years of the model output were considered to ensure that there was no dependence on the poorly constrained initial conditions.

We ran the model with four hypothetical ice-thickness histories since there are few constraints on the long-term thickness history of the ice cap (Fig. 8c). One history used constant thickness, a null hypothesis. One used ice twice as thick as present (~1200 m) until 11.7 ka, then constant thinning to the modern thickness at 4 ka; this history effectively assumes that substantial melt did not occur during the HCO, so Müller Ice Cap has not been thinner than present during the last 20 ka. This sequence follows available evidence for nearby Agassiz Ice Cap (Vinther and others, 2009) and Barnes Ice Cap 1200 km to the south (Koerner and Fisher, 2002; Sneed and others, 2008) and is consistent with geological evidence of the retreat of the Innuitian Ice Sheet after the LGM (Dalton and others, 2020). The final two simulations assume that the ice was thinner during the HCO. These simulations used thick (~1200 m) ice until 11.7 ka, relatively rapid thinning to 200 or 406 m thickness at 8 ka, then thickening to modern conditions from 4 ka to 2 ka. These latter histories compromise between evidence from Greenland ice caps like Flade Isblink and Hans Tausen that completely disappeared during the HCO (Hammer and others, 2001; Lemark and Dahl-Jensen, 2010), and Canadian ice caps that survived essentially at modern thickness like Agassiz Ice Cap (Vinther and others, 2009) while remaining consistent with the geological evidence (Dalton and others, 2020). We simulated the temperature for each of these thickness histories using four different values for modern accumulation (see Boundary conditions below), which resulted in 16 simulations of the temperature profile through time.

Figure 8. Paleoclimate data and hypothetical ice thickness used as forcings for model simulations. (a) Surface temperature, stacked from weather-station and ice-core data (Vinther and others, 2009; Kindler and others, 2014). (b) SMB assuming different modern values: 160 kg m⁻² yr⁻¹ (Koerner, 1979), 255 kg m⁻² yr⁻¹ as suggested by our snow pit, 273 kg m⁻² yr⁻¹ as suggested by the weather station on the ice cap, or 373 kg m⁻² yr⁻¹ (Müller, 1962). (c) Hypothetical ice-thickness histories.

Firn temperature

We used a similar model to the one above, but with smaller domain and shorter timesteps, to determine how representative the temperature at the bottom of the 6.62-m borehole is of the annual mean. The model was again implemented in Firedrake (Rathgeber and others, 2016). As opposed to the full-thickness temperature, this problem is diffusion dominated so we solved it using standard second order Galerkin finite elements without upwinding. We used a 50-m domain discretized in 10-cm increments. Time stepping used daily increments with a backwards Euler scheme. Downward advection took place at the velocity prescribed by the top 50 m of the Lliboutry profile described above. The surface temperature was forced to match that from the weather station, repeated 100 times for a 200-year record so that there is no sensitivity to initial conditions. At the bottom, the temperature gradient was pinned to 0.025 °C m⁻¹, the gradient found using the model described above. The temperature cycle at 6.62 m stabilized after two years, but the 200-year simulation ensured steady state even at the bottom of the domain. The final temperature (on May 26th, the date when the temperature was logged), and the variability during the second half of the simulation, were extracted to evaluate the logged temperature compared to mean annual. Since the firn temperature is only weakly sensitive to the accumulation rate and ice-thickness history (the latter of which only affects the firn temperature through the temperature gradient at the bottom of the domain), we used only one simulation to evaluate the firn temperature.

Depth-age

We modeled the depth-age profile of the ice to determine the resolution with which climate data are likely to be preserved. The evolution of age, τ, follows a simple advection equation,

(3)$$\displaystyle{{\partial \tau } \over {\partial t}} = w\displaystyle{{\partial \tau } \over {\partial z}} + 1, \;$$

where the 1 is a source term to cause to increase with time (e.g. Zwinger and Moore, 2009). As with the temperature modeling, we assume that w follows Eqn 2. As with the temperature model, we allow the ice-thickness to vary to capture possible changes to the ice cap after the last ice age (since 20 ka).

We solved Eqn 3 using discontinuous Galerkin finite elements with upwinding, again implemented in Firedrake. The model domain was discretized into 9600 increments (see Supplementary Figure 4 for a demonstration that this resolution is sufficient for numerical accuracy). Model simulations were run for 250 kyr using the four ice-thickness histories described above and five SMB histories described in the next section, leading to a total of 20 simulations of the depth-age scale. These long simulations ensure that results are insensitive to initial conditions, since they allow us to only interpret ages in ice advected in from the surface during the simulations.

Boundary conditions

For boundary conditions, the models needed geothermal heat flux, SMB, and temperature. Heat flux measurements on ice-free areas of Umingmat Nunaat range from 50 to 54 mW m⁻² (Jones and others, 1990; Majorowicz and Embry, 1998), while a gridded dataset suggests 47 mW m⁻² best represents conditions beneath the ice cap (Colgan and others, 2022). Of these, we use the highest value 54 mW m⁻², to ensure that the model does not erroneously predict a frozen bed.

For the temperature boundary condition, we stitched together different data sources to create a gapless 50-kyr record describing conditions at the top of Müller Ice Cap. To match local ice-cap conditions, we offset the Eureka record by subtracting the 2.55 °C mean difference between the Eureka weather station and the top of the ice cap during the past two years. From 11.7 ka to January 1953, we use a synthesized record of northern hemisphere ice cores with elevation-change effects removed (Vinther and others, 2009). The Vinther and others (2009) record includes data from the Agassiz Ice Cap, relatively close to Umingmat Nunaat, so we expect it to do a good job representing the temperature trends in the area. The Vinther and others (2009) record extends to calendar year 1960, allowing us to offset it to match the already-corrected Eureka record during that time. Before 11.7 ka, we use the temperature record from the NorthGRIP ice core (Kindler and others, 2014), since that core balances being relatively close among the Greenland cores (at 1200 km away) while providing relatively reliable temperature estimates. As with the other datasets, we offset the NorthGRIP record to match the already-corrected Holocene record during the period of overlap (11.7–10 ka). After concatenation, these sources provide a jump-free temperature record corrected to ice-cap conditions, using NorthGRIP data from 50 to 11.7 ka, data from multiple northern hemisphere cores for 11.7 ka through 1952, the Eureka weather station from 1953 to May 2021, and a weather station atop the ice cap from May 2021 to May 2023 (Fig. 8a).

We followed a similar procedure for SMB, but, due to the strong dependence of accumulation and melt on topography, we did not use data from the Eureka weather station. We do have a record from the weather station on top of the ice cap indicating relative change in surface height since May 2021, which, paired with an estimate of the density, provides an estimate of the SMB over the last two years. However, given the short record, this leaves significant uncertainty in the modern climatological SMB of Müller Ice Cap (i.e. the past two years may not be representative of average conditions). While we also obtained estimates of SMB from the weather station, snow pit, and firn cores, and there are two previously existing estimates of SMB on top of the ice cap (Müller, 1962; Koerner, 1979), transport of mass by melting and refreezing, uncertain densities, and uncertain dating of the firn prevent any measurement from being definitive. We thus consider a range of possible modern accumulations, from 160 to 373 kg m⁻² yr⁻¹. For a given modern SMB, we scale a multi-core Holocene SMB record (Vinther and others, 2009) to minimize the misfit between its 1940–1960 average SMB and the modern value. Like for temperature, we use NorthGRIP data before the Holocene (Kindler and others, 2014), scaling the NorthGRIP record to match the multi-core record during the period of overlap (11.7–10 ka). We assume a constant value matching the oldest SMB from the NorthGRIP core prior to 100 ka. For a given modern SMB, we ended up with a record comprised of a constant value from 250 ka to 100 ka, NorthGRIP data from 100 ka to 11.7 ka, multi-core data from 11.7 ka to 1960, and we assume that the SMB has been constant at the modern value since (Fig. 8b). While this record leads to very small SMB during the last glacial period, low surface temperatures would have prevented melt during that period, thus allowing the ice cap to persist.

Model results

Temperature and melt

In all model simulations, modern basal temperatures are found to be −9 °C or colder (Fig. 9a). Basal temperature is weakly sensitive to which of the SMB and thickness histories are used as forcing, with modern values ranging from −15 to −10 °C. Past thinning leads to colder basal temperatures (since the bed was less insulated from cold surface temperatures during periods of thinner ice), while thicker ice in the past leads to warmer basal temperatures at present. Lower accumulation leads to warmer basal temperatures since cold ice is then advected downward more slowly. In all model simulations, basal temperatures remained well below freezing (below −8 °C) throughout the entire period of simulation (Fig. 9c). In all cases, the minimum temperature is ~ 150 m below the surface, where there is a ‘cold plug’ that has not yet been affected by warming surface temperatures or geothermal heat from below. This cold plug suggests that if surface temperature were to remain steady at modern values, englacial ice temperatures would rise to re-equilibrate over the course of centuries, and basal temperatures would eventually rise, too.

Figure 9. Temperature modeling. (a) Full-depth temperature modeling. (b) Firn-temperature modeling. Line color indicates date, as shown by color bar. The simulation ended on model date May 26, 2023, corresponding to the day when we measured the temperature at 6.62 m depth on the ice cap. On that date, the observed temperature was −21.1 °C and the simulated temperature was −21.21 °C. (c) Basal temperature in different simulations, with colors matching a; none approach the pressure melting point.

The temperature model indicates that at 6.62-m depth, there is a seasonal swing of about 4 °C, lagged by 128 days from the surface temperature cycle (Fig. 9b). The model found a temperature of −21.21 °C on May 26, which can be compared to the measured −21.1 ± 0.2 °C on that day. The excellent agreement between the model and measured borehole temperature indicates that there has not been measurable heating of the shallow firn because of heat transfer by melting and refreezing and suggests that the past two years of data observed at the weather station are representative of those 3 km away at the borehole during the last few years.

Age

The modeled depth-age scale varies significantly based upon the thickness and SMB histories (Fig. 10a). In the top 250 m, the depth-age scale is determined solely by the SMB, likely a result of all the thickness histories converging during the last 2 kyr (Fig. 8c). At intermediate depths (250–400 m), the shape of the depth-age curve is determined by a combination of the thickness and SMB histories, though the spread in modeled age at a given depth only reaches about 2 kyr. In the deepest third of the ice, ice gets older rapidly and spread between model runs reaches tens of thousands of years at a given depth. While not truly asymptotic, the rapid increase in age creates quasi-asymptotes at different depths. For a given thickness history, lower accumulation leads to a shallower age asymptote. For a given SMB, thinning at the beginning of the Holocene leads to a deeper asymptote, while thinning during the HCO leads to a shallower asymptote. With early Holocene thinning and low accumulation, ice ages reach 11.7 ka 149 m above the bedrock, while if the ice cap was reduced to 200 m thick during the HCO and accumulation is high, ice does not reach 11.7 ka until 8.1 m above bedrock. 56 m above bedrock, the deepest area where stratigraphy is visible in the radar, ages varied from 3.8 to 110 ka depending on thickness and accumulation history. However, it is likely that well-layered ice continues below this depth but was undetectable with the radar. The age at 20 m above bedrock varies from 5.3 to 222 ka between simulations.

Figure 10. Age modeling of the potential core site. (a) Modeled depth-age profiles. Line color indicates thickness history, as in Figure 9, while line style indicates SMB history. Note that the x-axis is truncated, as modeled ages in increase toward infinity at the bedrock due to the lack of basal melt. (b) Modeled height above bedrock of 11.7 ka ice vs the age resolution at 11.7 ka. Note the log scale for age resolution. Colors indicate thickness history, as in a, while marker shape indicates modern accumulation rate. Shaded gray areas show undesirable conditions (11.7 ka less than 20 m above bedrock or worse than 10 kyr m⁻¹ age resolution at 11.7 ka).

The age resolution at 11.7 ka also varies with SMB and thickness history, with most combinations leading to measurable age resolution and others leading to more than 10 kyr m⁻¹ at 11.7 ka (Fig. 10b). There is a pseudo-linear relationship between the height above bedrock of 11.7 ka ice and the resolution at that age for a given thickness history, suggesting a tradeoff between the time span recorded at the site and the resolution. However, this tradeoff does not hold between thickness histories; the relationship between resolution and height above bedrock of 11.7 ka varies substantially depending on thickness history for a given SMB.

Discussion

Temperature and melt

Modeling of firn temperature indicates that the measurements from our weather station and from the bottom of the 6.62 m borehole are consistent with a mean annual temperature of −19.6 °C. This consistency suggests that there is not major transport of heat through meltwater since this would lead to discrepancies between firn and surface temperatures. For example, meltwater runoff from the site of the firn core would cause anomalously low firn temperatures, since latent heat would be absorbed at that site and carried elsewhere with the meltwater. While there is some transport of heat vertically through the percolation of meltwater, the consistency of the weather station and firn temperatures suggests that this effect is minor.

Despite the reasonably cold mean annual temperature, there is substantial surface melt experienced by the ice cap. Even 60 years ago, before several degrees of warming in the region (Fig. 6a), the top of the ice cap experienced some melt in most summers (Müller, 1962). Any ice core drilled on the ice cap will thus have to contend with melt layers during the recent portion of the record; given the millennial-scale temperature trends (Fig. 8a), the ice cap likely has experienced surface melt ever since the surface thinned to its present elevation after the end of the last glacial. Since the period spanned by our shallow firn core is likely the warmest during the last 50 ka, the modern fraction of isolated melt layers places an approximate upper bound on the refrozen portion of the ice. While the melt redistributes mass in the shallow firn, mixing a few years of accumulation, it is unlikely to contaminate multi-year average paleoclimatic signals and still allows preservation of gas bubbles that could be analyzed in a core.

Given the modern conditions on Müller Ice Cap, our full-thickness temperature modeling indicates that there is unlikely to be any basal melt at present. This is consistent with measurements from other Canadian ice caps. At Agassiz Ice Cap, 350 km from Müller Ice Cap on Umingmak Nuna, basal temperatures were found to be colder than −16 °C in multiple boreholes (Fisher and others, 1983). At Prince of Wales Ice Cap, 280 km away on Umingmak Nuna, basal temperature below a 178 m core drilled at 1630 m a.s.l. was −19.6 °C (Kinnard and others, 2008). At Devon Ice Cap, 550 km south of Müller Ice Cap on Tallurutit (Devon Island), temperature at the bottom of a 300 m borehole drilled to bedrock was found to be −18.4 °C (Paterson, 1976). While Müller Ice Cap is much thicker than these other ice caps, and is thus likely warmer near the bed, the consistently cold temperatures found beneath these ice caps suggest Müller Ice Cap is frozen at the bed, too. Without substantial warming, basal melt would be unlikely during periods with a thinner ice cap as well. Available data suggest ~2 °C warmer temperatures during the HCO (Vinther and others, 2009), which when combined with a lapse rate of 6 °C km⁻¹ is insufficient to cause basal melt. It is also unlikely that the ice cap was sufficiently thick to experience basal melt during the Last Glacial Maximum, since our modeling indicates much colder basal temperatures during that period.

Likely age of the ice

Accuracy of the age scale for the Müller Ice Cap is limited by the considerable uncertainty in past SMB and ice thickness in the area. Estimates of accumulation range from 160 to 373 kg m⁻² yr⁻¹ averaging across various timespans. We estimated accumulation in three ways, finding values varying from 255 kg m⁻² yr⁻¹ for the last year to 309 kg m⁻² yr⁻¹ for the last ~15 years; attempts to use weather data from Eureka to improve this estimate were not successful. None of these estimates of accumulation give a definitive average accumulation over a long window, and all could be accurate for the time span that they measure. However, given the quantity of measurements now available, we think it likely that the climatological (30-year mean) SMB lies within the range spanned by these measurements, and that ~273 kg m⁻² yr⁻¹ is the most likely value. With these varying thickness histories and 273 kg m⁻² yr⁻¹ of accumulation, the age at 20 m above the bed is at least 17.8 ka except for the simulation where only 200 m of ice survived the HCO, in which case the age is just 6.1 ka. In that simulation, 11.7 ka is reached 11.9 m above the bed, with 12.6 kyr m⁻¹ resolution, closer to the bedrock and more compressed than desired, but potentially still sufficient to retain some useful climate information. 56 m above the bedrock, where stratigraphy is clearly visible, ages are expected to be younger even with 273 kg m⁻² yr⁻¹ accumulation (potentially only 3.8 ka); however, since this is thought to be a detection limit rather the true end of stratigraphy, a useful paleoclimate record extending beyond this age is very likely.

Though we can loosely determine a best estimate for the accumulation, there is little objective basis by which to distinguish among thickness histories. However, we can separate the effects of pre-Holocene thickness and thickness during the HCO on the depth-age scale. The model suggests that a factor of 2 variation in pre-Holocene ice thickness changes the height above bedrock and age resolution at 11.7 ka ice by 30% and 100% respectively (compare the constant thickness and Holocene thinning runs in Fig. 10). Though significant, this suggests that if the pre-Holocene ice cap were between 0 and 1000 m thicker than present, the climate record would still have desirable resolution and length. On the other hand, the model is extremely sensitive to thickness during the HCO. A 100% change in the thickness at the HCO results in ~200% change in the height above bedrock of 11.7 ka ice and up to a 1000% change in its age resolution (Fig. 10). Uncertainty in the depth-age scale is thus dominated by uncertainty in the Holocene ice thickness. The cores from Agassiz Ice Cap (Fisher and others, 1983) provide the best analogy for what can be expected at Müller Ice Cap, since they were drilled relatively close by (350 km east) and at similar elevation (up to 1730 m a.s.l.). Since the highest elevation core from Agassiz Ice Cap recovered a full-Holocene record in 176-m thick ice (Fisher and others, 1995), it seems unlikely that Müller Ice Cap was entirely absent at the HCO. 280 km southeast of Müller Ice Cap, Prince of Wales Ice Cap also dates to at least the Wisconsin glaciation (Kinnard and others, 2008), as does Devon Island Ice Cap 550 km to the south (Paterson and Waddington, 1984; Koerner and Fisher, 2002). Although Meighen Ice Cap, 150 km west of Müller Ice Cap, does not extend to the Wisconsin (Koerner, 1997), its much lower surface elevation (<270 m a.s.l.) suggests it is a poor analog. We thus consider it most likely that a full-Holocene record (>11.7 kyr), possibly with reduced quality as ages approach 11.7 ka, is preserved in Müller Ice Cap, but without direct age constraints on HCO thickness large uncertainty in the depth-age scale is unavoidable.

Conclusions

Using a variety of field and remote-sensing datasets, we identified a potential ice-core site on Müller Ice Cap at 91.795°W, 79.874°N. This site has surface elevation 1963 m a.s.l., 606-m thick ice, and visible radar stratigraphy through the top 550 m. It is flowing slowly without significant thinning at present. The mean-annual surface temperature is about −19.6 °C, as measured by a weather station and in the firn. Surface melt occurs annually at present, although recent warming suggests that deeper ice is likely to be less affected by melt. Accumulation in the area is between 160 and 373 kg m⁻² yr⁻¹, with a best estimate of ~273 kg m⁻² yr⁻¹ based upon measurements from a firn core, snow pit, and weather station. Modeling indicates that basal melt is very unlikely. If the ice cap was not more than 200 m thinner than present during the HCO, then ice older than 11.7 ka is likely to be preserved more than 20 m above bedrock and with better than 1000 yr m⁻¹ age resolution. These conditions suggest that a useful paleoclimate archive, likely spanning the Holocene, is preserved in Müller Ice Cap.