Late Pleistocene 100-kyr glacial cycles paced by precession forcing of summer insolation

Climate variability over the past 800,000 years has long been described as being dominated by ~100-kyr glacial cycles, and researchers have debated whether these glacial cycles were driven by Earth’s orbital cycles of eccentricity, obliquity and precession. Some recent studies have suggested that these ~100-kyr glacial cycles are best characterized as groupings of two or three 41-kyr obliquity cycles; however, age uncertainties have made it difficult to distinguish whether the dramatic changes in ice-sheet size were more associated with 41-kyr obliquity or ~23-kyr precession cycles. We compare the impacts of obliquity and precession on glacial cycles using improved age estimates to analyse orbital phases during the onset of glacial terminations. Terminations are dated using a 640-kyr multiproxy stack of eight North Atlantic benthic δ18O records with well-constrained probabilistic age estimates derived from correlating North Atlantic ice-rafted debris to instances of abrupt Asian monsoon variability in high-resolution 230Th-dated speleothems. Rayleigh’s R statistics for the precession and obliquity phases of terminations demonstrate that, although both have statistically significant effects, the precession phase is more predictive of termination onset, particularly for the largest events. Thus, we conclude that Late Pleistocene ice sheets were sensitive to the precession forcing of Northern Hemisphere summer insolation intensity. Orbital precession played a more important role than obliquity during Late Pleistocene swings in ice-sheet extent, according to an analysis of benthic oxygen isotope records with precise age constraints.

