Introduction

Sea ice is a crucial component of the Arctic climate system, exerting a highly dynamic and highly seasonal influence over the energy balance, atmospheric and oceanic circulation patterns as well as biological and biogeochemical processes. In the present day, the rapid reduction in summer sea-ice area and the associated positive feedbacks are a source of tremendous concern. In particular, sea-ice loss plays a key role in Arctic amplification through its effects on surface albedo, surface heat fluxes and longwave radiation1,2. However, such feedback processes are certainly not unique to the present day. Previous interglacial periods, characterized by relative warmth and low global land-ice volume compared to the glacial periods that have dominated the Pleistocene, likely also saw corresponding reductions in global sea-ice area. Due to the limited quality and quantity of proxy records, it remains difficult to assess the magnitude or timing of any such changes in Earth’s history. Given persistent indications from both proxies and models that summer Arctic temperatures were considerably warmer than present day for at least brief periods during the Last Interglacial (LIG; ca. 130,000–115,000 years ago; also referred to as Marine Isotope Stage 5e, or MIS-5e) and the MIS-11c interglacial (ca. 425,000–395,000 years ago), it is very likely that large reductions in Arctic sea-ice extent also occurred during these periods. Noteworthy is also the fact that idealized modeling and some marine sediment core evidence suggests a role for sea-ice changes in helping to constrain the timing of glacial/interglacial transitions3,4, lending an element of criticality to understanding sea-ice changes in these periods.

The primary control mechanisms for the warmth and duration of the MIS-5e and MIS-11c interglacials include the Earth’s orbital parameters and changes in greenhouse gas (GHG) levels5,6. MIS-5e was characterized by a highly eccentric orbit, leading to an amplified seasonal cycle during its early phase (when precession was low) featuring a shortened but intensified boreal summer compared to preindustrial. MIS-11c, on the other hand, had much lower-amplitude insolation changes and maintained relatively warm conditions across almost two precessional cycles, something which appears to be unique among the last eight glacial-interglacial cycles7,8. Consistently high GHG levels helped maintain warm conditions throughout a long period of this interglacial8,9.

The specifics of the climate forcing mechanisms at work during the MIS-5e and MIS-11c interglacials notwithstanding, the overall dominant driver of Arctic climate is expected to be insolation levels. A frequently used metric for this external forcing is 21 June or summer solstice top-of-atmosphere insolation at 65°N, which is representative for high northern latitudes. However, this fails to consider the relation between seasonal intensity and length. According to Kepler’s second law, the square of Earth’s distance from the Sun is inversely proportional to its angular velocity, which implies that summer intensity and duration counterbalance one another10. Therefore, a more all-encompassing metric for insolation forcing is the integrated summer energy10,11 (ISE), a metric that integrates the insolation over all days exceeding a given insolation threshold throughout the year. This concept has thus far seen limited application to Arctic sea-ice levels, with one study12 using a similar approach to suggest some interglacials may have achieved sufficient whole-summer insolation levels for ice-free summers, and others not.

Evidence for a nearly sea ice-free summer in the Arctic is quite limited for the MIS-11c interglacial12, despite evidence that a higher sea level was achieved13 and that a much greater portion of the Greenland ice sheet melted9,14,15. In contrast, an earlier study16 argued for MIS-11 likely being the last interglacial period in which seasonally open ocean conditions predominated before sea ice became fully perennial in the central Arctic for the remainder of the Pleistocene. Thus, the limited sea-ice proxy data available for the Arctic dating back to MIS-11c is still highly uncertain and open to interpretation.

Arctic sea ice during MIS-5e is only slightly better understood. Multiple LIG simulations from the Climate Model Intercomparison Project 6 (CMIP6) suggest large reductions in summer Arctic sea ice. The multi-model mean minimum simulated sea-ice area reduces from 6.46 × 106 km2 in the preindustrial climate to 3.20 × 106 km2 during the LIG, with a few models indicating a near-complete loss of summer sea ice17. Another model-proxy data comparison study correlated the modeled Arctic summer temperature and ice-area changes across numerous models and compared those that best match the proxy temperature estimates (approximately 4.5 ± 1.7 K compared to preindustrial), suggesting an approximately three-quarters loss of summer Arctic sea ice in MIS-5e compared to preindustrial18. So far, proxy studies have presented conflicting evidence regarding changes in sea-ice cover in the central Arctic Ocean for MIS-5e. While sea-ice biomarker proxies indicated the presence of perennial sea ice in the central Arctic Ocean19, a recent study suggested a seasonally ice-free Arctic Ocean during the LIG based on the presence of sub-polar planktic foraminifera20. However, recent studies cast doubt on whether those records from the central Arctic Ocean can be assigned to MIS-5e21,22.

