Northern Hemisphere sea ice variability in a transient CGCM simulation of the past 2.6 Ma

Mandal, Gagan; An, Soon-Il; Park, Jae-Heung; Yun, Kyung-Sook; Liu, Chao; Paik, Seungmok

doi:10.1038/s41467-024-55327-2

Download PDF

Article
Open access
Published: 02 January 2025

Northern Hemisphere sea ice variability in a transient CGCM simulation of the past 2.6 Ma

Nature Communications volume 16, Article number: 39 (2025) Cite this article

2789 Accesses
9 Altmetric
Metrics details

Subjects

Abstract

The recent sea ice changes in the Northern Hemisphere (NH), necessitate elucidating the sea ice variability over the past 2.6 million years (Ma), when the Earth’s glacial cycles transitioned from ∼41 to ∼100 kyr periodicity, following the Mid-Pleistocene Transition (MPT) period (0.7–1.2 Ma). Here, we analyze a coupled general circulation model (CGCM) simulation to understand how the NH sea ice responds to changes in the transient orbital, greenhouse gas (GHG), and ice-sheet forcings. We find that the Earth’s axial tilt (obliquity) and axial wobble (precession) strongly influence the variability in high-latitude (>70° N) perennial sea ice and mid-latitude (35° N–70° N) seasonal sea ice, respectively, by modifying the net surface shortwave radiation. Meanwhile, the GHG forcing affects the glacial–interglacial sea ice predominantly in the Labrador Sea, Irminger–Iceland basin sector, and Central North Pacific regions during the MPT and post-MPT (0.0–0.7 Ma) periods by modulating the downwelling longwave radiation. Additionally, we confirm that variability with longer periodicity (∼100 kyr) from GHG and ice-sheet forcings is most pronounced in NH sea ice during the MPT and post-MPT periods.

Antarctic sea ice multidecadal variability triggered by Southern Annular Mode and deep convection

Article Open access 08 November 2024

Synergistic atmosphere-ocean-ice influences have driven the 2023 all-time Antarctic sea-ice record low

Article Open access 03 August 2024

Simulated Antarctic sea ice expansion reconciles climate model with observation

Article Open access 07 January 2025

Introduction

Sea ice is sensitive to climate variations and is critical for understanding past and future climate changes. Sea ice amplifies climate change in high–latitude regions through its high albedo and minimal thermal conductivity, together with other processes (e.g., temperature feedbacks, changes in water vapor, cloud cover, and atmospheric and ocean heat transport)^1,2 that alter the radiative budget in the Arctic. Additionally, sea ice dynamics influence the mid-latitude jet stream³, oceanic salt and freshwater fluxes⁴, global intermediate/deep water circulation⁵, and North Atlantic Deep Water formation⁶. Over the past several decades, the Northern Hemisphere (NH) sea ice has been significantly reduced, in response to the increasing greenhouse gas (GHG) concentrations and global warming⁷. However, the responses of the NH sea ice to internal processes and external climate forcings in the past remain largely unknown. To comprehend historical and recent changes in NH sea ice^8,9,10,11, understanding the responses of sea ice to external climate forcings and their feedbacks during the Quaternary (past 2.6 million years (Ma)) is crucial.

The Quaternary period is characterized by glacial–interglacial cycles. The Earth’s orbital parameters, namely (1) orbital eccentricity (∼100 thousand years (kyr)), (2) axial tilt (obliquity; ∼41 kyr), and (3) axial wobble or longitude of perihelion (precession; ∼23/19 kyr), significantly vary the latitudinal and seasonal distribution of incoming solar radiation¹², which is the primary external forcing of the glacial–interglacial cycles^13,14,15. Additionally, GHG concentration changes and variations of the extent and topography of the NH ice-sheets, together with their associated climate feedback mechanisms, amplify changes initiated by orbital variations¹⁶ especially by altering the global mean surface temperature (GMST)¹⁷ and influencing the NH mid-latitude westerlies¹⁸. The NH ice-sheets experienced cyclic waxing and waning, indicating glacial–interglacial cycles in response to GHG concentration changes¹⁹. These glacial–interglacial cycles transitioned from a relatively low amplitude and shorter periodicity (\(\sim\)41 kyr cycle) during the Early Quaternary (1.2–2.6 Ma; pre-Mid-Pleistocene Transition (MPT)) to a high amplitude and longer periodicity (∼100 kyr cycle) in the Late Quaternary (0.0–0.7 Ma; post-MPT)^19,20,21,22, following the MPT (0.7–1.2 Ma)^{23,24,25,26,27}. Therefore, in this study, we divided the Quaternary into three distinct periods—post-MPT, MPT, and pre-MPT—to comprehend the evolution of the NH sea ice.

Paleoclimate modeling has shown that incoming solar radiation has driven the NH sea ice variations during interglacials over the past 0.8 Ma¹⁷. Additionally, various proxy records and modeling studies that are primarily confined to the post-MPT period identified ∼23/19 kyr (precession) and ∼41 kyr (obliquity) signals in the NH sea ice variability, with some emphasizing the role of precession^28,29,30, and others highlighting the role of obliquity^17,31,32,33. Recently, a transient paleoclimate simulation covering the past 3 Ma presented a general overview of the Northern and Southern Hemisphere sea ice variability, revealing precession, obliquity, and CO₂/ice-sheet (∼60–100 kyr) signals and suggesting the combined effect of the GHG and ice-sheet forcings on sea ice variability³⁴. Despite these findings, uncertainties remain regarding the spatiotemporal variability and mechanisms underlying the response of the NH sea ice to the orbital, GHG, and ice-sheet forcings. Furthermore, whether precession or obliquity strongly influences the sea ice remains unclear.

Here, we analyze the 3 Ma experiment and explain the spatiotemporal variability and mechanisms of NH sea ice change during the Quaternary. The 3 Ma experiment is a coupled atmosphere–ocean–sea ice simulation performed by Yun, et al.³⁴ employing the Community Earth System Model version 1.2 (CESM1.2), forced by transient prescribed orbital³⁵, GHG¹⁹ (including CO₂, CH_4, and N₂O concentrations), and NH ice-sheet forcings¹⁹. Although the GHG concentrations and ice-sheets are not external forcings in the natural system, they are treated as such in the 3 Ma experiment. Therefore, in this study, we describe them from an external forcing perspective. The 3 Ma experiment provides an unprecedented view of the global Quaternary climate variability, as it is able to reproduce consistent glacial–interglacial variability³⁶ and key paleo–proxy–based estimates^34,36. Nevertheless, the 3 Ma experiment design has a few caveats, such as the absence of prescribed ice–sheet melted freshwater forcing, a constant Last Glacial Maximum (∼21 ka) bathymetry, and an acceleration technique with a factor of five applied to the external (orbital, GHG, and ice-sheet) forcings.

Results

Evolution of Northern Hemisphere sea ice during the Quaternary

The prescribed orbital forcing³⁵ (Fig. 1a–c) influences the intensity of summer (June–July–August; JJA) incoming solar radiation at 65° N³⁷, which is the latitude at which the NH ice-sheets most sensitive to orbital forcing^15,38. The incoming solar radiation at 65° N consistently varies throughout the Quaternary period (Fig. 1d) with predominant ∼23/19 kyr (precession) and ∼41 kyr (obliquity) signals (Supplementary Fig. 1). The Quaternary climate shows gradual cooling towards the post-MPT; however the trend of JJA incoming solar radiation at 65° N does not change considerably during the past 2.6 Ma (0.0001 \({{{{\rm{Wm}}}}}^{-2}\)/kyr) (Fig. 1d). Meanwhile, the prescribed CO₂ concentration¹⁹ (GHG forcing) (Fig. 1e) (Supplementary Fig. 2a–c) and prescribed global ice volume in sea level equivalent¹⁹ (ice-sheet forcing) (Fig. 1f) (Supplementary Fig. 2j–l) show a dominant ∼41 kyr periodicity during the pre-MPT period, ∼60–100 kyr (CO₂/ice-sheet signal) periodicity during the MPT period, and a longer ∼100 kyr periodicity during the post-MPT period, consistent with the deep sea proxy record estimates²¹ (Supplementary Fig. 3).