Climate variability over the past 800,000 years has long been described as being dominated by ~100-kyr glacial cycles, and researchers have debated whether these glacial cycles were driven by Earth's orbital cycles of eccentricity, obliquity and precession. Some recent studies have suggested that these ~100-kyr glacial cycles are best characterized as groupings of two or three 41-kyr obliquity cycles; however, age uncertainties have made it difficult to distinguish whether the dramatic changes in ice-sheet size were more associated with 41-kyr obliquity or ~23-kyr precession cycles. We compare the impacts of obliquity and precession on glacial cycles using improved age estimates to analyse orbital phases during the onset of glacial terminations. Terminations are dated using a 640-kyr multiproxy stack of eight North Atlantic benthic δ 18 O records with well-constrained probabilistic age estimates derived from correlating North Atlantic ice-rafted debris to instances of abrupt Asian monsoon variability in high-resolution 230 Th-dated speleothems. Rayleigh's R statistics for the precession and obliquity phases of terminations demonstrate that, although both have statistically significant effects, the precession phase is more predictive of termination onset, particularly for the largest events. Thus, we conclude that Late Pleistocene ice sheets were sensitive to the precession forcing of Northern Hemisphere summer insolation intensity.
Nearly a century ago, Milutin Milanković proposed that glacial cycles were caused by changes in Northern Hemisphere insolation at high latitudes during the summer months, as affected by orbital cycles in precession, obliquity and eccentricity 1 . Although statistically significant climate impacts have been identified for all three orbital cycles [2][3][4][5] , the cause of the predominantly 100-kyr cyclicity in Late Pleistocene ice volume remains unresolved. Although orbital eccentricity has an ~100-kyr cycle, the insolation forcing over ice sheets contains almost no 100-kyr power 5,6 . Strong 100-kyr power in climate records emerged ~800,000 yr ago during the mid-Pleistocene transition (MPT), before which glacial cyclicity was dominated by the 41-kyr obliquity cycle. Thus, debate has focused on whether Late Pleistocene ~100-kyr glacial cycles are more inherently associated with 41-kyr obliquity cycles 2,3,7,8 or amplitude modulation of the precession index by the 100-kyr eccentricity cycle 4,9,10 .
Studies favouring eccentricity and precession forcing note that glacial maxima occur during weak precession cycles (associated with eccentricity minima) and propose that these extended times of low summer insolation intensity provide an opportunity for large, unstable '100-kyr' ice sheets to accumulate 4,[9][10][11] . Thus, the eccentricity modulation hypothesis posits that the amplitude modulation of precession by the ~100-kyr eccentricity cycle sets the timing of weak insolation variability that allows enough time for large volumes of ice to accumulate. Larger Northern Hemisphere ice sheets increase the amount of ice at lower latitudes, where precession has relatively more influence on summer insolation, and decrease the elevation of the ice margin by glacio-isostatic adjustment 10 , which may explain why the increase in sensitivity to precession and, thus, the 100-kyr power of climate variability coincided with an increase in Northern Hemisphere ice volume Article https://doi.org/10.1038/s41561-023-01235-x age model for the past 640 kyr (Fig. 2). Asian monsoon variability in speleothems has been established as nearly synchronous with northern high-latitude millennial-scale events recorded in ice cores, attributed to fast atmospheric propagation of North Atlantic climate changes to monsoon regions 16,17 . Throughout the Late Pleistocene, low-latitude monsoons are shown to be highly sensitive to North Atlantic meltwater forcing 17 .
Chinese speleothem δ 18 O records exhibit responses to both abrupt meltwater-induced atmospheric perturbations and precession-driven insolation change 17 . However, when correlating specific IRD and AAMV events, we focused on abrupt isotopic shifts rather than the slower, large-scale orbital variability. We primarily identified AAMV events using the detrended record from ref. 16 that removed orbital-scale power from the speleothem record and evaluated their similarity to millennial-scale IRD events; thus, the orbital signals in the speleothem records should not introduce bias to the revised marine sediment core age models (Extended Data Fig. 2). Furthermore, rather than indirectly identifying the termination ages by the initiation of weak monsoon intervals 8 , we define the onset of a glacial termination using the rate of benthic δ 18 O change in cores near the source of glacial meltwater.
Bayesian age modelling software, BIGMACS 18 , was used to create multiproxy age models for each of the eight North Atlantic cores and a stack of their benthic δ 18 O records. BIGMACS probabilistically combines information provided by δ 18 O alignment, absolute age estimates (for example, IRD-AAMV tie points), and a model of sediment rate variability for marine sediment cores. The output of this software is a continuous regional benthic δ 18 O stack of eight North Atlantic cores [19][20][21][22][23][24][25][26] , seven of which have ages constrained by IRD tie points (Fig. 3). during the MPT from 1.25-0.7 Myr ago 11 . Alternatively, the Laurentide ice sheet may have been sensitive to precession throughout the Pleistocene, but its response before the MPT was masked by anti-phased precession responses in the Southern Hemisphere 12 .
In contrast, the obliquity skipping hypothesis 2,8 argues that the MPT is most easily understood if obliquity serves as the main driver of ice volume variability throughout the Pleistocene. This hypothesis posits that, during the Late Pleistocene, glacial terminations occurred every two or three obliquity cycles, with an average duration of 100 kyr. However, it is important to note that Late Pleistocene palaeoclimate records demonstrate climate sensitivity to both obliquity and precession, and some models suggest that both obliquity and precession forcing are needed to accurately reproduce the timing of Late Pleistocene terminations 7,10,11,13,14 .