The present paper investigates Arctic sea-ice variability on long-term (multi-millennial) and short-term (interannual) time scales during the interglacials of MIS-5e and MIS-11c as simulated by a state-of-the-art coupled climate model. In particular, the following questions are addressed: If MIS-5e and/or MIS-11c did experience seasonally ice-free conditions, when and for how long might they have occurred? What is the main driver of multi-millennial sea-ice variability during the interglacials? How does interannual variability of the Arctic sea-ice cover depend on the mean state? In addressing these questions, we strive towards understanding the various sea-ice regimes that can emerge under different forcing conditions in the Arctic.

Results

The basis of this study is the analysis of a series of time-slice simulations that represent a chronological sequence of climate states over MIS-5e and MIS-11c using the coupled climate model CESM v1.2.223. Five-time slice simulations were performed for each interglacial in order to capture their long-term evolution, starting at 127 ka (kiloyears before present) and ending at 116 ka for MIS-5e, and starting at 423 ka and ending at 403 ka for MIS-11c. Details of the experimental design have been described elsewhere24. As the highly eccentric orbit of MIS-5e skews the relative lengths of seasons in comparison to today25,26, results are presented after calendar-adjustment of model output data (see section “Methods”).

Sea-ice area seasonality in the two interglacials and the role of insolation forcing

The climatological annual cycles in Northern Hemisphere sea-ice area (calculated as the sum of the sea-ice concentration times the area of a grid cell for all cells that contain some sea ice) for the various time slices of MIS-5e and MIS-11c are shown in Fig. 1. Particularly large variations in sea-ice area are apparent during the minimum ice-area months of August-September-October. With a total sea-ice area of much less than 1 million square kilometers, the September reaches ice-free conditions at 127 ka. At 124 ka, the sea-ice area partly recovers, but remains much lower than during preindustrial. All later time slices of MIS-5e show summer ice areas larger than preindustrial. For MIS-11c, minimum ice conditions of ~1.8 million square kilometers are reached at 413 ka, but only the latest time slice (403 ka) shows a larger-than-preindustrial Arctic sea-ice area.

Fig. 1: Comparison of Northern Hemisphere sea-ice area seasonal cycles.
figure 1

The monthly mean sea-ice areas (>45° N) for the final 100 years of each MIS-5e time-slice simulation are depicted in the left panel and MIS-11c on the right. The black curve in both panels is the climatology from the preindustrial simulation.

Full size image

Also very apparent is the rapid nature of change between time slices in MIS-5e as compared to those in MIS-11c. Despite having shorter time intervals between the simulations, the MIS-5e simulations rapidly progress from a high-amplitude seasonal cycle with virtually no September sea ice at 127 ka to a low-amplitude seasonal cycle with greater ice coverage than preindustrial in all months already by 121 ka. MIS-11c, on the other hand, begins with similar ice evolution to preindustrial at 423 ka, gradually evolves towards a higher-amplitude cycle with low September sea-ice minima at 413 ka, then shifts from a low summer sea-ice regime at 408 ka to having greater sea-ice area than preindustrial by 403 ka. There are, however, no drastic alterations of the seasonal cycle between the two sets of interglacial simulations, with the seasonal-cycle amplitude and rate of change between time slices standing out as the primary differences.

A closer examination of the time series of September months across the simulations shows the trend in sea-ice area between simulations even more clearly. Figure 2 is constructed by stitching together the time series of the final 100 years of each time-slice simulation, with spacing between the time series proportional to the length of time between the simulations. Therefore, while the time series themselves are not to scale with the times indicated on the abscissa, they proportionally show the evolution of September sea ice throughout each interglacial as well as the interannual variability in the September ice area. The extremely low sea-ice area levels previously discussed in 127 ka are clear to see, as well as the broad similarities between 124 ka and the 413–408 ka period in MIS-11c. Also evident is the rapid recovery of September sea-ice levels in MIS-5e and the drastically shorter duration of full-interglacial conditions compared to MIS-11c.

Fig. 2: Relationship between sea-ice area and integrated summer energy.
figure 2

Pseudo-time series of Northern Hemisphere sea-ice area for the month of September (left axis) and ISE at 65° N. Sea-ice area time series consist of the final 100 years of each quasi-equilibrated time-slice simulation, used here to represent each period designated on the abscissa. While time series are not proportional, the spacing between the time slices does correspond to the interval between forcing conditions, and thus the MIS-5e time series (dark blue, top axis) is shorter than the MIS-11c one (dark red, bottom axis). The dashed black line shows the mean September sea-ice area from the preindustrial control simulation.

Full size image