**Fig. 1: 3 Ma experiment forcings and simulated climate using the CESM1.2 model.**

We investigated the NH (35° N–90° N) sea ice extent (ocean surface covered with more than 15% sea ice concentration) during the post-MPT (0.0–0.7 Ma), MPT (0.7–1.2 Ma), and pre-MPT (1.2–2.6 Ma) periods. Evidently, the simulated sea ice extent increases throughout the MPT and post-MPT periods (Fig. 1g) and expands equatorward, particularly in the West Pacific and Atlantic Ocean sectors (Fig. 2a, b). Here, we define > 70° N as high-latitude region and 35° N–70° N as mid-latitude region. Furthermore, perennial sea ice, i.e., sea ice that survives during summer and persists for over a year, is defined by the region north of the 90% contour during the MPT and post-MPT (purple line in Fig. 2a, b) and also corresponds to the region north of the 80% contour (purple line in Fig. 2c) during pre-MPT (cf. Supplementary Fig. 4). Evidently, the high-latitude region is primarily covered with perennial sea ice (north of the purple line in Fig. 2) throughout the Quaternary, whereas, the mid-latitude region is covered with seasonal sea ice (south of the purple line in Fig. 2), i.e., sea ice that forms in winter and melts in summer.

**Fig. 2: Climatology of anomalous sea ice during the post-MPT (0.0–0.7 Ma) and MPT (0.7–1.2 Ma), and climatological mean during the pre-MPT period (1.2–2.6 Ma).**

Evidently, the dominant periodicity of sea ice extent changes from ∼41 kyr during the pre-MPT period to a longer ∼100 kyr during the post–MPT period (Fig. 3a–c), which is aligned with climate shifts observed across the MPT^21,39 (Supplementary Fig. 3). Spectral analysis of the sea ice extent indicates a latitudinal dependency of sea ice variability throughout the Quaternary (Fig. 3a–c). Meanwhile, the simulated perennial (Fig. 1i) and seasonal sea ice exhibit a strong significant correlation during the Quaternary \(({{{\rm{r}}}}=0.86)\); however, variations are evident in their periodicity modes. To understand the sea ice variability, remembering that the ∼41 kyr (obliquity) and ∼23/19 kyr (precession) signals do not necessarily represent direct responses to orbital forcing is essential. This is because the prescribed CO₂ (GHG) (Supplementary Fig. 2a–c) and ice-sheet (Supplementary Fig. 2j–l) forcings also have a major ∼41 kyr signal during the pre-MPT period and a minor ∼23/19 kyr signal throughout the Quaternary. However, longer periodicity of ∼100 kyr and ∼60–100 kyr (CO₂/ice-sheet signal) would indicate a response to the GHG/ice-sheet forcings, as obliquity and precession have distinct ∼41 kyr and ∼23/19 kyr signals, respectively.

**Fig. 3: Periodicity–variance spectrum of NH sea ice and perennial sea ice.**

Our analysis shows that the pre-MPT high-latitude region covered with perennial sea ice (north of purple lines in Figs. 2c and 3c) has a dominant ∼41 kyr (obliquity) signal (Fig. 3c, f). Meanwhile, the pre-MPT mid-latitude region covered with seasonal sea ice (south of purple lines in Figs. 2c and 3c) has a dominant ∼23 kyr (precession) signal (Fig. 3c). These pre-MPT sea ice signals indicate strong influences of the orbital and GHG forcings, as ice-sheet forcing during the pre-MPT period is weaker than that during the post-MPT period. During the MPT period, perennial sea ice has a dominant ∼41 kyr (obliquity) signal, indicating the influence of the orbital forcing, and a major ∼63 kyr (CO₂/ice–sheet) signal (Fig. 3b, e), indicating the influences of the GHG and ice–sheet forcings. The MPT seasonal sea ice shows varying ∼63–56 kyr (CO₂/ice–sheet), ∼41 kyr (obliquity), and ∼23/19 kyr (precession) signals (Fig. 3b), indicating the influences of the orbital, GHG, and ice-sheet forcings. Further, the post-MPT period has a primary ∼100 kyr (eccentricity) signal throughout the NH latitudes (Fig. 3a), suggesting strong influences of the GHG and ice-sheet forcings. Meanwhile, the post-MPT perennial sea ice (Fig. 3a, d) also shows a secondary ∼41 kyr (obliquity) signal and seasonal sea ice shows a secondary ∼23 kyr (precession) signal, indicating the influence of the orbital forcing. Thus, our results indicate that NH sea ice variability with longer periodicity (∼100 kyr) during the MPT and post-MPT comes from the GHG and ice-sheet forcings which amplified the eccentricity forcing. Additionally, obliquity plays a stronger role in high-latitude perennial sea ice variability, and precession plays a stronger role in mid-latitude seasonal sea ice regions throughout the Quaternary.

Furthermore, the sea ice extent variability over the Atlantic and Pacific Ocean sectors is similar to the NH sea ice signals, particularly that of perennial sea ice (Fig. 4). However, significant differences are evident in the Atlantic sea ice extent. For example, the Atlantic sea ice extent has a dominant ∼23/19 kyr (precession) signal during the pre-MPT and MPT periods (Fig. 4b, c) and a ∼23/19 kyr (precession) signal (∼40% variance) during the post-MPT period (Fig. 4a) with respect to the NH sea ice variability. Consequently, compared to that of the Pacific Ocean sector, the Atlantic Ocean sector shows a stronger influence of the precession forcing.

**Fig. 4: Periodicity–variance spectrum of sea ice in the Atlantic and Pacific Ocean sectors.**

One indicator of the spatiotemporal variability of sea ice is the first Empirical Orthogonal Function (EOF) of sea ice concentration during the post-MPT, MPT, and pre-MPT periods. The first EOF of sea ice shows the major glacial–interglacial sea ice changes (Fig. 5a–c). During the pre-MPT period, the major leading pattern shows sea ice variability in the Labrador Sea, especially along the coast of Newfoundland (Fig. 5c). During the MPT period, an additional major leading pattern is evident in the Irminger–Iceland basin sector, especially along the coast of Norway (Fig. 5b). These two regions have open-ocean deep convection sites that are also evident in this model simulation. Following the MPT (post-MPT) period, three major leading patterns are identified, with an additional sea ice variability emerging over the Central North Pacific (CNP) (Fig. 5a). Such a sequential increase in the number of major leading patterns of sea ice variability is related to the southward expansion of sea ice extent (black solid line in Figs. 2 and 5a–c) associated with a gradual decrease in CO₂ concentration (GHG forcing) (Fig. 1e).

Fig. 5: Glacial–interglacial sea ice variability and periodicity–variance spectrum of sea ice extent in the Labrador Sea (Lr), Irminger–Iceland basin (Ir–Ic), and Central North Pacific (CNP) sectors.

