Study area

The study area comprises the northern part of the East Greenland Caledonides and the bedrock is primarily composed of crystalline basement except for a few places on Lambert Land and north of Nioghalvfjerdsfjorden, where it is overlain by Paleoprotozoic or Proterozoic sediments42. Topographic relief ranges between 0 and 500 m but locally mountains are up to ~1000 m high. The landscape is characterized by selective linear erosion; signs of glacial erosion are significant, particularly at lower elevations. During the LGM, the ice sheet advanced on to the continental shelf in Northeast Greenland43,44,45, but it has been contentious whether it reached the shelf edge ~250–350 km from the present ice margin or remained on the inner shelf (Fig. 1a)8. High-resolution multibeam swath bathymetry and shallow seismic data from the shelf offshore NEGIS show a number of glacial landforms including mega-scale glacial lineations, suggesting that the ice extended all the way the shelf edge in Northeast Greenland9,10,11,12. According to the existing 14C-based deglaciation chronology, the outer coast was deglaciated ~9.7–9.5 ka27. Following the deglaciation, the land was inundated with marine limits of 40–60 m above sea level46.

Cosmogenic exposure dating

In 2015 and 2016 we conducted fieldwork using helicopter and twin otter plane. We selected field sites using aerial and satellite imagery. Most samples were collected from boulders perched on ice scoured bedrock, except for three samples from Kap Amelié collected from boulders on drift and on Lambert Land where we collected two boulder samples on a moraine outside the Little Ice Age moraine (Fig. 1b). We aimed at sampling glacially abraded boulders on bedrock (Supplementary Fig. 1), with boulders >1 m in height and diameter to reduce the influence of snow cover on our resulting ages47. We measured shielding, and recorded the latitude and longitude and elevation using a handheld GPS with an uncertainty of <10 m. The boulder samples were collected using a rock saw. All samples were prepared using carrier “PHE1601” and were measured using the beryllium standard 07KNSTD48 at Aarhus AMS Centre (AARAMS). We used the CRONUS-Earth online calculator49, the Arctic production rate50, and time invariant scaling of Lal/Stone51,52 to calculate sample ages (Supplementary Table 1, Table 2). In addition, we used a rock density of 2.7 g cm−3 and made no correction for potential snow cover, and surface erosion. The study area has undergone glaciostatic uplift since the deglaciation of ~40 m at Blåsø and ~60 m at Hovgaard Ø46, and the sample elevation at the time of collection does not reflect its time-averaged sample elevation history. We calculated sample-specific elevation corrections25 and found that the uplift corrections are within 1-sigma of the AMS sample uncertainties similar to or lower than what have been calculated for West Greenland where the postglacial uplift was larger24,25. As the uplift history in Northeast Greenland is less well constrained compared to West Greenland46 and the vertical uncertainty of the GPS measurement is of the same order as the uplift correction, we present 10Be ages without correcting for glaciostatic uplift, similar to most other 10Be studies from Greenland19,20,22,23,53,54. We note that the lack of this correction does not significantly change our 10Be ages or our interpretations. Individual 10Be ages are presented with their 1-sigma analytical uncertainties, which include the uncertainty in the blank correction. When we compare our 10Be ages with other climate records we include the production rate uncertainty using the following formula:

$$\begin{array}{l}{\mathrm{Uncertainty}} = \\ \sqrt {\left( {1\sigma \,{\mathrm{std}}\,{\mathrm{error}}\,{\mathrm{of}}\,{\mathrm{mean}}} \right)^2 + \left({\,}^{10}{\mathrm{Be}}\,{\mathrm{age}} \times {\mathrm{production}}\,{\mathrm{rate}}\,{\mathrm{uncertainty}}\right)^2} \end{array}$$

$$\begin{array}{l}{\mathrm{Uncertainty}}\,{\mathrm{of}}^{10}{\mathrm{Be}}\,{\mathrm{ages}},{\mathrm{outer}}\,{\mathrm{coast}} = \\ \sqrt {(0.4)^2 + \left( {11.7 \ast 0.037} \right)^2} = 0.6\,{\mathrm{kyr}},\end{array}$$

where 0.4 = 1σ standard error of the mean (in kyr), 11.7 = mean 10Be age (in ka); 0.037 = the uncertainty associated with the Arctic production rate50 and “St” scaling51,52.

$$\begin{array}{l}{\mathrm{Uncertainty}}\,{\mathrm{of}}\,^{10}{\mathrm{Be}}\,{\mathrm{ages}},{\mathrm{present}}\,{\mathrm{ice}}\,{\mathrm{margin}} = \\ \sqrt {(0.2)^2 + \left( {9.3 \ast 0.037} \right)^2} = 0.4\,{\mathrm{kyr}},\end{array}$$

where 0.2 = 1σ standard error of the mean (in kyr); 9.3 = mean 10Be age (in ka), 0.037 = the uncertainty associated with the Arctic production rate50 and “St” scaling51,52.