To investigate the forcing of the long-term sea-ice evolution, the high-latitude ISE throughout MIS-5e and MIS-11c is also plotted in Fig. 2. The values of ISE are defined as the integrated energy flux of all days in which the top-of-atmosphere insolation exceeds a given threshold at a given latitude (in this case, 325 W m-2 at 65° N; see section “Methods”)10. It is obvious that the sea-ice area in both interglacials is inversely correlated with the ISE, which undergoes a rapid transition from very high levels at 127 ka to much lower than present day beyond 121 ka (present-day value approximately 4.42 GJ/m2) and peaks at around 412 ka during MIS-11c. While the magnitude of the insolation decrease throughout MIS-5e is much higher, both interglacials see a sharp increase in September sea-ice area to above-preindustrial levels (dashed black line in Fig. 2) as ISE decreases towards and below preindustrial levels, suggesting the possibility of threshold behavior as simulated in a previous study27. While the ISE appears to be a good metric to explain the multi-millennial sea-ice evolution in the Arctic, the summer solstice (21 June) insolation is poorly correlated with the total sea-ice area during MIS-11c, and hence cannot explain the long-term sea-ice behavior (Supplementary Fig. 1).

Periods of particularly low ice levels are confined to the early stages of MIS-5e (127 and 124 ka) and the middle stages of MIS-11c (413 and 408 ka). The spatial contrasts in the September ice area for these periods are therefore comparatively illustrated in Fig. 3. At 127 ka, the simulated ice of at least 15% concentration has retreated to a small cluster along the coasts of northern Greenland and the northern reaches of the Canadian Arctic Archipelago. By 124 ka, however, the September mean sea-ice extent is already substantially recovered towards preindustrial levels (given by the bold red line), extending across the entire Arctic from Greenland and the Canadian Arctic Archipelago to the central Siberian coast.

Fig. 3: Spatial depiction of mean September sea-ice area for the two periods with the lowest ice per interglacial.
figure 3

Mean September sea-ice concentrations are given for the final 100 years of each indicated time slice simulation: 127 ka (top-left), 124 ka (top-right), 413 ka (bottom-left), and 408 ka (bottom-right). The bold black line indicates the 15% ice-concentration contour, below which conditions are regarded as open water. The red line shows the 15% ice-concentration contour for the preindustrial simulation. Colored icons in the MIS-5e simulation panels correspond to selected sediment core locations from other studies. Black filling indicates sites where proxy data have been interpreted to represent year-round sea-ice cover through the LIG19, gray filling indicating partial seasonal coverage, and white indicates fully open water conditions (see text). Green frames denote well-dated records, while orange frames indicate cores with poor chronology.

Full size image

A number of LIG high Arctic sediment records of various origins were compiled for comparison with Paleoclimate Model Intercomparison Project (PMIP) models17, a selection of which have also been plotted in the MIS-5e panels of Fig. 3. Among these are three poorly dated biomarker records (orange square frames in Fig. 3), which could provide insights into the probable presence or absence of sea ice on the basis of biological productivity of a number of species19. However, the lack of quality timing information strongly limits the utility of the cores to providing semi-quantitative estimates of consistently closed or (seasonally) open water conditions. These poorly dated cores are PS2757 (just poleward of the Laptev Sea), PS2200 (just north of Greenland) and PS51/038 (central Arctic), all of which provide biomarker records suggesting predominant perennial sea-ice conditions19. This result stands in stark contrast to a more recent study that analyzed microfossil records from across the central Arctic Ocean. Based on high abundances of the subpolar planktic foraminifer Turborotalita quinqueloba it was concluded that the Arctic Ocean became seasonally ice-free at some point during the LIG20. Our simulations are therefore consistent with this finding, with predominantly open-water conditions simulated at 127 ka and ice concentrations of 55-85% by 124 ka (and therefore a return to predominantly perennial ice conditions).

Among the better-dated core locations (indicated by green square frames in Fig. 3) are the PS2138 cores, located northeast of Svalbard19,28. These rely on dinocysts and the IP25 ice proxy, can be dated to MIS-5e, and are indicative of seasonally ice-free conditions during most of MIS-5e. This region is characterized by June-to-August-mean sea-ice concentrations of less than 5% in our 127 ka simulation and approximately 15–35% in our 124 ka simulation, providing ample opportunity for open-water phytoplankton species to be productive and thus in strong agreement with these proxies. Another well-dated biomarker record from core PS93/006-1 west of Svalbard is indicative of seasonally ice-free conditions during early MIS-5e, while proxy data from core PS92/039-2 north of Svalbard have been interpreted to represent near-ice margin conditions29. The southernmost record shown in Fig. 3 is M23455, which indicates year-round ice-free conditions in MIS-5e30. We note that further south, in the Iceland Sea, the model results also match proxy data from sediment cores M23352 and PS1247 that indicate nearly ice-free conditions all year round17. In sum, it appears plausible that the peak forcing conditions of early MIS-5e were sufficient to cause largely ice-free summers in the Arctic, but these were constrained to a short period of no more than a few thousand years.