The EOF result reveals three major regions with strong variability in sea ice throughout the Quaternary. The Labrador Sea, Irminger–Iceland basin, and Central North Pacific regions, which are covered with seasonal sea ice, show a dominant ∼23/19 kyr (precession) sea ice signal during the pre-MPT and MPT periods (Fig. 5e, f, h, i, respectively), whereas the Labrador Sea and Irminger–Iceland basin show a ∼23 kyr (precession) signal (∼40% variance) during the post-MPT period (Fig. 5d), which closely resembles the post-MPT Atlantic sea ice variability (Fig. 4a). Our results show that perennial sea ice is sensitive to obliquity (Fig. 3). Meanwhile, as sea ice expands southward (Figs. 2 and 5a–c), the mid-latitude seasonal sea ice shows greater sensitivity to precession. Proxy studies based on marine sediment core MD01–2414, from the mid-latitude region in the central Okhotsk Sea (53° N, 149° E; Fig. 5a), reveal a dominant precession sea ice variability during the last 130 ka³⁰. Meanwhile, sediment core HLY–6, from the high-latitude region in the western Arctic Ocean (78° 18´ N, 177° W; Fig. 5a), shows a significant obliquity sea ice variability and episodes of perennial sea ice coverage²⁸. Additionally, the presence of detrital iron oxide mineral grains in the Arctic Coring Expedition (ACEX) composite core from the Lomonosov Ridge (87.5° N and 138.3° W) indicates that the Arctic has been covered with perennial sea ice over the past 14 Ma⁴⁰. Therefore, the proxy records agree with our findings that the high-latitude Arctic is covered with perennial sea ice and is strongly influenced by obliquity. However, a secondary effect of precession is also evident in the model results during the MPT and post-MPT periods. In contrast, the mid-latitude region is covered with seasonal sea ice and is more sensitive to precession.

Atmospheric and oceanic temperature variability during the Quaternary

Present-day observations, reanalysis data, and the global climate model have shown a tight coupling (negative correlation) between annual mean surface air and ocean temperatures and sea ice⁴¹. The simulated surface air temperature and sea ice extent are significantly (negatively) correlated during all three time periods and the surface air temperature leads sea ice by 1 kyr during the pre-MPT and MPT periods and shows no lag during the post-MPT period (Fig. 6d–f). Notably, the simulated pre-MPT zonal mean surface air temperature shows a dominant ∼41 kyr (obliquity) signal in the high-latitude region, while the mid- and low-latitude regions (0° N–30° N) show both ∼41 kyr (obliquity) and ∼23/19 kyr (precession) signals (Fig. 6c). Also, the cospectrum between surface air temperature and sea ice during the pre-MPT period shows strong ∼41 kyr and ∼23 kyr signals (Fig. 6i), indicating strong influences of the orbital and GHG forcings, as the pre-MPT ice-sheet forcing is weaker than the post-MPT ice-sheet forcing.

**Fig. 6: Periodicity–variance spectrum of NH surface air temperature and its cross-correlation and cospectrum with sea ice extent.**

The MPT surface air temperature shows a major ∼41 kyr (obliquity) signal in the high- and low-latitude regions and a dominant ∼100–56 kyr (CO₂/ice-sheet) signal in the mid- and high-latitude regions (Fig. 6b). Meanwhile, the cospectrum shows a strong CO₂/ice-sheet (∼100–56 kyr) signal alongside ∼41 kyr and ∼23/19 kyr signals (Fig. 6h). The presence of a dominant CO₂/ice-sheet signal during the MPT period suggests a more robust response of the surface air temperature to the GHG and ice-sheet forcings than that of sea ice. Moreover, similar to the case of sea ice variability, the post-MPT period has a dominant ∼100 kyr (eccentricity) signal throughout the NH latitudes, a ∼41 kyr signal in the high-latitude region, and a ∼23 kyr signal in the mid- and low-latitude regions (Fig. 6a), indicating the influences of the orbital, GHG and ice-sheet forcings.

Furthermore, the simulated surface air and ocean temperatures have similar latitudinal periodicity variations (Supplementary Fig. 5a–c) and are significantly (positively) correlated with no phase lag (Supplementary Fig. 5d–f) throughout the Quaternary. Our results suggest that the surface air and ocean surface temperatures are tightly coupled at orbital and millennial timescales. Consequently, oceanic heat released to the atmosphere sustains a strong link between sea ice and surface air temperature variability.

Response to the orbital forcing during the Quaternary

The net surface shortwave radiation over the ocean (positive downward) in the mid-latitude region shows a predominant precession (∼23/19 kyr) signal (Fig. 7a–c), whereas the high-latitude region shows a distinct obliquity (∼41 kyr) signal throughout the Quaternary (Fig. 7a–c). Consequently, these changes in the net surface shortwave radiation signals influence the surface air temperature and sea ice variability. Additionally, the response to the obliquity forcing in the high-latitude region is further amplified by the ice–albedo feedback¹⁷.

Fig. 7: Periodicity–variance spectrum of the mid-latitude and high-latitude net surface shortwave radiation and cross-correlation of summer incoming solar radiation at 65° N, surface air temperature, and sea ice with obliquity and precession indices from 0.0–2.6 Ma.

Next, we applied obliquity and precession bandpass filters (30–70 and 5–30 kyr, respectively) on the simulated summer incoming solar radiation, surface air temperature, and sea ice extent to investigate the responses of sea ice and surface air temperature variability to the orbital forcing (obliquity and precession). The bandpass-filtered summer incoming solar radiation at 65°N is significantly correlated with the obliquity (positively) and precession (negatively) indices. Furthermore, summer incoming solar radiation at 65° N does not lag the obliquity index (purple line in Fig. 7d) during the Quaternary; however, it lags the precession index by 4 kyr (purple line in Fig. 7e). Therefore, an increase in the obliquity index or a decrease in the precession index increases the NH incoming summer solar radiation at 65° N.

Further, the bandpass-filtered high- and mid-latitude surface air temperature are significantly correlated with the obliquity (positively) and precession (negatively) indices, respectively. High-latitude surface air temperature lags the obliquity index by ∼6 kyr (red line in Fig. 7d) during the Quaternary, whereas mid-latitude surface air temperature lags the precession index by ∼8 kyr (red line in Fig. 7e). Meanwhile, the bandpass-filtered sea ice extent is significantly correlated with the obliquity (negatively) and precession (positively) indices and the sea ice extent lags both the obliquity and precession indices by ∼9 kyr (blue line in Fig. 7d, e) during the Quaternary. Therefore, our results show a tight coupling between the simulated summer incoming solar radiation, surface air temperature, and sea ice extent with the orbital forcing throughout the Quaternary. Additionally, high obliquity and/or low precession forcing increases the NH summer incoming solar radiation at 65° N, which increases the surface air/ocean temperature. Consequently, sea ice formation decreases, causing a reduction in sea ice extent and vice versa, consistent with previous studies^17,29,30,42.

Specifically, perennial and seasonal sea ice are significantly correlated with the obliquity (negatively) and precession (positively) indices, respectively. Perennial sea ice lags the obliquity index by ∼5 kyr (green line in Fig. 7d), and seasonal sea ice lags the precession index by ∼10 kyr (green line in Fig. 7e). Hence, the NH mid-latitude region covered with seasonal sea ice is strongly influenced by the precession forcing, whereas the high-latitude region covered with perennial sea ice is more strongly influenced by the obliquity forcing throughout the Quaternary.

Furthermore, the Earth’s orbital forcing has a ∼100 kyr (eccentricity) signal which modulates the strength of the precessional cycle. Naturally or in a more realistic model setup with dynamical carbon cycle and ice-sheets, the ∼100 kyr signal in post-MPT sea ice and surface air temperature would have been explained by eccentricity, amplified by feedbacks involving GHG and ice-sheets. However, without the amplification of these physical quantities through feedbacks within the Earth’s climate system, the variability of sea ice due to the orbital forcing alone is likely to be relatively weaker. Meanwhile, this experiment treats GHG and ice-sheet as external forcing. Therefore, a more reasonable interpretation of the model results is that the signal of the orbital forcing was amplified through feedbacks within the climate system, and these consequent GHG and ice-sheet changes more effectively influenced sea ice variability.

Response to the GHG forcing during the Quaternary