Challenges in resolving obliquity versus precession pacing
Identifying the orbital configurations associated with the onset of deglaciation would resolve the relative impacts of obliquity and precession in pacing 100-kyr glacial cycles 2,8,13 and help characterize the sensitivity of glacial climate to different patterns of insolation forcing. Late Pleistocene glacial terminations are often identified in ocean sediment cores as large, rapid decreases in benthic δ 18 O records, which are measured from the carbonate shells of benthic foraminifera and serve as a proxy for ice volume and correlative changes in deep-sea temperatures.
The biggest impediment to resolving the orbital configurations associated with the onset of deglaciation is age uncertainty, because orbital tuning cannot be used for age models that test orbital forcing hypotheses 2,4 . Beyond the radiocarbon dating limit of ~50 kyr ago, untuned age estimates for sediment core records typically rely on sparse magnetic reversals 15 . Assuming a constant mean sedimentation rate between magnetic reversals generates large age uncertainties, with an average standard deviation of 9 kyr (approximately half a precession cycle) for glacial terminations of the past 800 kyr (refs. 2,4).
Analysis using such an age model suggested that Late Pleistocene glacial terminations were clustered with respect to the obliquity phase, but not precession, leading to the development of the obliquity skipping hypothesis 2 . However, accurately evaluating the precession phase of terminations requires smaller uncertainties, because precession is a shorter cycle. Studies that analyse the amplitude of insolation forcing associated with terminations, instead of their relative timing, find that obliquity and precession contribute approximately equally to the forcing associated with glacial terminations 3,7 .
Recently, ages for two terminations (TX and TXII) during the MPT were estimated using tie points between marine sediment cores and a 230 Th-dated speleothem record 8 . Similar obliquity phases for these and subsequent terminations were interpreted as evidence that Late Pleistocene terminations are predominantly obliquity paced, but that study did not statistically evaluate precession pacing. Another study found that precession did pace Northern Hemisphere ice ablation over the past 1.7 Myr (although not necessarily global ice volume) and that the phase of ice-sheet responses to obliquity has changed across the MPT 11 . Because the MPT is associated with changes in the duration of glacial cycles, the orbital sensitivities of 100-kyr glacial cycles are best characterized without including MPT terminations in the analysis.