Besides proxy data from high Arctic marine sediment cores, model studies and proxy data from outside the central Arctic Ocean may provide further inside into the state of Arctic sea ice during the LIG. A multi-model analysis found that LIG (127 ka) simulations that best match reconstructed temperatures over the Arctic region show a mix of ice-free and near-ice-free summers18, while studies with an isotope-enabled climate model found the signature of summer Arctic sea-ice loss in Greenland ice cores31. Furthermore, a recent compilation of multiple proxy data based on well-dated records from the Norwegian Sea indicates low-salinity conditions during the early LIG (128–124 ka), which has been attributed to enhanced melting of sea ice in the central Arctic Ocean32. Thereafter, surface salinity in the Norwegian Sea gradually increased throughout the LIG. Indeed, enhanced sea-ice melting reduces sea-surface salinity in the Arctic Ocean-Nordic Seas system. Our simulations reproduce these processes and show a similar temporal progression of surface salinity in the Nordic Seas as a consequence of a gradually increasing Arctic sea-ice cover during the LIG (Supplementary Fig. 2). As the sea-surface salinity has a direct impact on surface density in the Nordic Seas, these processes result in changes in convection, deep-water formation and hence the Atlantic Meridional Overturning Circulation (AMOC) as shown below.

The selected MIS-11c time slices demonstrate a less extreme evolution of sea ice, with 413 ka ice extent (15% boundary) slightly smaller than that at 124 ka, but with notably lower ice concentrations across the central Arctic basin. The 413 ka mean September ice area also broadly resembles the low ice pack of recent warmer summers since the year 2000, with ice fractions of 15–50% across much of the central Arctic and a primary ice center surviving near Greenland and the Canadian Arctic Archipelago. September sea ice in the 408 ka simulation already shows a recovery towards higher sea-ice concentrations and bears a very strong resemblance to the simulated 124 ka ice cover.

Interannual sea-ice variability and characteristics

Fields of September sea-ice concentration are subject to substantial interannual variability. For instance, during the 127 ka time slice individual years show larger total sea-ice area than some years during the 124 ka time slice (Fig. 2). To obtain further insight into sea-ice variability, standard deviations of the time series of September ice concentration from the last 100 simulation years per time slice were calculated (Fig. 4). Unsurprisingly, the 127 ka simulation stands out as having the least overall interannual variability, due to the very low ice concentrations across virtually the entire Arctic. In this case, the highest variability is actually concentrated along the northern coasts of Greenland and Ellesmere Island. Nevertheless, considerable variance in the central Arctic Ocean indicates the presence of individual years with substantial September ice cover.

Fig. 4: Interannual variability of September sea ice.
figure 4

Spatial comparison of the standard deviation of September sea-ice fraction across selected time-slice simulations. The black contour represents the 15% mean sea-ice boundary as in Fig. 3.

Full size image

A quite different pattern emerges for the other three time slices examined in Fig. 4: all of 124 ka, 413 ka, and 408 ka have local standard deviation minima in the Nares Strait and surrounding coastal regions, corresponding to the most stable area of summer ice in these simulations. Instead, maximum interannual variability occurs in the eastern hemisphere within the mean 15% sea-ice concentration contour. The pattern at 413 ka is slightly different, with a center of variability directly over the central Arctic, suggestive of increased melt sensitivity here near the MIS-11c peak in ISE. Overall, regions with moderate levels of mean ice concentration (approximately 20–60%) are associated with the greatest interannual variability.

This relation is also reflected in a comparison of the total sea-ice area means and standard deviations from entire time slices to each other (Fig. 5). Depicted is the mean September sea-ice area per time-slice simulation against the standard deviation of that 100-year time series. The relationship bears strong resemblance to a skewed bell curve with a long right tail. The greatest sea-ice standard deviations occur with a mean ice area of around 3.0 million km2 (the 124 ka and 408 ka time slices), tailing off quickly at higher and lower mean areas. Such a pattern can be both physically and mathematically explained: the cluster of simulations with sea-ice levels greater than the preindustrial simulation are characterized by sea ice across the entire central Arctic, extending to both the Siberian and Canadian coasts, and are thus geographically constrained. The very low-ice state at 127 ka simply does not have enough ice cover to produce a high standard deviation, which is partly due to the effects of self-reinforcing Arctic amplification helping to maintain consistently low ice levels year after year. The intermediate states, however, have moderate concentrations of sea ice across the Arctic, which are subject to rapid export, breakup, or melting, depending on the dominant weather patterns of a given summer. On the other hand, Arctic summers that skew cold, cloudy, or relatively stable could provide for less favorable melt conditions across a wide expanse of these vulnerable ice regions.