The Earth’s surface and lower atmosphere are warming in response to increased GHG concentrations. As the temperatures increase, more water vapor is released into the atmosphere, intensifying the greenhouse effect by causing increased emissions of longwave radiation from the atmosphere to the surface, resulting in additional warming⁴³ (i.e., the water vapor feedback⁴⁴). Cloud cover also affects downwelling longwave radiation⁴³. Here, we show that the changes in prescribed CO₂ concentration can explain the major regions of glacial–interglacial sea ice variability (Fig. 8a–c), as shown in Fig. 5a–c. Furthermore, changes in downwelling longwave radiation (clear-sky conditions) caused by GHG and low-level water vapor changes explain the glacial–interglacial major sea ice regional changes in the (a) Labrador Sea, Irminger–Iceland basin and Central North Pacific during the post-MPT period (Fig. 8d), (b) Labrador Sea and Irminger–Iceland basin during the MPT period (Fig. 8e), and Labrador Sea during the pre-MPT period (Fig. 8f) as shown in Fig. 5a–c. Additionally, changes in downwelling longwave radiation explain the intensification of sea ice extent during the MPT and post-MPT periods (Fig. 8d, e).

**Fig. 8: Relationship between the CO₂ (GHG) forcing and sea ice during the quaternary.**

In contrast to the orbital forcing, the GHG forcing has a uniform zonal distribution. The prescribed CO₂ (GHG) concentration (Supplementary Fig. 2a–c) and the simulated GMST (Supplementary Fig. 6) have similar periodicity signals throughout the Quaternary. Meanwhile, the pre-MPT prescribed CO₂ concentration (Fig. 1e) and the simulated GMST (Fig. 1h) show gradual (significant) downward trends, which intensify across the MPT by approximately two and three times, respectively. In contrast, the simulated sea ice extent shows an upward (significant) trend during the pre-MPT (\(1.5\) thousand km²/kyr), which also intensifies by a factor of approximately 2.2 (\(3.4\) thousand km²/kyr) during the MPT period (Fig. 1g), indicating the influence of the GHG forcing on sea ice. Proxy sea ice records from sediment core U1343 (57° 33.4′ N, 176° 49.0′ W; 1,950 m) located in the Bering Sea shelf region are in agreement with our result, highlighting a significant increase in sea ice extent throughout the MPT period^27,45.

The prescribed CO₂ concentration (Fig. 1e) and simulated GMST (Fig. 1h) are significantly (positively) correlated with no lag during the Quaternary (Fig. 8g). Consequently, this indicates that as GHGs and GMST increase simultaneously, warmer summers and enhanced sea ice melting occur. Furthermore, the NH sea ice extent with GMST and high-latitude/perennial sea ice with high-latitude surface air temperature are significantly (negatively) correlated with no lag during the Quaternary (Fig. 8g). Therefore, our results show that compared to that during the post-MPT period, the higher (26%) prescribed GHG forcing during the pre-MPT increased GMST by 2.2 °C (21%), causing warmer seasons. Consequently, the pre-MPT period would have enhanced summertime sea ice melting and reduced wintertime sea ice formation, causing a comparably smaller extent of sea ice (Fig. 2c). In contrast, during the post-MPT period, the relatively lower GHG forcing and GMST facilitate a comparably greater extent of sea ice (Fig. 2a) because of somewhat colder seasons.

Furthermore, our analysis shows that the GHG forcing can largely explain the simulated NH sea ice variability, while a secondary influence of precession forcing is also evident during the Quaternary (Table 1). Notably, the obliquity forcing strongly influences high-latitude/perennial sea ice, and a secondary influence of the precession forcing cannot be ignored. Additionally, the orbital indices lead to sea ice change, while GHG changes lag behind sea ice change, suggesting that changes in the orbital forcing initiate sea ice change, which is then amplified by the GHG forcing.

Table 1 Standardized regression coefficients from a lag multiple linear regression of NH sea ice extent and high-latitude/perennial sea ice with the GHG forcing, obliquity, and precession indices during the Quaternary

Full size table

Discussion

Previous proxy^{28,30,33,46,47,48,49} and modeling^{29,30,31,32,42,50} studies mostly recorded or simulated sea ice concentrations limited to the post-MPT period. A recent study (i.e., the 3 Ma experiment) has presented a broad outline of the Northern and Southern Hemisphere sea ice variability, showing a ∼21 kyr (precession), ∼41 kyr (obliquity), and ∼60–120 kyr (CO₂/ice-sheet) periodicity, suggesting the combined effect of the GHG and ice-sheet forcings³⁴ on sea ice variability. Here, we analyze the 3 Ma experiment to investigate the spatiotemporal and latitudinal sea ice variability as well as the mechanism of NH sea ice variability during the post-MPT, MPT, and pre-MPT periods. Net surface shortwave radiation over the mid-latitude region has a prominent precession signal in response to the precession forcing, as precession strongly influences daily incoming solar radiation in the mid-latitudes and tropics^35,51,52,53. Meanwhile, net surface shortwave radiation over the high-latitude region has a dominant obliquity signal in response to the obliquity forcing, as obliquity strongly influences daily incoming solar radiation in the high-latitude region¹⁷. In addition, the geographical features of the NH high-latitude region, which is mostly limited to the Arctic Ocean, influence the sea ice³⁴. Therefore, our study highlights the latitudinal distribution of net surface shortwave radiation (also affected by the sea ice–albedo feedback), which is strongly influenced by the orbital forcing and influences the tightly coupled surface air temperature and sea ice variability (Fig. 9).

**Fig. 9: Schematic showing how the orbital, GHG, and ice-sheet forcings influence the NH sea ice.**

Furthermore, our study shows that the glacial–interglacial sea ice varies predominantly in the mid-latitude region covered with seasonal sea ice (Fig. 2) and is strongly influenced by the precession (Fig. 3a–c). Meanwhile, the high-latitude region is covered with perennial sea ice (Fig. 2), which is strongly influenced by the obliquity (Fig. 3d–f). Our results agree with those of previous studies, that showed that sea ice and surface air temperature are strongly correlated with precession^28,29,30 and obliquity^{17,28,31,32,33}. Proxy studies indicate the presence of perennial sea ice in the Arctic over the past 14 Ma⁴⁰, in the high-latitude region after ∼0.4 Ma²⁸, and in the central Arctic Ocean during the Late Quaternary (penultimate glacial and Last Interglacial) period⁴⁶. Meanwhile, other records indicate seasonal sea ice during the pre-MPT period⁴⁹ and before \(\sim\)0.35 Ma⁴⁸. Our results mostly agree with available sediment proxy records, which generally suggest that perennial sea ice covered the central Arctic basin (i.e., high-latitude region in this study) during most—if not all—of the Quaternary⁵⁴.

Yin and Berger¹⁷ identified a linear relationship between simulated GMST and the GHG forcing over the past 0.8 Ma. Here, we also show that the GHG forcing and GMST have been strongly correlated over the past 2.6 Ma. Consequently, an increase (decrease) in CO₂ concentration increases (decreases) the GMST, which decreases (increases) the NH sea ice extent (Fig. 8g). Moreover, northward oceanic heat transport impacts the NH sea ice extent, albeit to a smaller extent compared to the Southern Hemisphere^34,55. Further, longer periodicity signal (∼100 kyr and ∼60–100 kyr) from the GHG and ice-sheet forcings is most noticeable in the NH sea ice during the MPT and post-MPT periods (Fig. 3a–c). Evidently, the GHG forcing together with the water vapor feedback⁴⁴ influence the downwelling longwave radiation, which explains the majority of the glacial–interglacial sea ice variability in the Labrador Sea, Irminger–Iceland basin and Central North Pacific regions during the MPT and post-MPT periods (Fig. 8d–f). The Labrador Sea and Irminger–Iceland basin are deep convection regions where deep water is formed; hence, sea ice could influence the Atlantic Meridional Overturning Circulation (AMOC) strength and variability through the sea ice–ocean feedback⁵⁶. Additionally, sea ice loss results in increased evaporation from the ocean surface, resulting in increased atmospheric water vapor levels⁵⁷.

