Simulating the glacial cycles from Marine Isotope Stage (MIS) 53 to 36
We first conduct a numerical experiment using a combination of realistic variations of astronomical forcings and a constant CO2 forcing (pCO2 = 230 ppm; a representative of the estimated pCO2; Fig. 1C) 33–38 to identify the role of combined astronomical forcings in driving the glacial cycles during the EP. The time series of the simulated ice volume is shown in Fig. 2 (black line). The most striking feature of the result is that, even with constant pCO2, the 41-kyr glacial cycles of the NH ice volume are clearly simulated between MIS 53 and 36. The result indicates that the response of the ice-sheet dynamics to the variations in astronomical forcings in the NH high latitudes primarily drives the 41-kyr glacial cycles.
Our simulation exhibits the timing of each deglaciation and the shape of each glacial cycle that are very similar to the δ18O signature during this period (Fig. 2F), despite the δ18O signature reflecting not only the ice volume but also deep-ocean temperature 3,8. For example, the lengths of interglacial periods from the longest (e.g., MIS 49, 47, and 37) to the shortest (MIS 41) are consistent with those in the δ18O record (Fig. 2F, Supplementary Fig. 2C and 2D) (See Material and Methods for the definition). The power spectrum of the simulated ice volume shows a strong 41-kyr signal with weaker 23- and 19-kyr signals (Fig. 2G; black line), which is consistent with the high-resolution data from the North Atlantic (U1308 core site; purple line in Fig. 2F; Supplementary Fig. 2D) 8. The simulated deglaciations occur with a relatively short interval, i.e., less than 15 kyr (Supplementary Fig. 2E), which is also consistent with the observed trend for the rapid deglaciation before the MPT 26–28.
To compare our results with reconstructions of sea level, we add estimates of the variations in the volume of the SH ice sheet, namely the Antarctic ice sheet, by Raymo et al. 21 and de Boer et al. 39 (Red and orange lines, respectively, in Fig. 3C) to those of our NH ice sheet (black line in Fig. 3D). The work by Raymo et al. 21 proposed an offset between the NH and SH ice volumes, with a large variation in SH ice volume of up to ~ 30 m at MIS 48 (Fig. 3C red line). On the other hand, recent simulations of the SH ice sheet using a sophisticated ice-sheet model suggest that the amplitude of the variation of SH ice sheet was at most ~ 15 m sea level equivalent during the deglaciations from MIS 50 and 48 39–41. Despite the non-negligible amplitude of SH ice-sheet variation, the amplitude of sea-level variation is almost as large as that of the NH ice sheet alone (Fig. 3C, red and blue lines). The amplitude of the simulated sea-level change ranges between ~ 50–95 m, which is consistent with the range shown in the reconstructions (~ 60–90 m; Fig. 2A green and blue lines) 6,7. These results reinforce that the 41-kyr glacial cycles are primarily caused by the variation in the NH ice sheet owing to the variations in astronomical forcings. Considering the variation of the SH ice sheet simulated by these studies, the power spectrum of the variation in the modeled global sea level shows a clear dominance at 41-kyr periodicity even with a precession signal weaker than that of the NH ice-volume alone, because of the anti-phased precession signal in the SH (Fig. 3B and 3D) 21,22.
Effect Of Orbital-scale Variability Of Co Forcing
To elucidate the role of orbital-scale variability in pCO2 (i.e. the carbon cycle) in establishing the 41-kyr glacial cycles, we conduct further experiments using the published pCO2 variations (Fig. 2E, red and blue lines) 35,36. We use two reconstructions based on an inverse modeling approach from the δ18O records 35 and a method based on the δ13C gradient between the deep Pacific and intermediate North Atlantic 36. The results with these non-constant pCO2 reconstructions exhibit little difference from the simulation with constant pCO2 (230 ppm) in terms of the shape and the amplitude of the glacial cycles. The timing of the deglaciation and the dominance of the 41-kyr periodicity are also almost identical (Fig. 2G). Thus, it could be inferred that the establishment of the 41-kyr periodicity, shape, and amplitude of the glacial cycles at ~ 1.6–1.2 Ma are determined primarily by the variations in astronomical forcing (Fig. 2A–D and 2G). Below, we further analyze the model results with constant pCO2 (230 ppm) to understand the fundamental role of astronomical forcings in shaping the 41-kyr glacial cycles.
Characteristic Shape Of The 41-kyr Glacial Cycles
To characterize the shape of the simulated glacial cycles, we summarize the time series of the ice volume calculated under constant pCO2 (230 ppm) (black line) and the lengths of the glacial and interglacial periods in Fig. 2H (blue and red lines, respectively). For this analysis, the boundaries between the interglacial and glacial periods are defined as the mid-point of the maximum and minimum ice volume for each deglaciation that is observed in the δ18O signature (see Methods). Results indicate that the lengths of the interglacial and glacial periods are inversely correlated. In other words, when an interglacial period is long, the following glacial period tends to be short, and vice versa. The difference in the lengths of interglacial and glacial periods among the 41-kyr glacial cycles are also shown in Supplementary Fig. 2. The shape of the glacial cycles appears “plateau-like” when an interglacial period is long (cycles from MIS 49, 47, 39, and 37), whereas it takes a “U-shaped” form when an interglacial period is short (cycles from MIS 43). The intermediate stages of the U-shaped and plateau-like cycles exhibit a “saw-toothed” form (cycles from MIS 51, 45, and 41). Different time series of glacial cycles could be caused by different characteristics of the astronomical forcing. We will discuss more of this later.
Understanding The 41-kyr Glacial Cycles Based On The Threshold Mechanism
To understand why the glacial cycles of the NH ice sheet terminate every 41 kyr, we compare the three consecutive cycles from MIS 47 to 41 with the steady-state response of an ice sheet to astronomical forcing (Fig. 4A and 4E; See also Supplementary Fig. 3). The steady-state response of an ice sheet to the surface temperature anomaly due to insolation and CO2 was derived in our previous study 12. The cyan and magenta lines represent the steady-state achieved if the model is initiated from a small and large NA ice volume, respectively, showing a hysteresis behavior within the approximate range of − 4 to 2 K (gray-hatched area). The ice sheet gains mass if the temperature anomaly falls below the lower branch of the hysteresis curve (cyan line) and it loses mass if the temperature anomaly rises above the upper branch (magenta line). These lower and upper branches correspond to the thresholds of ice-sheet growth and retreat, respectively. In order to understand the establishment of the 41-kyr glacial cycles, we compare the results of the transient variation before the MPT for the period from MIS 47 to 41 shown as trajectories to the steady-state response.
For the typical 41-kyr cycles with a sawtooth shape (black lines in Figs. 4A and 4D), ice-sheet retreat occurs once in every two precession cycles. At the first insolation maximum (~ 1380 ka) after the onset of the glacial cycle, the intensity of insolation is insufficient to exceed the threshold of termination (upper branch of the hysteresis curve, magenta line in Fig. 4A). This is because the effects of precession and obliquity forcings cancel each other (Figs. 4B and 4C), so that the ice sheet continues to grow. The second local maximum of insolation (~ 1360 ka) is strong because it is compounded by the contributions from both obliquity and precession. Besides the strong insolation, there is also sufficient time for the ice sheet to retreat. After crossing the upper threshold, the ice sheet starts to lose mass at ~ 1364 ka because of the larger ice-sheet size than at the first local insolation maximum, in addition to the strong insolation. The combined effect of ice-sheet size and insolation allows the temperature to exceed the threshold, and thus triggers the onset of deglaciation (black lines in Figs. 4A and 4D). Once deglaciation has commenced, feedback from the glacial isostatic rebound assists in complete deglaciation until an interglacial state is reached (Supplementary Fig. 4) 12. As a result, the deglaciation occurs once in two precession cycles through a threshold mechanism that results from the hysteresis response of the ice sheet, as in the LP.
For the glacial cycle from MIS 47 to 45, the interglacial period is longer than one precession cycle (red lines in Figs. 4A and 4D). The ice sheet does not develop at the first minimum in precession forcing, but it does develop rapidly at the second. In contrast, the interglacial period is short in the case of the cycle from MIS 43 to 41 (blue lines in Figs. 4A and 4D). The ice sheet starts to develop soon after the previous deglaciation at the first minimum in precession forcing, and the glacial period continues for more than one precession cycle. In all three cases, the deglaciation occurs once in two precession cycles through the threshold mechanism. The differences between the shapes of the 41-kyr glacial cycles (red and blue lines in Fig. 4D) reflect the phase relation between obliquity and precession.
Timing Of Deglaciation And Length Of Interglacials
To quantify the contribution of each astronomical forcing (obliquity and precession) that determines the timing of deglaciation in our simulation, we conduct a phase analysis of obliquity and precession forcings at the deglaciations (see Methods) 27. In Figs. 5A and 5B, the so-called “phase wheels” for obliquity and precession, respectively, are shown. These phase wheels summarize the phase angle of the astronomical forcing at each deglaciation, defined in Fig. 2H (midpoint of the maximum and minimum ice volume; see also Supplementary Fig. 5), in a unit circle (marked by diamonds in Figs. 5A and 5B) 27,42. The average phase angle of both precession and obliquity, represented by the direction of the thick black line, is around zero, showing that the deglaciation occurs mostly at precession minima and obliquity maxima. This result infers that both obliquity and precession forcings pace the deglaciations of the 41-kyr glacial cycles, consistent with the discussions on the 100-kyr glacial cycles during the LP 10,12,13,43. While the phase angles of precession are tightly concentrated around 0˚ for all deglaciations investigated in this study, i.e., from − 12° to + 23° (Fig. 5A), those of obliquity are spread over a wider range, i.e., from − 67° to + 44° (Fig. 5B), suggesting that precession is crucial in determining the precise timing of deglaciation. Moreover, the dispersion of the phase angles at the deglaciations can be represented by the Rayleigh’s R value, which is denoted by the length of the black lines in Figs. 5A and 5B (1 corresponds to the same phase angles at all deglaciations). In our simulation, the Rayleigh’s R values for our deglaciations are 0.98 for precession and 0.80 for obliquity, which quantitatively indicates that precession primarily controls the timing of the deglaciations of the 41-kyr glacial cycles during the EP.
The phase relation between obliquity and deglaciation appears to control the length of the following interglacials. In order to analyze the relationship between the length of interglacial periods and the phase angles of precession and obliquity, the beginning and end of all simulated interglacials relative to the variations in precession and obliquity are shown in colored lines in Figs. 5C and 5D, respectively. The onset of the interglacials occurs at maxima of precession forcing (Figs. 5A and 5C). The duration of the interglacial periods ranges from 10 to 30 kyr (i.e., 0.5–1.5 precession cycles). For MIS 45, 43, and 41, the interglacial periods end within one precession cycle ( ≲ 20 kyr), while, for MIS 51, 49, 47, 39 and 37, they continue beyond one precession cycle ( ≳ 20 kyr). As shown in Fig. 5D, the duration of an interglacial that starts earlier than the obliquity maximum is long (i.e., MIS 49, 47, 39 and 37), whereas it is short when the interglacial starts after the obliquity maximum (i.e., MIS 43 and 41). The end of interglacial periods occurs around the minima of the obliquity and precession forcings, while the timing of onset of the interglacials is determined strongly by the maxima of precession forcing (Fig. 5C). Thus, when the maxima of precession forcing occurs before (after) the obliquity maxima, the interglacials tend to be longer (shorter) than one precession cycle. In other words, the length of an interglacial period is controlled by the lead–lag relationship between the maximum of obliquity and precession forcings.
The relative phase angle of the obliquity and precession forcings during a deglaciation determines the shape of the subsequent glacial cycle. The phase angle of obliquity at maximum of precession forcing, ΦPRE–OBL, is denoted by the direction of cyan arrows in Figs. 5E–G). On the right-hand side of Figs. 5E–G, three typical glacial cycles (Plateau-like, saw-toothed, and U-shaped) are shown schematically as red lines, together with idealized obliquity (black line) and precession forcings (cyan line). If ΦPRE–OBL is close to zero, a saw-toothed glacial cycle is formed (cycles from MIS 51 and 45 in Fig. 5F). If ΦPRE–OBL is negative, a plateau-like glacial cycle is formed (cycles from MIS 49, 47, 39, and 37 in Fig. 5E), while if it is positive, a U-shaped glacial cycle is formed (cycles from MIS 43 in Fig. 5G). Note that, although we categorize MIS 41 as saw-toothed, ΦPRE–OBL at that time is positive, i.e., corresponding to a U-shaped cycle. For this specific case at MIS 41, rapid glaciation is not achieved owing to the small eccentricity (Fig. 2A).
Effect Of Amplitude Of Eccentricity And Obliquity
In order to understand the effect of the amplitude of obliquity and eccentricity (i.e., the amplitude of precession forcing) on the 41-kyr glacial cycles at ~ 1.6–1.2 Ma, we conduct additional sensitivity tests. We use a set of artificially modified astronomical forcings (amplitudes of eccentricity reduced to 60% and 20%, and amplitude of obliquity reduced to 60%) (Fig. 6A and 6D). When the amplitude of eccentricity is reduced, the 100-kyr periodicity of the ice volume intensifies (blue and red lines in Fig. 6A). For the case of 20% eccentricity amplitude (red line in Fig. 6A), only deglaciations from MIS 50, 46, and 40 are reproduced. Owing to the small eccentricity, other deglaciations are not realized and the ice sheet continues to develop. Conversely, with a decreased obliquity amplitude (60% as shown by the red line in Fig. 6D), the glacial cycles also show 100-kyr periodicity (red line in Fig. 6F).
Dependence on the mean levels of atmospheric pCO2
To clarify the relative importance of the astronomical forcings to the mean levels of pCO2 for the period ~ 1.6–1.2 Ma, we conduct a set of calculations with different constant values of atmospheric pCO2 (210, 230, 250 and 270 ppm) (Supplementary Fig. 6). The amplitude of the glacial cycle is smaller for the case with high pCO2. Even with different atmospheric pCO2 down to 210 ppm, the dominant 41-kyr periodicity of the glacial cycles during this period remains. The relative power of the ~ 20-kyr component remains almost unchanged.
Changes Of Ice-sheet Thickness And Extent
The ice sheet distribution for each glacial maximum for the experiment with constant pCO2 (230 ppm) is shown in Figs. 7A, 7B, and Supplementary Fig. 7. In most glacial maxima, the Laurentide and Cordilleran ice sheets are separated and the southern margin does not extend beyond the Great Lakes. In one case (MIS 48 in Fig. 7B), the Laurentide and Cordilleran ice sheets are merged and the glacial ice sheet covers a large area of the continent, including the Great Lakes. This glacial cycle (MIS 48) corresponds to a case categorized as a long interglacial and plateau-like glacial cycle. The ice sheet at the glacial maximum of MIS 48 is geographically extensive and thin, because it expands quickly within half of one precession cycle and retreats before the NA ice sheet reaches equilibrium (See Supplementary Movie 1). This ice sheet is much thinner than the Last Glacial Maximum (LGM) ice sheet (Figs. 7A and 7C, and Supplementary Fig. 8). We conjecture that the relative phase angle at a deglaciation that produces the plateau-like cycle together with high eccentricity is one of the conditions contributing to a thin and extensive North American ice sheet during the EP (Figs. 2A and 2D). The southern edge of the NA ice sheet at MIS 48 reaches the Mississippi drainage basin 18,44, as suggested by geological reconstruction 45,46. If a thick regolith layer had covered the North American continent before the MPT, the North American ice sheet might have extended over an even larger area of the continent 9, 45–48.