Fig. 5: Scatterplot of mean September sea-ice area against the standard deviation.
figure 5

Area values represent the means for the final 100 years of each time-slice simulation and the standard deviations are based on this same 100-year time series.

Full size image

Atlantic Ocean circulation

While we have attributed the differences in sea-ice evolution during the two interglacials to the insolation forcing (ISE), previous work has suggested that the LIG Arctic sea-ice cover might also be sensitive to changes in poleward ocean heat transport by the AMOC19. A stronger AMOC transporting more heat poleward has also been suggested to explain the long-lasting and anomalous (relative to insolation forcing) high-latitude warmth of the MIS-11c interglacial9. An expected consequence would therefore be an inversely correlated AMOC strength and Arctic sea-ice area. As depicted in Fig. 6, however, this does not appear to be the case in our simulations. Though the variations in the strength of the AMOC between time slices are relatively small and the interannual variability relatively large, the depicted trends would actually suggest a weak positive correlation between the AMOC and sea-ice area. During the minimum ice periods at 127 ka and 413 ka, the AMOC is at its weakest in our simulations, which can be attributed to low surface salinity in the Nordic Seas reducing surface density and hence deep-water formation (see above; Supplementary Fig. 2). The high-ice periods in the latter stages of MIS-5e are characterized by a generally strong but highly variable AMOC. While it is true that the AMOC in early MIS-11c (423 ka) appears marginally stronger than in the early stages of MIS-5e (127-124 ka), there is considerable overlap in their interannual variability. These observations come with the notable caveat that our simulations utilize a static Greenland ice sheet, and therefore the only meltwater or runoff contributed by Greenland in the model is from precipitation and seasonal snowpack. However, it appears that establishing a clear link between AMOC strength and sea-ice levels would require substantially more investigation.

Fig. 6: Relationship between the AMOC and September sea-ice area.
figure 6

Pseudo-time series comparing the September sea-ice area as in Fig. 2 for MIS-5e (left) and MIS-11c (right) simulations, but including a comparison to the mean annual strength of the AMOC.

Full size image

Our results are largely consistent with a recent model study suggesting that the AMOC plays a modest role for the ocean heat transport into the Arctic Ocean33. Instead, that study showed that ocean heat transport changes into the Arctic are dominated by changes in the Nordic Seas polar gyre strength. To assess a potential influence of the Nordic Seas circulation on the interglacial Arctic sea-ice cover we evaluated the polar gyre strength at 70° N for all time slices of MIS-5e (Supplementary Fig. 3) and MIS-11c (Supplementary Fig. 4). However, no correlation is found between the gyre strength and the Arctic sea-ice area such that we rule out a relevant role of the Nordic Seas gyre on the overall sea-ice evolution during MIS-5e and MIS-11c.

Surface temperature response

Disentangling external forcing and feedback from one another is notoriously difficult, but the robust seasonal temperature anomaly patterns in our simulations grant us some insight into the effects of the orbital forcing versus sea-ice feedbacks. The favorable alignment of precession and obliquity early in MIS-5e results in very high boreal summer insolation and thus, high ISE. Analysis of calendar-adjusted mean 2-meter air temperatures for all four seasons from our 127 ka simulation show the response to be a somewhat uniform increase in temperatures of 2–5 °C during the boreal summer season (June-July-August; bottom-left panel of Fig. 7). Over the central Arctic Ocean, anomalies are capped by the large uptake of latent heat from melting sea ice and the heat capacity of open waters and are generally in the 2–3 °C range. Far more extreme temperature anomalies emerge in the boreal fall, with +7–10 °C anomalies spanning the entire Arctic Ocean and much more modest warm anomalies dominate the adjacent land regions.

Fig. 7: Seasonal 2-meter air-temperature anomalies for early MIS-5e.
figure 7

Surface air temperatures given are seasonal averages from the 127 ka simulation: Winter (top-left), spring (top-right), summer (bottom-left), and fall (bottom-right). Absolute temperatures are depicted by the solid (≥0 °C) and dashed (<0 °C) black contours, while anomalies relative to the preindustrial control simulation are given by shading, with red depicting where interglacial temperatures were warmer than preindustrial and blue colder.

Full size image