The strong association between sea ice and GHG forcing observed in the 3 Ma experiment could be because the GHG forcing anomalies relative to preindustrial conditions are amplified by a factor of 1.5 to compensate for the relatively weak climate sensitivity in the applied CESM1.2 model³⁴ (see Methods for details). Moreover, as mentioned in the Introduction, because the GHGs and ice-sheets are treated as external forcings, the model cannot simulate changes resulting from interactions among the orbital, GHG, and ice-sheets. Therefore, a more reasonable interpretation of the model findings is that the orbital forcing signal was intensified through feedbacks within the climate system, and these associated alterations in GHGs and ice-sheets more effectively influenced sea ice variability.

Additionally, changes in the prescribed NH ice-sheet topography and extent would influence surface albedo, alter the meridional temperature gradient by influencing the zonal distribution of temperature change⁵⁸, and modulate the NH westerlies¹⁸ and Intertropical Convergence Zone⁵⁸. Previous studies have shown that the intensification of NH glaciation (ice-sheet) occurred across the MPT (Fig. 1f) and Pliocene-Pleistocene Transition (PPT; 2.5–3.2 Ma)^20,37,39. The Northern Annular Mode, which characterizes the influence of ice-sheet on the atmospheric circulation, and the simulated high-latitude/perennial sea ice shows an intensification across the MPT and PPT (Supplementary Fig. 7). Meanwhile, the glacial–interglacial variability of the simulated NH net surface shortwave radiation indicates a strong ice–sheet influence through the albedo feedback and shows an intensification in magnitude and coverage over northern North America during the post-MPT and MPT periods (Supplementary Fig. 8a, b). Furthermore, the prescribed sea level change, which is equivalent to ice-sheet forcing¹⁹, influences sea ice variability across the Arctic continental shelves, particularly the Labrador Sea region, during the pre-MPT period (Supplementary Fig. 9c). This is because sea ice expanded less during the pre-MPT period than during the MPT and post-MPT periods, causing it to vary strongly along the continental shelves. Therefore, our study documents that the NH ice-sheet forcing together with orbital and GHG forcings strongly influenced the NH sea ice variability throughout the Quaternary (Fig. 9).

The 3 Ma experiment has a few caveats, considering its unprecedented length. For instance, the experiment employs a relatively small acceleration technique with a factor of five which is applied to the external (i.e., prescribed orbital, GHG, and ice-sheet) forcings, compared with the acceleration factor of 10 used in previous studies^59,60. Consequently, it reduces the 3 Ma orbital year forcing to 0.6 Ma model years. Acceleration-induced biases affect the NH high-latitude ocean surface temperature and sea ice^59,60,61, with the biases depending on the acceleration factors⁶¹. For instance, Lunt, et al.⁶¹ examined numerous aspects of the ocean–atmosphere system over the past 30,000 years with different acceleration factors (1, 2, 5, 10). They found that most of the bias is in the Southern Ocean, whereas the NH is relatively unaffected even with an acceleration factor of 10. Additionally, an acceleration factor of five, which is employed in the 3 Ma experiment, is a good way to minimize the resulting biases in ocean surface temperature and sea ice variability in the high-latitude region where the surface climate is closely associated with the deep ocean via deep convection. Next, the AMOC evolution shows relatively small changes due to the acceleration factors^60,61. An acceleration factor of five would allow the deep ocean to adjust on advective timescales. However, this would still lead to biases in deep ocean variables for the precession timescale but less bias for obliquity and eccentricity timescales³⁴. Furthermore, Yun, et al.³⁴ compared several 3 Ma simulated variables with a range of proxies, showing that, at least for near-surface variables, the biases in the 3 Ma experiment are likely very small.

Another caveat is that the 3 Ma experiment does not consider the ice-sheet meltwater flux into the ocean. Also, the 3 Ma experiment uses a constant ocean bathymetry and land–sea boundary from the Last Glacial Maximum (∼21 kyr), resulting in a sea level ∼135 m lower than the present sea level¹⁹. Yun, et al.³⁴ note that ignoring the ice-sheet meltwater flux likely causes an underestimation of the AMOC variability, which has a strong precession signal and would change the contribution of the ocean heat transport to the meridional heat transport. A constant Last Glacial Maximum sea level also causes some unrealistic land–sea distributions in the high-latitude region during interglacials, such as a closed Bering Strait and subaerial Barents Shelf (Fig. 2c). A closed Bering Sea prevents ocean heat and water transport between the Arctic and the Pacific via the Bering Strait, which could lead to the emergence of stronger hysteresis behavior in the AMOC during glacial periods⁶². Furthermore, a closed Bering Sea would hinder ocean heat transport into the Pacific high-latitude region compared with the Atlantic and influence regional sea ice variability. Therefore, sea ice variability in the Atlantic Ocean sector has a stronger AMOC-influenced precession signal than that in the Pacific Ocean sector (Fig. 4). Additionally, the simulated sea ice extent in the Pacific Ocean sector south of the Bering Strait (Supplementary Fig. 10) and the entire Pacific Ocean (Fig. 4d–f) have similar seasonal sea ice variability akin to that of the NH sea ice (Fig. 3), possibly because the simulated NH meridional heat is transported equatorward, mainly through the atmosphere, for high obliquity and high atmospheric CO₂ levels³⁴. Additionally, a subaerial Barents Shelf reduces sea ice formation during glacial–interglacial periods, influences surface albedo changes⁵⁸, and warms the high-latitude region⁵⁸. Then again, using a volume-conserving ocean model with time-evolving ocean volume and sea level would increase computational costs and likely result in unidentified biases and ambiguities that have not been fully explored³⁴.

The 3 Ma experiment employed orbital forcing from Berger³⁵, which underestimated the magnitude of June incoming solar radiation at 65° N but strongly (positively) correlated with later numerical estimates by Laskar, et al.⁵³ during the pre-MPT period (Supplementary Fig. 11). Since this study focuses on sea ice variability over a longer timescale, this discrepancy in the orbital estimates is unlikely to significantly impact glacial–interglacial sea ice variability during the pre-MPT period. Nonetheless, future studies should assess changes in the climate system more comprehensively. Furthermore, the experiment does not explicitly simulate aerosols, dust, and climate–carbon cycle feedbacks, while it includes a prescribed vegetation component where the vegetation changes only in the land ice-covered areas. Vegetation cover affects albedo and evapotranspiration and may significantly influence the Arctic surface energy budget^63,64. Therefore, to attain a more realistic response to longwave radiation and consider the radiative impacts of CH₄ and N₂O concentrations (Supplementary Fig. 2d–i), vegetation, and dust, the GHG forcing anomalies relative to preindustrial conditions are amplified by a factor of 1.5 (see Methods for details). Additionally, the lack of an interactive carbon cycle fails to resolve the interaction between sea ice and the outgassing of CO₂ from the ocean to the atmosphere on a glacial–interglacial timescale⁶⁵.

Methods

Model and simulation

The 3 Ma experiment³⁴ is a transient paleoclimate model simulation covering the past three million years, employing the coupled CESM version 1.2 (CESM1.2). The experiment comprises fully coupled atmosphere (Community Atmosphere Model version 4; CAM4), ocean (Parallel Ocean Program version 2; POP2), sea ice (Los Alamos Sea Ice Model; CICE4), and land (Community Land Model version 4; CLM4.0; prescribed vegetation) components. The atmospheric model has latitude–longitude horizontal resolution (T31) of approximately 3.75° × 3.75° with 26 vertical levels. The ocean and sea ice models have a 100 × 116 horizontal grid (gx3v5), with the ocean model comprising 25 vertical levels. The ocean and sea ice output data are regridded and analyzed at a latitude–longitude horizontal resolution of 3.75° × 3.75°, identical to that of the atmospheric outputs. NCAR Command Language (NCL)⁶⁶ version 6.6.2 was used to analyze data and generate all figures.