A more precise age model
Here we use speleothem-based age models, similar to Bajo et al. 8 , to provide improved termination age estimates suitable for evaluating the precession phase of Late Pleistocene terminations. Specifically, we correlate ice-rafted debris (IRD) layers in North Atlantic cores to abrupt Asian monsoon variability (AAMV) recorded in Chinese speleothem δ 18 O records 16 (Fig. 1, Extended Data Fig. 1 16 . Grey shading highlights millennial-scale events that could potentially be used as tie points. Vertical dashed lines mark tie points that were used to generate the sediment core age models, usually at the start and/or end of millennial-scale events. Not all cores recorded the same set of events.
Article https://doi.org/10.1038/s41561-023-01235-x Our North Atlantic benthic δ 18 O stack differs slightly from the global LR04 benthic δ 18 O stack 27 because our stack only contains cores from the North Atlantic and its ages are derived from IRD-AAMV tie points instead of orbital tuning. Because we are interested in the orbital configurations associated with termination onset, we identify the age of each termination (Table 1) based on when the stack's rate of change first exceeds a specified threshold (Methods), similar to Huybers and Wunsch 2 . This method identifies nine events, which we categorize as either full terminations (transitioning from full glacial to full interglacial conditions) or partial terminations (Methods).  Supplementary Table 3. The shaded grey region is the 95% confidence interval for this age model (n = 1,000). b, The relative age uncertainty throughout the stack's composite age model (shaded grey).

Termination ages and orbital phases
Our termination age estimates generally differ from recently published ages 8,13 by only a few kiloyears (Table 1) and with much less uncertainty than ages based on constant sedimentation rates 2,4 . Our age uncertainty estimates represent the combined effects of uncertainties in tie-point identification, speleothem ages, benthic δ 18 O alignment during stack construction and stack δ 18 O values (Methods). The 95% confidence interval (CI) half-width for all termination ages averages 3.6 kyr (or 3.1 kyr if TVIIa is excluded). Although the shifts in individual termination ages from previous studies are modest, their precession phases change more than obliquity phases (Methods) due to precession's shorter cycle length and, thus, higher sensitivity to age uncertainty. The orbital phases of all nine identified terminations ( Fig. 4) were used to calculate a Rayleigh's R statistic for obliquity (R ob ) and precession (R pr ), which quantifies how tightly clustered the phases are (Methods) 2 . We test the null hypotheses that the nine glacial terminations occur independently of obliquity or precession using the same stochastic model as Huybers and Wunsch (Methods) 2 . The critical R value identified using this model is a function of both the number of terminations and the orbital cycle examined. For nine terminations, the critical R values for obliquity and precession are 0.56 and 0.57, respectively (Methods).
The obliquity phases of terminations in our stack produce R ob = 0.65, comparable to previous studies 2, 8 , and the precession phases yield R pr = 0.76, notably higher than previously calculated. Whereas Huybers and Wunsch 2 failed to reject the null hypothesis for precession pacing, we find that R ob and R pr are both statistically significant, exceeding the critical R values for obliquity and precession, respectively (Methods). Thus, the null hypothesis may be rejected for both orbital cycles, signifying that both obliquity and precession both affect the timing of terminations.
To evaluate the effect of age uncertainty on these R values, we repeated the above phase calculations for termination ages identified in Monte Carlo samples of the probabilistic stack (Methods).
The Rayleigh's R samples for precession have a 95% CI range of 0.39-0.81, and the obliquity 95% CI is 0.53-0.75 (Fig. 4). The higher sensitivity of R pr to age uncertainty is evident from its wider CI compared to R ob . Of the 100,000 R samples, 94% of R ob samples and 69% of R pr samples exceed the critical value. However, the termination age estimates derived from the median of our probabilistic stack, which are more reliable than any individual sample, produce an R pr value of 0.76, which is higher than the upper limit of the 95% CI for R ob .
We further categorize each termination event as either a full or partial termination. In agreement with a previous study 7 , we define full terminations as those that transition from a fully glacial state to a fully interglacial state. Partial terminations either do not start from a glacial state (TIIIa and TVIIa) or do not reach interglacial conditions (TVI). This distinction allows for the evaluation of how orbital phases differ for large-amplitude terminations compared to weaker, partial terminations.
When including only the six full terminations, R pr increases to 0.83 and R ob decreases slightly to 0.61. Thus, the climate feedbacks that produce the complete loss of large Laurentide and Eurasian ice sheets appear to be particularly sensitive to precession. Interactions between obliquity and precession, such as the combined or competing influences of both orbital cycles, may explain why the precession phases differ between terminations. The variable ice-sheet size at each glacial maximum may also affect the precession phase of termination onset based on the hypothesis of Barker and colleagues that the

Precession sensitivity of Late Pleistocene ice sheets
In summary, our multiproxy North Atlantic stack with well-constrained termination ages shows that, although obliquity and precession both play statistically significant roles in termination timing, precession appears more important than obliquity in predicting the onset of Late Pleistocene glacial terminations. Precession phase appears particularly important for larger terminations. Furthermore, the disproportionate effects of age uncertainty on precession phase may yield a continued underestimation of precession's relative importance compared to obliquity. Additionally, we find no evidence that Late Pleistocene glacial cycles have typical durations of ~80 and ~120 kyr as predicted by the obliquity skipping hypothesis. However, we emphasize that obliquity also plays a fundamental role in the 100-kyr world as demonstrated by the statistical significance of Rayleigh's R for obliquity.
Collectively, our findings support the eccentricity modulation hypothesis that precession forcing of summer insolation intensity plays a crucial role in pacing 100-kyr glacial cycles 4,[9][10][11][12][13]28,29 and validates descriptions of the Late Pleistocene as a '100-kyr world'. Models of the Laurentide ice sheet demonstrate that, as its volume increases and its southern extent reaches lower latitudes, it becomes more sensitive to precession summer insolation changes due to differences in the seasonal and latitudinal forcing between obliquity and precession and the effects of glacial isostatic adjustment 10 . Thus, the increased sensitivity of large Late Pleistocene ice sheets to precession-driven insolation intensity combined with the eccentricity modulation of precession amplitude 9 explains the ~100-kyr glacial cycles in the Late Pleistocene.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-023-01235-x.

Identifying abrupt Asian monsoon variability
The timings of AAMV were identified using the δ 18 O record from a composite of Chinese speleothems 16 that are among the best dated palaeoclimate records extending back to 640 kyr BP, with an average δ 18 O resolution of ~120 yr. We calculated median age estimates and age uncertainties for the speleothems that make up the composite record using the speleothem age modelling software StalAge 30 . StalAge uses individual U-series ages, their corresponding uncertainties, and stratigraphic information to generate age constraints. The algorithm uses a Monte Carlo simulation to produce an age model with 95% confidence limits that include the uncertainties for the age data, identifies major and minor outliers, and removes age inversions. There are no adjustable parameters. We used StalAge to create individual age models (Extended Data Fig. 3) for each speleothem in the composite record, which come from the Dongge (D-3, D-4 and D-8) [31][32][33] , Hulu (H) 31,34 , Linzhu (LZ15) 35 and Sanbao caves (SB-3, SB-11, SB-12, SB-14, SB-32 and SB-58) 35,36 . After identifying the start and end of each AAMV event in the composite record, we determined the age of the event and its uncertainty based on the StalAge results for the individual speleothem in which the event occurred. When multiple speleothems recorded the event, we used the age and uncertainty from the record with smaller age uncertainty and higher δ 18 O resolution for that event.
Oxygen isotopes of Chinese speleothems reflect both changes in the strength of the Asian monsoon as well as changes in the precipitation source regions. For this reason, we make a distinction between what we define as abrupt Asian monsoon variability AAMV and weak monsoon intervals (WMI) 16 . The onset of AAMV was determined by the identification of abrupt δ 18 O increases in the composite speleothem record (Fig. 1). More specifically, the beginning of an AAMV was generally identified as the first point after which there is a noticeable increase in the rate of change of δ 18 O in the composite age model. Because the ends of each instance of AAMV tend to show more gradual changes in δ 18 O value, they are often not as distinct as initiations. The end of each of these intervals is defined by the point at which the δ 18 O value is similar to the value at which the AAMV began. The detrended speleothem record removes the influence of changes in insolation due to orbital forcing 16 . Although we used both the detrended and non-detrended composite record when deciding which δ 18 O changes may have been rapid enough to constitute an AAMV, we placed greater emphasis on the detrended record for the identification of potential tie points (Fig. 1).
The earliest and latest possible dates for both the start and end of each instance of AAMV were also identified, providing error bars for our tie-point identification. The uncertainty range for the start of each AAMV begins with the oldest δ 18 O data point that could reasonably define its onset, in most cases this is when the δ 18 O values first begin to increase. The latest age for the start of the AAMV was identified as the point where the rate of change of δ 18 O begins to decrease but the δ 18 O value has not yet peaked. The range of possible timings for the end of AAMV is defined similarly, but the confidence intervals for the end of the AAMV were more difficult to identify because the δ 18 O values tend to decrease more gradually.

Identifying North Atlantic ice-rafting events
North Atlantic IRD events were identified using previously published data from seven marine sediment cores that also have published benthic δ 18 O records [19][20][21][22][23][24][25] . These cores, located in or near the IRD belt in the North Atlantic Ocean (Extended Data Fig. 1 Fig. 4), but it was not used to generate the stack's age model due to a lack of IRD data.
The initiation of IRD events was determined by abrupt changes in IRD proxies ( Fig. 1 and Supplementary Table 1). Most records measured IRD grains per gram 19,21,22,25 , the sample ratio of IRD grains (250 µm to 2 mm in size) to planktonic foraminifera 20 , or the ratios of dolomite and quartz to calcite 24 . For site U1308, the indirect IRD proxies used are bulk δ 18 O (a proxy for biogenic versus terrestrial in the carbonate source), the ratios of calcium and silicon to strontium, bulk density and magnetic susceptibility 23 . Because core U1308 has multiple indirect proxies for ice-rafting, all proxies were compared simultaneously when identifying IRD events. IRD peaks in U1308 were only identified if a response was observed in more than one IRD proxy. The highest-resolution proxy was typically used to determine the exact start/end date for ice-rafting, unless the peak was more clearly resolved in a lower-resolution proxy (for example, due to a low signal-to-noise ratio in the high-resolution proxy). The core depths for the beginning and end of each IRD event were defined similarly to AAMV, including identification uncertainty.

Tie points and age uncertainty
Age-model tie points were created when IRD events could be linked to a corresponding AAMV event in the speleothem δ 18 O record. Not every IRD event had a corresponding AAMV and vice versa. Potential shifts between the original sediment core age models and revised age models were evaluated based on whether they were within the age uncertainty windows for previous benthic δ 18 O age estimates and based on similarity between the size and grouping of identified IRD and AAMV events. We expect most IRD-AAMV tie-point ages to differ from ages in the LR04 benthic δ 18 O stack by less than ~4 kyr, based on the age uncertainty estimates for the stack 27,34 . We also do not expect a large IRD event to correlate to a small shift in speleothem δ 18 O or vice versa. It was often important to consider all three criteria: age shift, size and grouping. If there were two instances of AAMV occurring near each other, and two IRD events of similar sizes, it was assumed that those two pairs of events should be correlated even if the shift in age was more than 4 kyr. If a single IRD event could potentially correlate to two different AAMVs, a tie point was recorded for each of them, with one serving as the primary tie point that would be used in calculations, and the other as an alternative tie point. Alternative tie points were incorporated into the age uncertainty for the tie point that was used in the stack. For example, when approaching the older part of the records where the uncertainty is higher, there are three instances of AAMVs between ~625 and 635 kyr. In core U1308, there are two possible IRD events that could reasonably be tied to any combination of the three AAMV events. This broad window of possible tie points was reflected in larger identification uncertainties.The next section describes how these tie points were used during construction of the North Atlantic benthic δ 18 O stack and its age model.
Multiple sources of uncertainty were combined to reflect accurately the total uncertainty associated with each IRD-AAMV tie point ( Supplementary Tables 1-3). The earliest and latest start and end ages for each identified AAMV event in the speleothem record from the StalAge age models were assumed to represent a 95% CI for the age of the event, which we described using a Gaussian probability density function (PDF) with a standard deviation equal to one-quarter of the range of the 95% CI. To calculate the combined age uncertainty for each AAMV event, we also account for the estimated standard deviation of uncertainty for AAMV identification. We assume the two sources of uncertainty (event identification and age model) are independent, and thus calculate the standard deviation of the total AAMV age uncertainty as the square root of the sum of squared standard deviations (Supplementary Tables 2 and 3).
The age uncertainty for IRD events was calculated in a similar way (Supplementary Table 1). The difference between the lowest and highest depth for each identified instance of IRD was assumed to represent the 95% CI for the depth of that tie point in the sediment core. Again, we assume a Gaussian uncertainty distribution for the tie-point

Nature Geoscience
Article https://doi.org/10.1038/s41561-023-01235-x identification and approximate the uncertainty standard deviation as one-quarter the 95% CI width. Finally, the standard deviation for depth uncertainty is converted to an age uncertainty by dividing by the core's estimated sedimentation rate at the depth of the tie point, which approximates the one-standard-deviation age uncertainty for each tie point in the IRD records. These sedimentation rate estimates were based on the core's median age model generated during creation of an initial stack. We do not incorporate the impacts of sedimentation rate uncertainty in this step because the tie-point uncertainty must be defined before the final age model can be generated. However, sedimentation rate uncertainty should have a negligible impact over the small depth ranges considered for each IRD event.
The IRD and AAMV tie-point identification uncertainties were assumed independent and were combined using the square root of the sum of their squared standard deviation. This yields a standard deviation for the full tie-point age uncertainty, which is listed in the 'additional_ages.txt' files for each core input to BIGMACS to describe the tie points' Gaussian age distributions.
Creating sediment core age models and a benthic δ 18 O regional stack Revised age-depth models for all cores and a stack of their δ 18 O records were created using the recently developed Bayesian age modelling software BIGMACS 18 . BIGMACS, which stands for 'Bayesian inference Gaussian process regression and multiproxy alignment for continuous stacks', creates age-depth models with quantified uncertainties for multiproxy alignment. Multiproxy alignment probabilistically combines information provided by δ 18 O alignment, absolute age estimates (in this case, IRD-AAMV tie points), and a model of sediment rate variability for marine sediment cores. The age uncertainties for the IRD and speleothem records were combined for the total tie-point uncertainty (as described in the previous section) and input to BIGMACS as a standard deviation for normally distributed age uncertainty at the start and/or end depth for each IRD event.
BIGMACS uses Gaussian process regression to model a benthic δ 18 O stack, and iteratively optimizes the alignment of all input cores to updated versions of the stack until converging to a stable, locally optimal solution. The first alignment is done to an initial stack, in this case the LR04 stack 27 . After the first alignment iteration, the software continues to optimize the parameters and estimated stack values to improve alignments until convergence. For example, BIGMACS estimates and applies average shift and scale factors for the δ 18 O values of each core during stack construction (Extended Data Fig. 4 and Supplementary Table 4). The software outputs a regional benthic δ 18 O stack of the North Atlantic cores with ages constrained by our tie points. Although the Gaussian process regression creates a continuous time series for the mean and standard deviation of benthic δ 18 O, we sampled and analysed the stack at a resolution of 0.1 kyr. The algorithm also outputs 100 different samples of the benthic δ 18 O stack. Additionally, sample age models are generated for each individual core with corresponding 95% CIs that reflect uncertainty in the alignment of the core's δ 18 O record to the stack.
To generate uncertainty estimates for the absolute age of the stack, we constructed an age model for a single core (IODP U1308) using only IRD-AAMV tie points. Because IODP U1308 is both continuous and consistently the highest-resolution record during the time interval examined, IRD-AAMV tie points from all other cores were mapped onto the composite depth scale of core U1308 using the median age models from the previous stack construction step (Supplementary  Table 3). After mapping all IRD-AAMV tie points to the U1308 depth scale, the age uncertainty for each tie point was increased to include the alignment uncertainty between the original core and site U1308. The output of this age-only BIGMACS run is 1,000 sample age models for the stack constrained only by the IRD-AAMV tie-point age estimates (Fig. 2 and Extended Data Fig. 5). The median stack age model is used for our primary analysis of results, and the 1,000 sample age models are used for age uncertainty estimates.

Identifying termination ages
Glacial terminations, defined as the rapid loss of large continental ice sheets, produce abrupt decreases in seawater δ 18 O that are recorded by benthic δ 18 O (refs. 5,6,35). We identified termination onsets based on the rate of change of δ 18 O for our median stack using a slope threshold of −0.02‰ per century (Extended Data Table 2), selected to agree with our analysis of sample stacks. The uncertainty in the identification of termination onsets is analysed using the slopes of the stack samples measured over a 1-kyr moving window, which still provides a slope estimate every 0.1 kyr. We again use a slope threshold of −0.02‰, which corresponds to the 95th percentile of the sample slopes (Extended Data Figs. 6 and 7). Because there are several instances where the smoothed slope exceeds our threshold throughout the record, we define identification windows of ±5 kyr from the termination onset ages identified in the median stack.

Calculating orbital phases of termination onset
The phase of obliquity and precession forcing for each glacial termination was calculated using the same technique as ref. 2, which is linearly based on the fraction of time between the nearest minimum and maximum (that is, a half cycle) of orbital obliquity and precession, respectively. We calculate phases using the half cycle, because not all obliquity or precession cycles are equal in length, and the time between maxima and minima is not always symmetric. The orbital phases of the nine glacial terminations were plotted on phase wheels for obliquity and precession ( Fig. 4 and Extended Data Table 3).
The null hypothesis for each cycle is that the timing of glacial terminations is independent of obliquity or precession forcing and terminations thus occur with equal probability for all phase values. To evaluate the extent to which termination phases vary for each orbital cycle, we use the Rayleigh's R statistic, which is defined as where N is the numbers of terminations and φ n is the orbital phase at the nth termination 2 . In words, this equation is describing the mean vector of the termination phases plotted on the obliquity or precession phase wheel. The critical R statistic represents the threshold at which the clustering of phases, or mean of the vectors, exceeds 95% of realizations for the null hypothesis and depends on the number of termination events identified within the fixed time interval. The critical R value is calculated using the null-hypothesis stochastic model of ref. 2, which consists of a biased random walk with a threshold to trigger glacial terminations set so that the average length of time between terminations is ~100 kyr. This simplified random walk ice-volume variability model is described by the formula where V t is the ice volume in 1-kyr time steps and µ t represents a random change in ice volume drawn independently from a normal distribution with mean µ = 1 and standard deviation σ = 2 (ref. 2). If V t exceeds a threshold T o = 90, a glacial termination will be triggered, thus resetting the ice volume to zero linearly over a 10-kyr interval. In this model, the duration between glacial cycles is set to have a relatively normal distribution of 100 ± 20 kyr. The critical R value is calculated from the precession or obliquity phases for the appropriate number of consecutive terminations, in this case, nine. A total of 100,000 Monte Carlo simulations of the null hypothesis stochastic model of Huybers and Wunsch 2 are run until each generates nine terminations. Evaluating the https://doi.org/10.1038/s41561-023-01235-x R values for obliquity and precession phases from these stochastically simulated terminations yields a critical value for Rayleigh's R of 0.56 for obliquity, and 0.57 for precession.

Uncertainty of termination ages and Rayleigh's R
We evaluated the effects of stack uncertainty on termination ages and the R statistics for obliquity and precession. The age modelling software BIGMACS provides samples of the stack's δ 18 O values and age model in proportion to their estimated probabilities. Specifically, we use 100 samples of δ 18 O time series from the stack's Gaussian process regression combined with 1,000 age model samples consistent with the tie-point constraints. The nine termination events are first identified for each of the 100 stack samples based on a rate-of-change threshold of −0.02‰ per century (as described in the section Identifying termination ages). These 100 samples of each termination event are then interpolated to depths on the composite depth scale of core U1308. The resulting 100 termination depths each have 1,000 age samples, yielding a total of 100,000 age samples for each termination (Extended Data Figs. 6 and 7). Similar to the method used by Khider and colleagues, these samples provide 95% CIs for the ages of termination onset 37 . Obliquity and precession phases for each set of nine terminations are used to generate 100,000 samples of the Rayleigh's R values (Fig. 4) and thus CIs for the R statistic as well.

Extended Data Fig. 2 | Orbital phases of identified AAMV events. (a, b)
Obliquity and precession phases for all AAMV events that were used as age constraints for the multiproxy North Atlantic stack. The wide spread of orbital phases for these AAMV events affirm that these tie points were unlikely to have biased our benthic δ 18 O age model due to orbitally forced changes in monsoon strength. (c, d) Obliquity phases for precession minima and maxima, respectively. Filled circles indicate the precession minima/maxima which occur a near termination onset (note that two precession maxima appear relatively close to TIIIa*), while open circles represent those which do not occur near a glacial termination. (e) Precession phases of obliquity maxima. Black R values for Extended Data Fig. 2c-e are calculated for all orbital maxima/minima that occur from 0-640 kyr, while the green and pink R values are calculated only for those which occur near terminations.