Comparison of the fall temperature anomaly pattern with that from the 413 ka simulation (MIS-11c; Fig. 8) suggests that the extreme fall anomalies in both 127 ka and 413 ka are largely a result of sea-ice feedbacks in the Arctic. Whereas fall temperatures over Greenland and Eurasia at 127 ka are generally 1–3°C above preindustrial, they are at-or-below preindustrial levels over nearly all Northern Hemisphere land areas for the 413 ka fall season. We suggest that the extreme warmth over the Arctic Ocean largely results from heat released back to the atmosphere from an ocean that, in preindustrial and modern times, is already covered by early-winter sea ice for much of the season. This is consistent with previous analysis of the seasonal Arctic energy budget for MIS-5e, in which the largest Arctic temperature anomaly occurs in fall due to the anomalously-warmed Arctic Ocean returning stored summer heat and releasing water vapor to the atmosphere, which leads to a positive longwave radiative feedback34.

Fig. 8: Seasonal 2-meter air-temperature anomalies for mid-MIS-11c.
figure 8

As in Fig. 7, but for the 413 ka simulation.

Full size image

In both interglacial periods, the impacts of the less extensive sea-ice cover continue to manifest in substantial winter surface air temperature anomalies of generally 1–5°C over the Arctic, despite the absence of insolation during the polar night. This effect even lingers into the following spring, with Arctic air temperatures remaining slightly warmer than preindustrial, potentially portending an earlier onset to the sea-ice melt season. Simulations of the mid-Holocene suggest that this pattern is consistent across interglacial insolation maxima, with summer sea-ice losses proving persistent through fall and winter and driving warm anomalies in high latitudes35. These interglacial periods therefore exhibit obvious signs of a sea-ice loss driven Arctic amplification.

Discussion

In the interglacial time-slice simulations presented here are numerous noteworthy features that correspond to other simulations of past warm interglacial conditions or hypothesized features of ice-free Arctic summers in the near future. The large increases in sea-ice levels between some consecutive interglacial time slices (especially 124–121 ka) despite relatively modest ISE changes seem to indicate that Arctic climate may be susceptible to abrupt changes at insolation forcing thresholds. It is unclear if a similar sea-ice collapse threshold exists early in interglacials36, as our time-slice simulations do not extend far enough back into early MIS-5e and MIS-11c to avoid effects from the presence of remnant ice sheets at these times.

Seeking to understand the response of the Arctic sea ice to warm interglacial conditions is of course motivated by comparisons to the present and near-future, in which sea ice is projected to rapidly decrease. MIS-5e in particular has been suggested to be a key case study for ice-free summers in the Arctic, with apparent Atlantification of central Arctic waters paralleling trends already observed in the present day20. Indeed, our 124 ka and 408 ka summer sea-ice simulations already bear a strong resemblance to recent record-low sea-ice summers such as 2007 and 2012, suggesting a regime shift is already in progress. The 127 ka simulation provides a potential picture of the future state of the Arctic sea ice after this regime shift is complete: the disappearance of multi-year sea ice, a transition towards a high-amplitude seasonal cycle of ice cover, and very low levels of late-summer sea ice with low interannual variability for virtually the entire high Arctic.

Our chosen modeling methodology employs time-slice simulations with prescribed present-day topography, bathymetry, and land cover. Freshwater input from melting land ice was neglected, as were local geographic features like the Eemian Sea and Karelian seaway37. A previous study38 found that alteration of Greenland’s topography alone in MIS-5e simulations would lead to additional warming across most of the high latitudes, although primarily in winter and as a secondary effect to insolation forcing. Others have highlighted the role of vegetation feedbacks in reproducing annual warming at high northern latitudes in simulations of the LIG39. It is therefore likely that the robust Arctic warmth in our simulations would be even more enhanced in simulations with an interactive ice sheet and vegetation. On the other hand, substantial meltwater input could produce a counteracting cooling effect in response to a weakened AMOC and reduced ocean heat transport40.

Indeed, based on multiple proxy records from the Nordic Seas and North Atlantic it has been argued that meltwater, particularly from the Greenland ice sheet, was an important forcing factor for LIG ocean surface temperatures, likely leading to a regional delay of peak interglacial warming41,42. However, temperature evolution during the LIG was spatially inhomogeneous and asynchronous, with peak warming intervals partly antiphased between the Nordic Seas and North Atlantic43. A new study suggests an open-ocean environment in the eastern Fram Strait at 127 ka despite relatively cool surface conditions44. Subsequently, the influence of Atlantic Water on local sea-surface temperature gradually increased towards the late MIS-5e. To what extent ice-sheet melt has impacted the central Arctic sea-ice cover remains elusive at present. However, given the dominant impact of ISE on Arctic sea ice, as demonstrated in our study, we argue that maximum insolation forcing around 127 ka induced seasonally ice-free conditions in the central Arctic, despite a meltwater-driven delay of the thermal optimum in the Nordic Seas region that is not explicitly considered in our experiments.