The 3 Ma experiment is forced by transient prescribed orbital³⁵ forcing, which changes the astronomical insolation, ice-sheet¹⁹ forcing, which changes the NH ice-sheet extent and topography, while the Southern Hemisphere ice-sheet is kept constant, and GHG (including CO₂, CH₄, and N₂O concentrations)^19,21 forcing, which changes the radiative forcing. The GHG forcing includes CO₂ concentration estimates from various datasets. From 0.8 to 3.0 Ma, the estimated GHGs consist of simulated CO₂ concentration from a previous CLIMBER-2 transient simulation¹⁹ containing an interactive carbon cycle. Meanwhile, from 0.0 to 0.8 Ma, the GHGs estimate includes CO₂ concentration measurements from air bubbles from the European Project for Ice Coring in Antarctica (EPICA) Dome C ice core⁶⁷. The GHG forcing also includes N₂O and CH₄ concentration estimates as a CESM forcing obtained by regressing the EPICA Dome C ice core N₂O and CH₄ concentrations on the LR04 stack²¹ deep sea proxy record post 0.8 Ma and extending the linear statistical model to 3 Ma using the LR04 stack deep sea record. The LR04 stack²¹ deep sea proxy record is constructed by the graphic correlation of an average of globally distributed 57 benthic foraminiferal calcite \({{{\rm{\delta }}}}{{{{\rm{O}}}}}^{18}\) records for the past 5.3 million years (unit is in per mil), and it measures the global ice volume and deep ocean temperature. The ice-sheet topography changes are obtained from the simulated CLIMBER-2 orography data from the previous transient simulation¹⁹. The 3 Ma experiment also uses constant ocean bathymetry from the Last Glacial Maximum (∼21 ka) and a prescribed vegetation component, where the plant functional types (i.e., vegetation) are changed only in the NH land ice-covered areas.

The 3 Ma simulation uses an acceleration factor of five for the external (orbital, ice-sheet, and GHG) forcings, which reduces the 3 Ma orbital year forcing to 0.6 Ma model years. The orbital, ice-sheet, and GHG forcings are updated every 100 orbital years, corresponding to 20 model years. The CESM1.2 model in the 3 Ma experiment resolution exhibits a comparatively low standard equilibrium climate sensitivity of 2.4°C warming for a two-fold CO₂ increase with respect to the preindustrial CO₂ concentration. This is approximately 1.5 times smaller than the likely range of estimates (3.7 ± 1.2 °C) obtained within 37 other climate model simulations conducted as part of the Coupled Model Intercomparison Project Phase 6 (CMIP6)⁶⁸. Therefore, to achieve a more realistic response to past longwave radiative forcings and indirectly account for the radiative impacts of other CO₂-correlated forcings from vegetation, dust, N₂O, or CH₄, the CO₂ forcing¹⁹ anomalies relative to preindustrial conditions are amplified by a factor of 1.5 in the 3 Ma experiment. The resulting effective equilibrium climate sensitivity is ~3.8 °C, which is quite consistent with the CMIP6 estimate. Additionally, it is in reasonable agreement with prior paleoclimate estimates (i.e., 3.2 °C)^69,70, which were derived from reconstructions of GMST and radiative forcings during the past 0.784 Ma. The re-scaling of the CO₂ forcing in CESM1.2 resulted in a more realistic representation of the magnitude of tropical, global mean, and Antarctic temperature changes^34,36.

Furthermore, the 3 Ma experiment does not explicitly simulate the climate–carbon cycle feedback. In reality, the GHG forcing is a response and byproduct of the feedback from the carbon cycle. Therefore, an interactive carbon cycle model would be more realistic. However, under the current model’s capability, including a carbon cycle model, may produce other model biases. Therefore, by prescribing the GHG forcing, which involves CO₂, CH₄, and N₂O concentrations, we can assess its impact on sea ice without invoking a complex causality in the Earth’s carbon cycle system. Additionally, because the GHG and ice-sheet forcings are treated as external climate drivers, the interactive feedback processes between the GHGs and ice-sheets were not incorporated. Incorporating these interactions would provide a more physically realistic approach. However, achieving a balanced representation of climate change associated with the natural carbon cycle and the climate feedback resulting from ice-sheet dynamics presents significant challenges at this stage. Therefore, the initial approach was to consider these factors as climate forcings to minimize potential biases in the 3 Ma experiment.

The 3 Ma simulation data were compared against large-scale climate reconstructions and various paleo–proxy datasets. The 3 Ma simulation has been shown to reproduce the magnitude and timing of glacial–interglacial variability^34,36 consistent with paleoclimate data constraints and estimates from the intermediate complexity Earth system model. The 3 Ma simulation agrees well^34,36 with key paleoclimate proxy records, such as the magnitude and timing of GMST variations, annual mean ocean surface temperatures over the tropics and elsewhere (western South and North Pacific, southern Indian, and North Atlantic), Greenland and Antarctic temperatures, the East Asian summer monsoon, and the East African hydroclimate.

The 3 Ma simulations are available in annual average millennial mean data. Yun, et al.³⁴ provide the complete details on the 3 Ma model set-up and experimental design.

Data availability

The 3 Ma experiment data are available on the OPeNDAP ICCP climate data server https://climatedata.ibs.re.kr (3 Ma transient simulation with CESM1.2; dataset https://doi.org/10.22741/iccp.20230001).

Code availability

All codes used for the analysis in this study are available upon request from G.M. (gagan_mandal@yonsei.ac.kr). NCAR Command Language (NCL) version 6.6.2 was used to generate all figures.

References