We excluded two outliers based on the most general knowledge of the regions glacial history that are older than the LGM (GL1519, GL1545) and most likely contain 10Be inherited from a previous period of exposure47.

Radiocarbon dating

A number of shell fragments were collected on the surface of a moraine outside the Little Ice Age moraine on Lambert Land. The shell fragments were identified to species level, when possible. Only large pieces from a single specimen were used for dating. In the laboratory, shell fragments were cleaned and leached using HCl removing c. 25% of the outer shell. Around 10 mg of material was used for the 14C analysis; all contained enough carbon for AMS radiocarbon measurement. All radiocarbon ages have been calibrated into calendar years using IntCal1355 and a reservoir age of 550 years (∆R = 150)56 in Oxcal 4.357. The 14C ages are reported with 2-sigma uncertainty (Supplementary Table 3).

Temperature reconstruction

Greenland ice cores provide detailed records on the timing and magnitude of past mean-annual temperature change. However, GrIS mass loss occurs during the summer months58, and therefore summer (JJA) temperatures are more relevant than mean-annual temperatures when considering past margin positions. For the last 22 kyr Greenland ice core δ15N-based temperature reconstructions59 were merged with climate model simulations60,61,62 to generate Greenland-wide, seasonally resolving temperature reconstructions31. Supplementary Fig. 2a shows the reconstructed mean-annual (ANN), summer (JJA), and winter (DJF) temperatures at our study location (79oN, 20oW). Coupled ocean atmosphere GCM simulations are not available through MIS 3; therefore the same approach cannot be used to investigate summer temperatures during that time period. Instead, we rely on a multiple regression approach in which it is assumed that three key forcings dominate the Greenland temperature evolution: AMOC strength, greenhouse gas radiative forcing, and local orbital forcing. Summer temperature at the site (T JJA ) is then given by:

$$T_{\mathrm{JJA}} = a_1 \times {\mathrm{AMOC}} + a_2 \times 5.35\ln \frac{{p_{{\mathrm{CO}}_2}}}{280} + a_3 \times {\mathrm{\Phi}}_{\mathrm{JJA}}^{79^\circ {\mathrm{N}}},$$ (1)

where AMOC denotes the estimated strength of the overturning in Sv, \(p_{{\mathrm{CO}}_2}\) is the atmospheric CO 2 dry mixing ratio in ppm, \({\mathrm{\Phi }}_{{\mathrm{JJA}}}^{79^\circ {\mathrm{N}}}\)is the average solar insolation at 79°N north during the months June, July, and August, and a 1 through a 3 are linear scaling coefficients. The CO 2 mixing ratio is converted to radiative forcing using the approach by ref. 63, with a pre-industrial reference concentration of 280 ppm. All forcings are shown in Supplementary Figs. 2B-D. In reconstructing T JJA , we use a multi-ice core \(p_{{\mathrm{CO}}_2}\) compilation64, and insolation values calculated following ref. 65. The AMOC strength is the most uncertain of the three forcings and is reconstructed as follows. We start from the Greenland Summit66,67 δ18O record (average of GRIP and GISP2 δ18O records) corrected for mean ocean δ18O68, and convert it to site (mean annual) temperature using an effective isotope sensitivity of α = 0.29 ‰ K−159. Using the logic underlying Eq. 1, we remove the effect of CO 2 forcing (b 2 = 3.05 Km2 W−1) and insolation (b 3 = 0.0481 Km2 W−1 sensitivity to local summer insolation) from the GISP2 temperatures, where the stated (annual mean) sensitivities were obtained from the single-forcing deglacial GCM simulations of ref. 61, in which greenhouse gas and orbital forcings were applied separately. It is then assumed that remaining temperature variability is due solely to AMOC variability, which we estimate using a sensitivity (mean annual) of b 1 = 1.07 ± 0.25 K Sv−1, which optimizes the correlation to the reconstruction by ref. 31, as shown in Supplementary Fig. 2B. The NEGIS MIS 3 summer temperature anomaly is then estimated with Eq. 1 using coefficients a 1 = 0.238 ± 0.1 K Sv−1, a 2 = 3.76 ± 1.5 Km2 W−1, and a 3 = 0.137 ± 0.03 Km2 W−1; these coefficients are found using multiple regression analysis on the 0-22ka NEGIS JJA reconstruction; the coefficients are in good agreement with those found from the single-forcing coupled GCM experiments of ref. 61. The MIS 3 JJA reconstruction is shown in Supplementary Fig. 2e, together with an uncertainty envelope that was found by adding all stated uncertainties in quadrature. We want to emphasize the large uncertainty inherent in trying to reconstruct both the AMOC and NEGIS summer temperature based on regression techniques—the resulting values should be considered an order-of-magnitude estimate.

Data availability

The data that support the findings of this study are available in the supplementary information or it can be acquired from the corresponding author on request.