In a larger context, Arctic sea-ice loss and the atmospheric circulation changes that it may induce can impact the degree to which the Greenland ice sheet experiences surface melt. Today, the increased presence of open water on the Greenland coast enables increased oceanic heat flux and increasing melt area, particularly over western portions of the ice sheet in August and September45. These satellite-based observations of the modern era are further supported by analysis of an ensemble of CMIP5 models, which indicate that although the effect is mostly localized and minor compared to overall atmospheric warming trends, loss of sea ice in the Baffin Bay and the Greenland Sea do indeed result in greater melt extent46. Idealized model experiments also indicate the potential for anomalous mid-tropospheric ridging to emerge over the Baffin Bay and western Greenland with sea-ice loss, potentially raising the frequency of blocking episodes over Greenland and enhancing the onshore flow of heat in specific regions47. The present balance of evidence would therefore indicate that Arctic sea-ice loss likely contributed to enhanced ice-sheet melting in the MIS-5e and MIS-11c interglacials and will do so again in the near future.

The time-slice approach applied in our study has obvious disadvantages compared to transient model simulations. For example, the timing of changes and their abruptness can only be determined vaguely. Moreover, this approach may overlook potential oceanic changes that could be important for the long-term evolution of the sea-ice cover over the course of the interglacials. However, due to their high computational demands, long-term transient simulations covering both the LIG and MIS-11c have so far only been carried out with coarse-resolution or intermediate-complexity models27. With increasing computing power, we expect transient simulations covering both the LIG and MIS-11c using more comprehensive models in the future.

Conclusions

Here we have sought to evaluate the relationship between Arctic sea ice, its interannual variability, and the ISE in two of the most critical interglacials of the Pleistocene, MIS-5e and MIS-11c. This has been done via a comparative time-slice modeling approach, although comparison with some sea-ice proxy reconstructions did provide indications that our MIS-5e simulations show what could be a realistic scenario for the summer sea-ice state during the LIG. A similarly direct evaluation of the model performance during MIS-11c is not possible due to the lack of suitable proxy data.

The two interglacials were subject to very different orbital forcing conditions, resulting in different patterns of sea-ice retreat that were strongly related to the ISE. Our model results suggest that the early stages of MIS-5e experienced extreme high Arctic warmth resulting from very high ISE and a consequent near-collapse of summer Arctic sea ice around 127 ka. However, a combination of (1) the rapid recovery of sea-ice levels as ISE dropped by the 124 ka simulation and (2) conflicting indications regarding seasonally ice-free conditions from the proxy data would suggest that any phase of ice-free summers during MIS-5e must have been confined to a brief period of a few thousand years or less. Such a period would have been characterized by the survival of virtually no multi-year ice, a high-amplitude seasonal cycle of sea-ice cover, and relatively low interannual variability in summer sea ice.

MIS-11c, on the other hand, did not achieve the same minimum summer ice area, but may have had 15,000 years or more (ca. 423-408 ka) with summer sea-ice levels below preindustrial levels, and without a highly pronounced insolation maximum. The evolution of sea-ice cover in MIS-11c is also suggestive of a high degree of sensitivity to the ISE rather than the peak summer insolation, which has traditionally been used attempting to explain sea-ice changes in previous studies19,29,30. The large increase of September sea-ice area from 124 to 121 ka and from 408 to 403 ka as ISE decreases below present-day values is also suggestive of a threshold response in summer sea-ice melt. In our simulations, the evolution of ISE would seem to explain virtually all of the sea-ice behavior, as ice-cover conditions are very similar between simulations from different interglacials with similar levels of ISE. Any additional contributions to Arctic sea-ice loss from ocean circulation changes appear to have been minor.

Finally, our simulations indicate a clear relationship between Arctic sea-ice area and its variability. Interannual variability in summer sea-ice levels maximizes for intermediate levels of ice cover, with very high and very low ice-cover states being geographically constrained. Intermediate levels of ice cover, however, are characterized by the presence of fractional ice cover across the central Arctic, where it is subject to export, mechanical breakup, and weather anomalies. On long-time scales, these more volatile summer sea-ice regimes are likely to signify transitions between relatively stable high-ice summers and low-ice summers (or vice versa), a transition that appears to be underway in the modern Arctic.

Methods

Model and experimental design

We use the Community Earth System Model CESM v1.2.223. The fully coupled system comprises model components for the atmosphere (CAM5), the ocean (POP2), the land surface (CLM4), and the sea ice (CICE4). POP2 and CICE4 share the same horizontal grid with a nominal resolution of 1°, where the North Pole is placed in the middle of Greenland. A horizontal grid of 1.9° × 2.5° is used for CAM5 and CLM4. CAM5 runs with 30 vertical layers and POP2 with 60 depth levels.