Pithan, F. & Mauritsen, T. Arctic amplification dominated by temperature feedbacks in contemporary climate models. Nat. Geosci. 7, 181–184 (2014).
Article ADS CAS MATH Google Scholar
Serreze, M. C. & Barry, R. G. Processes and impacts of Arctic amplification: a research synthesis. Glob. Planet. Change 77, 85–96 (2011).
Article ADS MATH Google Scholar
Kidston, J., Taschetto, A. S., Thompson, D. W. J. & England, M. H. The influence of Southern Hemisphere sea-ice extent on the latitude of the mid-latitude jet stream. Geophys. Res. Lett. 38, L15804 (2011).
Article ADS Google Scholar
Abernathey, R. P. et al. Water-mass transformation by sea ice in the upper branch of the Southern Ocean overturning. Nat. Geosci. 9, 596–601 (2016).
Article ADS CAS MATH Google Scholar
Goosse, H., Roche, D. M., Mairesse, A. & Berger, M. Modelling past sea ice changes. Quat. Sci. Rev. 79, 191–206 (2013).
Article ADS MATH Google Scholar
Bretones, A., Nisancioglu, K. H., Jensen, M. F., Brakstad, A. & Yang, S. Transient increase in Arctic deep-water formation and ocean circulation under sea ice retreat. J. Clim. 35, 109–124 (2022).
Article ADS Google Scholar
Notz, D. & Stroeve, J. Observed Arctic sea-ice loss directly follows anthropogenic CO₂ emission. Science 354, 747–750 (2016).
Article ADS CAS PubMed Google Scholar
Kinnard, C. et al. Reconstructed changes in Arctic sea ice over the past 1,450 years. Nature 479, 509–512 (2011).
Article ADS CAS PubMed MATH Google Scholar
Brennan, M. K., Hakim, G. J. & Blanchard‐Wrigglesworth, E. Arctic sea‐ice variability during the instrumental era. Geophys. Res. Lett. 47, e2019GL086843 (2020).
Article ADS Google Scholar
Cavalieri, D. J. & Parkinson, C. L. Arctic sea ice variability and trends, 1979-2010. Cryosphere 6, 881–889 (2012).
Article ADS MATH Google Scholar
Walsh, J. E., Fetterer, F., Scott stewart, J. & Chapman, W. L. A database for depicting Arctic sea ice variations back to 1850. Geogr. Rev. 107, 89–107 (2017).
Article Google Scholar
Masson-Delmotte, V. et al. Information from paleoclimate archives. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Doschung, A. Nauels, Y. Xia, V. Bex, and P. M. Midgley. (Cambridge University Press, 2013). https://doi.org/10.1017/cbo9781107415324.013.
Hays, J. D., Imbrie, J. & Shackleton, N. J. Variations in the Earth’s orbit: pacemaker of the ice ages. Science 194, 1121–1132 (1976).
Article ADS CAS PubMed MATH Google Scholar
Imbrie, J. et al. On the structure and origin of major glaciation cycles 1. Linear responses to Milankovitch Forcing. Paleoceanography 7, 701–738 (1992).
Article ADS MATH Google Scholar
Milankovitch, M. Kanon der Erdbestrahlung und seine Anwendung auf des Eiszeitenproblem. Special publication 132, section of mathematical and natural sciences. R. Serb. Acad. Sci. 33, 633 (1941).
MathSciNet MATH Google Scholar
Abe-Ouchi, A. et al. Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume. Nature 500, 190–193 (2013).
Article ADS CAS PubMed MATH Google Scholar
Yin, Q. Z. & Berger, A. Individual contribution of insolation and CO₂ to the interglacial climates of the past 800,000 years. Clim. Dyn. 38, 709–724 (2012).
Article MATH Google Scholar
Oster, J. L., Ibarra, D. E., Winnick, M. J. & Maher, K. Steering of westerly storms over western North America at the Last Glacial Maximum. Nat. Geosci. 8, 201–205 (2015).
Article ADS CAS Google Scholar
Willeit, M., Ganopolski, A., Calov, R. & Brovkin, V. Mid-Pleistocene transition in glacial cycles explained by declining CO(2) and regolith removal. Sci. Adv. 5, eaav7337 (2019).
Article ADS CAS PubMed PubMed Central Google Scholar
Ahn, S., Khider, D., Lisiecki, L. E. & Lawrence, C. E. A probabilistic Pliocene–Pleistocene stack of benthic δ18O using a profile hidden Markov model. Dyn. Stat. Clim. Syst. 2 https://doi.org/10.1093/climsys/dzx002 (2017).
Lisiecki, L. E. & Raymo, M. E. A. Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanography 20, PA1003 (2005).
ADS MATH Google Scholar
McClymont, E. L., Sosdian, S. M., Rosell-Melé, A. & Rosenthal, Y. Pleistocene sea-surface temperature evolution: Early cooling, delayed glacial intensification, and implications for the mid-Pleistocene climate transition. Earth-Sci. Rev. 123, 173–193 (2013).
Article ADS CAS Google Scholar
Panieri, G., Knies, J., Vadakkepuliyambatta, S., Lee, A. L. & Schubert, C. J. Evidence of Arctic methane emissions across the mid-Pleistocene. Commun. Earth Environ. 4, 109 (2023).
Article ADS Google Scholar
Clark, P. U. et al. The middle Pleistocene transition: characteristics, mechanisms, and implications for long-term changes in atmospheric pCO₂. Quat. Sci. Rev. 25, 3150–3184 (2006).
Article ADS MATH Google Scholar
Berends, C. J., Köhler, P., Lourens, L. J. & Wal, R. S. W. On the cause of the mid‐Pleistocene transition. Rev. Geophys. 59, e2020RG000727 (2021).
Article ADS Google Scholar
Elderfield, H. et al. Evolution of ocean temperature and ice volume through the mid-Pleistocene climate transition. Science 337, 704–709 (2012).
Article ADS CAS PubMed MATH Google Scholar
Detlef, H. et al. Sea ice dynamics across the mid-Pleistocene transition in the Bering Sea. Nat. Commun. 9, 941 (2018).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Cronin, T. M. et al. A 600-ka Arctic sea-ice record from Mendeleev Ridge based on ostracodes. Quat. Sci. Rev. 79, 157–167 (2013).
Article ADS MATH Google Scholar
Wu, Z., Yin, Q., Guo, Z. & Berger, A. Hemisphere differences in response of sea surface temperature and sea ice to precession and obliquity. Glob. Planet. Change 192, 103223 (2020).
Article Google Scholar
Lo, L. et al. Precession and atmospheric CO₂ modulated variability of sea ice in the central Okhotsk Sea since 130,000 years ago. Earth Planet. Sci. Lett. 488, 36–45 (2018).
Article ADS MATH Google Scholar
Mantsis, D. F., Clement, A. C., Broccoli, A. J. & Erb, M. P. Climate feedbacks in response to changes in obliquity. J. Clim. 24, 2830–2845 (2011).
Article ADS Google Scholar
Tuenter, E., Weber, S. L., Hilgen, F. J. & Lourens, L. J. Sea-ice feedbacks on the climatic response to precession and obliquity forcing. Geophys. Res. Lett. 32, L24704 (2005).
Article ADS Google Scholar
Turney, C. S. M. et al. Obliquity‐driven expansion of North Atlantic sea ice during the last glacial. Geophys. Res. Lett. 42, 10382–10390 (2015).
Article ADS MATH Google Scholar
Yun, K.-S. et al. A transient coupled general circulation model (CGCM) simulation of the past 3 million years. Climate 19, 1951–1974 (2023).
MATH Google Scholar
Berger, A. L. Long-term variations of caloric insolation resulting from the Earth’s orbital elements. Quat. Res. 9, 139–167 (1978).
Article MATH Google Scholar
Timmermann, A. et al. Climate effects on archaic human habitats and species successions. Nature 604, 495–501 (2022).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Willeit, M., Ganopolski, A., Calov, R., Robinson, A. & Maslin, M. The role of CO₂ decline for the onset of Northern Hemisphere glaciation. Quat. Sci. Rev. 119, 22–34 (2015).
Article ADS Google Scholar
Ditlevsen, P. The Pleistocene glacial cycles and millennial-scale climate variability. Atmos. Ocean 60, 233–244 (2022).
Article ADS MATH Google Scholar
Ruggieri, E., Herbert, T., Lawrence, K. T. & Lawrence, C. E. Change point method for detecting regime shifts in paleoclimatic time series: application to δ18O time series of the Plio-Pleistocene. Paleoceanography 24, PA1204 (2009).
Article ADS MATH Google Scholar
Darby, D. A. Arctic perennial ice cover over the last 14 million years. Paleoceanography 23, PA1S07 (2008).
Article MATH Google Scholar
Olonscheck, D., Mauritsen, T. & Notz, D. Arctic sea-ice variability is primarily driven by atmospheric temperature fluctuations. Nat. Geosci. 12, 430–434 (2019).
Article ADS CAS Google Scholar
Wu, Z., Yin, Q., Guo, Z. & Berger, A. Comparison of Arctic and Southern Ocean sea ice between the last nine interglacials and the future. Clim. Dyn. 59, 519–529 (2022).
Article Google Scholar
Stephens, G. L. et al. The global character of the flux of downward longwave radiation. J. Clim. 25, 2329–2340 (2012).
Article ADS MATH Google Scholar
Held, I. & Soden, B. Water vapor feedback and global warming. Annu. Rev. Energy Environ. 25, 441–475 (2000).
Article MATH Google Scholar
Teraishi, A., Suto, I., Onodera, J. & Takahashi, K. Diatom, silicoflagellate and ebridian biostratigraphy and paleoceanography in IODP 323 Hole U1343E at the Bering slope site. Deep Sea Res. Part II Topical Stud. Oceanogr. 125-126, 18–28 (2016).
Article ADS Google Scholar
Stein, R., Fahl, K., Gierz, P., Niessen, F. & Lohmann, G. Arctic Ocean sea ice cover during the penultimate glacial and the last interglacial. Nat. Commun. 8, 373 (2017).
Article ADS PubMed PubMed Central Google Scholar
Sha, L. et al. Solar forcing as an important trigger for West Greenland sea-ice variability over the last millennium. Quat. Sci. Rev. 131, 148–156 (2016).
Article ADS MATH Google Scholar
Cronin, T. M. et al. Enhanced Arctic amplification began at the mid-Brunhes event ~400,000 years ago. Sci. Rep. 7, 14475 (2017).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Polyak, L., Best, K. M., Crawford, K. A., Council, E. A. & St-Onge, G. Quaternary history of sea ice in the western Arctic Ocean based on foraminifera. Quat. Sci. Rev. 79, 145–156 (2013).
Article ADS Google Scholar
Yin, Q. Z. & Berger, A. Insolation and CO₂ contribution to the interglacial climate before and after the Mid-Brunhes event. Nat. Geosci. 3, 243–246 (2010).
Article ADS CAS MATH Google Scholar
Tigchelaar, M. & Timmermann, A. Mechanisms rectifying the annual mean response of tropical Atlantic rainfall to precessional forcing. Clim. Dyn. 47, 271–293 (2016).
Article MATH Google Scholar
Merlis, T. M., Schneider, T., Bordoni, S. & Eisenman, I. Hadley circulation response to orbital precession. Part II: subtropical continent. J. Clim. 26, 754–771 (2013).
Article ADS Google Scholar
Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. Astron. Astrophys. 428, 261–285 (2004).
Article ADS MATH Google Scholar
de Vernal, A. & Hillaire-Marcel, C. Quaternary Arctic sea-ice cover: Mostly perennial with seasonal openings during interglacials. Glob. Chang. Mag. 30, 94–95 (2022).
Google Scholar
Goosse, H. & Renssen, H. A two-phase response of the Southern Ocean to an increase in greenhouse gas concentrations. Geophys. Res. Lett. 28, 3469–3472 (2001).
Article ADS CAS MATH Google Scholar
Mandal, G., Hettiarachchi, A. I. & Ekka, S. V. The North Atlantic subpolar ocean dynamics during the past 21,000 years. Dyn. Atmos. Oceans 106, 101462 (2024).
Article Google Scholar
Ridley, J. K., Blockley, E. W. & Ringer, M. A. Arctic Sea Ice Causes Seasonal Differences in the Response of Arctic Water Vapor to Climate Warming in the CMIP6 Model, HadGEM3‐GC3.1. Geophys. Res. Lett. 50 https://doi.org/10.1029/2022gl102541 (2023).
Cao, J., Wang, B. & Liu, J. Attribution of the Last Glacial Maximum climate formation. Clim. Dyn. 53, 1661–1679 (2019).
Article MATH Google Scholar
Timm, O. & Timmermann, A. Simulation of the Last 21 000 Years Using Accelerated Transient Boundary Conditions. J. Clim. 20, 4377–4401 (2007).
Article ADS MATH Google Scholar
Varma, V., Prange, M. & Schulz, M. Transient simulations of the present and the last interglacial climate using the Community Climate System Model version 3: effects of orbital acceleration. Geosci. Model Dev. 9, 3859–3873 (2016).
Article ADS MATH Google Scholar
Lunt, D. J., Williamson, M. S., Valdes, P. J., Lenton, T. M. & Marsh, R. Comparing transient, accelerated, and equilibrium simulations of the last 30,000 years with the GENIE-1 model. Clim 2, 221–235 (2006).
ADS Google Scholar
Hu, A. et al. Role of the Bering Strait on the hysteresis of the ocean conveyor belt circulation and glacial climate stability. Proc. Natl Acad. Sci. 109, 6417–6422 (2012).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Otto-Bliesner, B. L., Lofverstrom, M., Bakker, P. & Feng, R. Arctic warming and the Greenland ice sheet during the last interglacial. Glob. Chang. Mag. 27, 22–23 (2019).
Google Scholar
Swann, A. L., Fung, I. Y., Levis, S., Bonan, G. B. & Doney, S. C. Changes in Arctic vegetation amplify high-latitude warming through the greenhouse effect. Proc. Natl Acad. Sci. USA 107, 1295–1300 (2010).
Article ADS CAS PubMed PubMed Central Google Scholar
Stein, K., Timmermann, A., Kwon, E. Y. & Friedrich, T. Timing and magnitude of Southern Ocean sea ice/carbon cycle feedbacks. Proc. Natl Acad. Sci. USA 117, 4498–4504 (2020).
Article ADS CAS PubMed PubMed Central Google Scholar
The NCAR Command Language (Version 6.6.2)[Software]. (2019). Boulder, Colorado: UCAR/NCAR/CISL/TDD. https://doi.org/10.5065/D6WD3XH5.
Jouzel, J. et al. Orbital and millennial Antarctic climate variability over the past 800,000 years. Science 317, 793–796 (2007).
Article ADS CAS PubMed Google Scholar
Meehl, G. A. et al. Context for interpreting equilibrium climate sensitivity and transient climate response from the CMIP6 Earth system models. Sci. Adv. 6, eaba1981 (2020).
Article ADS PubMed PubMed Central Google Scholar
Friedrich, T., Timmermann, A., Tigchelaar, M., Elison Timm, O. & Ganopolski, A. Nonlinear climate sensitivity and its implications for future greenhouse warming. Sci. Adv. 2, e1501923 (2016).
Article ADS PubMed PubMed Central Google Scholar
Friedrich, T. & Timmermann, A. Using Late Pleistocene sea surface temperature reconstructions to constrain future greenhouse warming. Earth Planet. Sci. Lett. 530, 115911 (2020).
Article CAS Google Scholar