Our experiments for past interglacials aim to reproduce the climate conditions of the entire MIS-11c and MIS-5e periods. However, transient CESM simulations of 20,000 years or more are impractical with the computational resources currently available. A time-slice method has therefore been used to develop climate states across MIS-11c and MIS-5e. The time-slice methodology can fundamentally be described as using shorter simulations with constant forcing to provide windows of insight into much longer climatic periods. The time slices chosen for our experiments are based upon the identification of precession minima and maxima. Additional intermediate times are selected to fill large time gaps. Using Laskar et al.48 orbital parameters and insolation, five key points in the orbital forcing trajectory within MIS-11c were identified. These simulations and their respective justifications are as follows: 423 ka (intermediate point in orbital forcing, near the beginning of warm interglacial conditions), 418 ka (precession index maximum), 413 ka (intermediate point), 408 ka (precession index minimum), and 403 ka (near time of estimated sea-level highstand during MIS-11). Simulations for MIS-5e similarly cover the entire span of interglacial conditions. However, due to the substantially shorter duration of interglacial warmth compared to MIS-11c, time steps between the various time slices are shorter. The representative time slices are as follows: 127 ka (precession index minimum), 124 ka (intermediate point), 121 ka (intermediate point in precession cycle), 118.5 ka (additional intermediate point), and 116 ka (precession index maximum, near beginning of glacial inception in the Northern Hemisphere).

The simulations for each time slice are branched from a preindustrial control run, as described previously24, and integrated for 1000 years with constant forcing, enabling statistical equilibration of all climate system components except for the deepest ocean. Forcing for each time slice comprises orbital parameters48 and atmospheric GHG levels49,50,51. Following from the experimental setup of PMIP4, a nominal adjustment of +23 ppb was made for methane52. Due to its short atmospheric lifetime and having a disproportionate number of sources in the Northern Hemisphere, atmospheric methane concentration is subject to a hemispheric gradient. Therefore, based on a 46 ppb difference observed between Greenland and Antarctic ice cores for the early Holocene, the +23 ppb adjustment represents a more appropriate global average value than the raw Antarctic number52. Supplementary Table 1 summarizes the orbital and GHG forcing values used for each time slice.

Perhaps the most limiting assumption employed in our simulations is the use of present-day topography, sea level, and ice-sheet configurations. Meltwater from the Greenland ice sheet is also neglected. However, this model configuration does utilize the sea-ice model CICE4, which contains an explicit melt-pond scheme. In multi-model and proxy comparison studies, this has been shown to produce much greater loss of sea ice than for models that lack this feature18,53 and is arguably more realistic due to the degree of summer sea-ice melt that is driven by melt ponds54.

All presented results are climatologies or time series produced from the final 100 years of the time-slice simulations. Given anomalies are presented relative to the preindustrial control simulation24.

Testing the sensitivity to initialization

In order to demonstrate that the climate state that develops during each time-slice simulation is independent of the specific initial conditions (within reasonable limits), a sensitivity test was conducted. Two versions of the 418 ka simulation were run: one initialized using preindustrial conditions and one initialized using the final conditions from the 423 ka simulation. A comparison of the equilibrated periods (last 100 years) of each run has shown that no statistically significant differences exist between them. The preindustrial control initialization was therefore used for each of the time slice simulations, as this enabled running numerous time slices in parallel to each other.

Calendar-adjustment of MIS-5e data

The calendar problem for comparing monthly or seasonal climatologies from paleo-datasets or model simulations has long been recognized25,26. An angular calendar is a more appropriate choice for comparing climatologies across eras in which orbital arrangements substantially differed. We therefore utilize for this purpose the PaleoCalAdjust Fortran program26. This program recalculates the new start and end dates of all months in today’s calendar based on the orbital conditions of the past. The monthly data from the simulation is by default naively timestamped using an invariable modern 365-day calendar. It therefore must be deconstructed into quasi-daily values, interpolated, and re-aggregated to reflect the appropriate date range. All MIS-5e simulations are therefore subjected to this treatment to remove prominent calendar effects that would otherwise contaminate seasonal averages26,55.

Integrated summer energy

For the modern climatological annual cycle and high latitudes (65°N) the mean surface temperature is near freezing point (0 °C) when insolation intensity is between 300 and 350 W/m2 (see Fig. 2c in ref. 10). Therefore, we have taken 325 W/m2 as a threshold value. However, taking values of 300 or 350 W/m2 instead has little influence on the ISE curves (cf. Fig. S1 in ref. 10). For example, the MIS-11c ISE maximum is at 412 ka for the 325 W/m2 threshold (Fig. 2), at 413 ka for a 300 W/m2 threshold, and at 412 ka for a 350 W/m2 threshold (see Table S1 in ref. 10).