Download references

Acknowledgements

We are grateful to Prof. Axel Timmermann and the group involved in the 3 Ma experiment. We also acknowledge Prof. Fei-Fei Jin for valuable input. This study was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean Government (MSIT) (NRF−2018R1A5A1024958). S.-I.A. was supported by a Yonsei Fellowship funded by Lee Youn Jae. K.-S.Y. was supported by the Institute for Basic Science (project code IBS-R028-D1). The CSEM1.2 simulations were conducted on the IBS supercomputer “Aleph”. We also acknowledge the support of KREONET.

Author information

Authors and Affiliations

Irreversible Climate Change Research Center, Yonsei University, Seoul, Republic of Korea
Gagan Mandal, Soon-Il An, Chao Liu & Seungmok Paik
Department of Atmospheric Sciences, Yonsei University, Seoul, Republic of Korea
Soon-Il An
Division of Environmental Science and Engineering, Pohang University of Science and Technology, Pohang, Republic of Korea
Soon-Il An
School of Earth and Environmental Sciences, Seoul National University, Seoul, Republic of Korea
Jae-Heung Park
Center for Climate Physics, Institute for Basic Science (IBS), Busan, Republic of Korea
Kyung-Sook Yun
Pusan National University, Busan, Republic of Korea
Kyung-Sook Yun

Authors

Gagan Mandal
View author publications
You can also search for this author inPubMed Google Scholar
Soon-Il An
View author publications
You can also search for this author inPubMed Google Scholar
Jae-Heung Park
View author publications
You can also search for this author inPubMed Google Scholar
Kyung-Sook Yun
View author publications
You can also search for this author inPubMed Google Scholar
Chao Liu
View author publications
You can also search for this author inPubMed Google Scholar
Seungmok Paik
View author publications
You can also search for this author inPubMed Google Scholar

Contributions

G.M. and S.-I.A. conceived this study, and S.-I.A. supervised this research. G.M. carried out the data analysis, investigation, visualization and wrote the original manuscript draft. K.-S.Y. performed the 3 Ma experiment. G.M., S.-I.A., and J.-H.P. revised the manuscript. G.M., S.-I.A., J.-H.P., K.-S.Y., C.L., and S.P. contributed to the discussion and revision of the manuscript.

Corresponding author

Correspondence to Soon-Il An.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature Communications thanks Christian Stepanek, who co-reviewed with Fernanda Matos, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary Information

Transparent Peer Review file

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Mandal, G., An, SI., Park, JH. et al. Northern Hemisphere sea ice variability in a transient CGCM simulation of the past 2.6 Ma. Nat Commun 16, 39 (2025). https://doi.org/10.1038/s41467-024-55327-2

Download citation

Received: 02 May 2024
Accepted: 06 December 2024
Published: 02 January 2025
DOI: https://doi.org/10.1038/s41467-024-55327-2