The role of Northern Hemisphere summer insolation for millennial-scale climate variability during the penultimate glacial

Previous glacial intervals were punctuated by abrupt climate transitions between cold (sta-dial) and warm (interstadial) conditions. Many mechanisms leading to stadial-interstadial variability have been hypothesized with ice volume being a commonly involved element. Here, we test to which extent insolation modulated stadial-interstadial oscillations occurred during the penultimate glacial. We present a replicated and precisely dated speleothem record covering the period between 200 and 130 ka before present from caves located in the European Alps known to be sensitive to millennial-scale variability. We show that the widely proposed relationship between sea level change and stadial-interstadial variability was additionally modulated by solar insolation during this time interval. We ﬁ nd that interstadials occurred preferentially near maxima of Northern Hemisphere summer insolation, even when sea level remained close to its minimum during peak glacial periods. We con ﬁ rm these observations with model simulations that accurately reproduce the frequency and duration of interstadials for given sea-level and insolation forcing. Our results imply that summer insolation played an important role in modulating the occurrence of stadial-interstadial oscillations and highlight the relevance of insolation in triggering abrupt climate changes.

M illennial-scale interstadial-stadial cycles were first identified in ice-core δ 18 O records from Greenland and named Dansgaard/Oeschger events 1,2 . Associated abrupt transitions from stadial to interstadial conditions were characterized by temperature increases of up to 16°C within a few decades 3 . Interstadial conditions typically lasted a few hundreds to several thousands of years. Although subsequently identified also in several other archives 4,5 especially in the Northern Hemisphere (NH), the cause of stadial-interstadial oscillations is still intensively debated. Proposed triggering mechanisms commonly involve variations in the strength of the Atlantic Meridional Overturning Circulation (AMOC), which is thought to react to meltwater discharge from large continental ice sheets [6][7][8][9] , CO 2 -induced changes in the North Atlantic tropical hydrology 10 , varying extent of sea ice possibly associated with modified wind stress due to changes in the elevation of the Laurentide ice sheet 11 , and ocean-sea ice interactions in the Nordic Seas 12,13 .
Data from Greenland ice cores 14 and North Atlantic sediment cores 15 as well as a synthetic Greenland isotope curve derived from Antarctic ice core data 16 consistently show that stadialinterstadial oscillations neither occurred during warm periods (interglacials) nor during the coldest periods of glacials, when sea level was close to its minimum [16][17][18] . Interstadials of the last glacial cycle occurred predominantly during intermediate climate conditions, e.g. Marine Isotope Stage (MIS) 3. NH summer insolation-here we use June-July-August (JJA) values at 65°Nis known to play a key role in dictating glacial-interglacial cyclicity 19,20 . Until now, however, it was difficult to scrutinize the role of the NH summer insolation on the occurrence of stadialinterstadial variability 21 , because changes in NH summer insolation 22 during most parts of the last glacial and especially during MIS 3 and 2 were relatively small. Records reaching further back in time commonly lack a precise chronology and/or are of insufficient resolution.
The goal of this study is to explore the influence of insolation on the occurrence of stadial-interstadial oscillations and the characteristics of interstadials during a glacial period that is characterized by larger variations in NH summer insolation than the last glacial period. We focus on speleothems covering the penultimate glacial because this combination offers sufficiently high precision to resolve millennial-scale events and to allow for an accurately dated chronological framework. We chose caves in the European Alps as this region has been shown to be highly sensitive to stadial-interstadial oscillations during the last glacial period 5,[23][24][25] . Oxygen isotope values in speleothems from this region have shown to respond sensitively to climate change with high (low) δ 18 O reflecting NH warm (cold) conditions during both orbital and millennial time scales 5,23-25 .

Results
We investigated four speleothems from the Melchsee-Frutt Cave system in central Switzerland (~1755 m above sea level, N 46°4 7.309, E 8°16.819, see Supplementary Notes) which cover the penultimate glacial. The speleothems formed between~202 and 133 ka and partly grew contemporaneously, providing a high degree of replication. In total, we performed 84 U-Th age determinations and measured more than 4200 stable isotope samples. A composite record based on the depth-age models of the four stalagmites was obtained using 'iscam' 26 (for details, see methods section "Age-depth modelling").
The stacked δ 18 O record shows a large decrease from −8 ‰ at the end of MIS 7 to −12 ‰ at the beginning of MIS 6 ( Fig. 1a). MIS 7 values are comparable to those of MIS 5 as recorded in other speleothems from the same cave site 27 . The transition into MIS 6 occurred abruptly at 187.6 ± 0.8 ka, while the penultimate termination (MIS 6 -5 transition) is not observed in an individual speleothem but can be inferred by comparing the specimen growing before and after the glacial termination. Our stacked record documents high-amplitude millennial-scale variability throughout MIS 6 ( Fig. 1, see methods section "Determination of interstadials" for details on the algorithm employed to detect interstadial onsets) that strongly resembles the iconic pattern of stadial-interstadial cycles recorded in Greenland ice cores for the last glacial interval. This includes the abrupt δ 18 O increase into interstadials, followed by a gradual decrease and a fast fall-back to stadial conditions. In total, we identified 16 interstadials in the penultimate glacial between~175 and 135 ka. This number is comparable to that of the last glacial period, for a comparable time window of about 40 ka. Note that our record lacks 7 ka of the penultimate glacial from~182-175 ka.
Millennial-scale variability was similarly observed in MIS 6 records in the North Atlantic realm [28][29][30][31] . A benthic δ 18 O record from the Iberian Margin 29 (MD01-2444) shows a stadialinterstadial succession comparable to our record ( Fig. 1), which reflects changes in the local deep water signature. The age model of this marine core was established by tuning to a Chinese speleothem record 32 . While the timing of the interstadials is not in full agreement with our record, clusters of interstadials at~170 and~150 ka are present in both records. Also, an alkenone-based sea-surface temperature record from the Iberian Margin 30 (MD95-2043) and pollen data from a sediment core from Lake Ioannina 31 in Greece show millennial-scale variability, with higher sea-surface temperatures and a larger contribution of arboreal pollen during interstadial conditions. Furthermore, periods of muted millennial-scale variability (i.e. around 185, 155 and 140 ka) as observed in our speleothem record agree well with planktonic and benthic δ 18 O records from the Iberian margin 29 (Fig. 1), as well as with pollen data from the same site (29, Fig. S9), Lake Ioannina pollen data 31 and an alkenone-based seasurface temperature reconstruction from another deep-sea core on the Iberian Margin 30 (Fig. 1). However, as the age models of these records are based only on few anchor points linked to orbital parameters or other climate records, it is difficult to compare them independently. Nevertheless, even without an independent chronology those records clearly indicate that interstadials were widely distributed events during the penultimate glacial, with pronounced changes in deep and surface water characteristics as well as in atmospheric circulation and associated impacts on vegetation.
Millennial-scale variability during MIS 6 was also reported from speleothems in China 32 and Peru 33 , where it is superimposed on a strong insolation-paced monsoon signal (Fig. 2b, c). In both regions, the change in δ 18 O is attributed to varying amounts of precipitation. In contrast to those records, the stable oxygen isotope records of the Alpine speleothems do not follow orbital cycles, as was also observed for the last glacial period 25 . A synthetic stable isotope record from Greenland based on the hydrogen isotope composition of an Antarctic ice core 16 also reveals millennial-scale variability, but shows distinct differences to our new speleothem record (Fig. 2d). Although there is some disagreement on the precise timing of interstadials, the different records generally agree regarding longer-term episodes with frequent occurrence of interstadials (e.g.,~175-159.2 ka and 151.1-146 ka), and more "quiet" episodes (186.5-178.6 ka, 159.2-151.1 ka and 146-138.4 ka). Some of those differences might be ascribed to dating uncertainties as well as to the method for estimating Greenland ice core δ 18 O from Antarctic ice core δD-reconstructed temperatures. However, clusters of interstadials are observed in all four records with approximately the same periods of occurrence between~175 and 161 ka as well as around 150 ka (Fig. 2).

Discussion
Two approximately 8 ka-long periods between~159.2 and 151.1 ka as well as between~146 and 138.4 ka with no interstadials can be observed in the Melchsee-Frutt record (Fig. 3a). The timing of these episodes in our record is consistent with the synthetic Greenland record 16 and the two speleothem records from the monsoon regions 32,33 with no stadial-interstadial transitions detected in the Greenland record and only very few in the two speleothem records during these periods. When combining information from our record and from the synthetic record a third~8 ka-long period without interstadials is apparent during the early part of the penultimate glacial between~186.5 and 178.6 ka (evidence from our record for the period between 186.5 and 182 ka). Within errors these three intervals coincide with maxima in global ice volume ( Fig. 3a) as derived from a high-resolution planktonic foraminifer δ 18 O record in the Red Sea 34 . This relationship is also consistent with other sea level reconstructions (Fig. S8). Furthermore, rapid successions of interstadials during the penultimate glacial are conspicuously common around maxima in NH summer insolation (at~175 and 151 ka). The same pattern is also seen in the last glacial period when an extended phase of frequent yet comparably short interstadials (i.e. MIS 3) followed a phase of low interstadial variability between~70 and 60 ka (Fig. 3b). Earlier studies [16][17][18] 17,18 ), necessary to precondition stadialinterstadial oscillations. This concept is not fully explaining the pattern of stadial-interstadial variability we observe in the precisely dated, high-resolution data for MIS 6 presented here. First, interstadials were markedly absent in the records shown in Fig. 2 during a phase of intermediate climate conditions between 186.5 and 178.6 ka when sea level was about −70 m, but contemporaneous with an interval of minimal NH summer insolation (Fig. 3). This period without millennial-scale variability came to an end as the NH summer insolation approached a local maximum. Second, between about 160 and 138 ka, sea level experienced a long lasting low-stand with only minor variations between −90 and −100 m. Nonetheless, between about 151 and 146 ka a series of interstadials occurred. This period of climate instability was accompanied by elevated NH summer insolation (Fig. 3). We note that, while NH summer insolation seems to support the occurrence of stadial-interstadial oscillation, it does not appear to be a necessary condition for their appearance as exemplified by some events around 187 and 163 ka. This suggests that different background climate configurations (e.g. ice sheet dimensions and sea level vs solar insolation) provide additional windows within which stadial-interstadial oscillations can occur.  16 . The speleothem records are on their own radiometric age scale, while the synthetic Greenland isotope record is tuned to the Chinese speleothems 16,32 .
These observations illustrate the importance of NH summer insolation in modulating the occurrence or suppression-and thus the frequency-of interstadials. As the NH summer insolation strongly affects the melting of sea ice and snow in the high latitudes, we hypothesise that higher NH summer insolation supports the initial warming at the onset of the interstadials, which lead to an initial reduction in sea ice 12,13,35,36 via the sea ice-albedo feedback, in turn resulting in a run-away reduction of sea ice. We argue that during peak glacial conditions, when sea level is at or close to its minimum (around 150 ka), an initial seaice retreat initiated by the triggering mechanism of stadialinterstadial oscillation is strongly enhanced under high NH summer insolation conditions, pushing the system over into the interstadial. In contrast, during intermediate climate conditions with respect to global sea level, e.g., around 183 ka, low NH summer insolation values lead to a less effective sea-ice albedo enhanced reduction in the sea ice extent, which might have prevented interstadials. This finding corroborates the importance of NH summer insolation. However, during other periods of low NH summer insolation, when sea level was higher than −70 m (e.g., around 115 ka or 85 ka), the triggering mechanism itself was sufficiently strong to induce full interstadial conditions without strong support of NH summer insolation.
Several other mechanisms possibly modulating millennial-scale variability as a result of changes in boreal summer insolation have been proposed. Variations in atmospheric CO 2 levels and freshwater fluxes from melting continental ice sheets into the North Atlantic, altered by NH summer insolation, have also likely contributed to the amplification of an initial warming, as shown in several modelling studies [8][9][10][11] . Margari et al. 29 proposed that European ice sheet melting and associated large seasonal meltwater discharge via the Fleuve Manche river into the North Atlantic may have led to AMOC disruptions, and extended stadial conditions around 155 ka. This could explain the absence of interstadials around this period in the Melchsee-Frutt Cave record and would be another possible mechanism which might have contributed to the pattern of millennial-scale variability we observe for MIS 6.
Furthermore, the role of boreal summer insolation on Asian monsoon strength during glacial maxima was previously investigated 32 . The occurrence of weak Asian monsoon events appears to be related to a slowdown of the AMOC, and thus cooling of the North Atlantic region, induced by the melting of the continental ice sheets due to increased solar insolation. However, an AMOC-related slowdown was often followed by enhanced stadial-interstadial variability 15,30 with strong interstadial conditions. Interstadials in Asia were accompanied by elevated amounts of rain. Both, more frequent interstadials and enhanced stadials can thus be related to increased boreal summer insolation, bringing together our explanation with Cheng et al. 32 .
There are no fully coupled climate models available that we could use to explore the combined influence of the climate background, as expressed by sea level, and summer insolation on the occurrence of interstadials during a full glacial-interglacial cycle. To investigate this influence nevertheless, we focus our analysis on a stylized model that comprises the key dynamical mechanisms involved in stadial-interstadial cyclicity 13 . This model has recently been shown to successfully reproduce observed features of the stadial-interstadial cycles of the last glacial, such as the sawtooth shape of interstadials, the shifted antiphase relationship between Greenland and Antarctic ice cores, and in particular the varying frequency of interstadial onsets in the course of the last glacial (for more details see methods section "Conceptual climate model experiments"). We use two model layouts, the model layout used by Boers et al. 13 and a slightly modified one. In the original model, interstadials are triggered when North Atlantic subsurface water temperatures surpass a constant threshold after which sea ice rapidly retreats. For the second model layout we replace the constant threshold by a dynamical one, which inversely follows NH summer insolation. For example, if NH summer insolation is high, the threshold for subsurface water temperature to trigger an interstadial is reduced accordingly. This can be justified as under a given small North Atlantic subsurface warming and associated sea-ice melting, high NH summer insolation is more likely to accelerate any initial perturbation leading into an interstadial by e.g., a more effective sea-ice albedo feedback. Under relatively low NH summer insolation conditions a similarly small subsurface warming might be not sufficient for the sea-ice-albedo-feedback to enhance the initial perturbation. A similar observation about a NH-summer insolation-dependent stadial-interstadial oscillation was provided by a modelling study 37 claiming that warmer boreal summers favour sea-ice reduction and hence increase the openwater area in the northern North Atlantic, which can contribute to enhanced ocean heat loss during winter. We compare the model results with the proxy data by means of the number of stadial-interstadial transitions per 10 ka periods for the last 200 ka. In addition, we also compare the number of years which are warmer than average stadial conditions per 10 ka period, which provides a measure for the duration of both interstadials and interglacials. To calculate the number of transitions and warm years from the climate proxy data, we use a combination of NGRIP δ 18 O data 14 (last 120 ka) and the δ 18 O data from the Melchsee-Frutt cave region (period between 200 and 120 ka; Fig. 1; this study and MIS 5 data 27 ). To fill the gap in the speleothem dataset we use the number and length of interstadials as obtained from the synthetic Greenland record 16 . The beginnings and ends of interstadials in the NGRIP ice core record are taken from Rasmussen et al. 14 . The events from the speleothem study are identified as explained in detail in the methods section "Determination of interstadials".
While both model versions depict the stadial-interstadial characteristics during the last glacial approximately equally well, the results differ largely during MIS 6. Good agreement with the proxy data for the period from around 200 ka to 120 ka is found for the number of stadial-interstadial transitions as well as for the number of warm years between the simulations obtained with a variable temperature threshold (Fig. 4). The model results, which are obtained from the constant threshold approach, show less convincing agreement with respect to the proxy data (Fig. S10). Especially, the number of interstadials during MIS 6 are better represented in the results obtained with a variable insolation threshold. The number of warm years also seems to fit better with the variable threshold approach (Fig. 4), both within the glacial, where the minima of the number of warm days are better represented, and interglacial, where the maxima are more realistically expressed than in the constant threshold model version (Fig. S10). Additionally, the performance during and at the end of the MIS 5 interglacial (especially between~110 and 80 ka) is more convincing in the variable threshold approach.
The comparison between the proxy data and model results thus supports the idea that the NH summer insolation played a more important role in triggering stadial-interstadial oscillation than previously thought. When summer insolation is low and cold conditions prevail, the likelihood of a stadial-interstadial oscillation is small, which is reflected in the muted or reduced occurrence of interstadials according to the proxy data during these intervals. This highlights the fact that neither climate background conditions nor NH summer insolation alone can fully explain these interstadial-free phases. We conclude that millennial-scale climate variability during glacial periods although likely directly triggered by internal feedbacks of the ocean-climate-cryosphere system, needs the right combination of background climate (i.e., sea level) and summer insolation to occur. We encourage to investigate the millennial-scale variability under various combinations of background climate and boreal summer insolation in more detail using climate models. Our results suggest that northern hemisphere summer insolation is not only a pace maker for glacial-interglacial transitions, but also plays a key role for millennial time scale climate variability, directly modifying the occurrence of stadial-interstadial oscillations.

Methods
Analytical methods. All stalagmites were analysed for their U and Th concentrations and respective isotopic ratios. Sub-samples of~0.2 g were hand-drilled along growth layers. In total 39 samples were measured for stalagmite M39_1_62; 23 for M39_1_61; 16 for M39_1_61B and 6 for M6_1_14. Isotopic measurements were conducted using a MC-ICP-MS (ThermoFisher Neptune Plus) at Heidelberg University, using the analytical protocols of Wefing et al. 38 .
Chemical preparation followed sample dissolution, spiking with 229 Th, 233 U, and 236 U and subsequent purification of U and Th using wet-column chemistry with Eichrom UTEVA resin. Isotope ratios were used to calculate U-series ages according to the decay equations 39 and error propagation of statistical uncertainty. The 234 U and 2 30 Th half-lives published in Cheng et al. 40 were used. Correction for detrital contamination assumes a 232 Th/ 238 U weight ratio of 3.8 ± 1.9 and 230 Th, 234 U and 238 U in secular equilibrium, and proved to be insignificant for most samples.
Stable C and O isotopes of samples were taken along the growth axis using a vertical milling machine (Sherline) at different resolution, ranging from 1 mm for the M6_1_14 to 0.2 mm for M39_1_62 and the bottom half of M39_1_61B. Samples were measured using an automated carbonate-extraction system (KIEL IV) interfaced with a MAT253 isotope ratio mass spectrometer (IRMS; ThermoFisher Scientific) at the Helmholtz Centre Potsdam GFZ. Samples of 60-90 µg were dissolved in 103 % H 3 PO 4 at 72°C and the isotopic composition was measured on the released and cryogenically purified CO 2 . All isotopic ratios were expressed in the delta notation relative to VPDB using the international reference material NBS19. Replicate analysis of the internal reference material (C1) yielded 1 sigma precision of 0.06‰ for both δ 13 C and δ 18 O.
Samples for oxygen isotope measurements of M6_1_14 were sampled at a resolution of 0.5 mm, which refers to one sample per~8 years, and were measured using a Delta plus XL IRMS at the University of Innsbruck. The mass spectrometer is linked to an on-line Gasbench II carbonate preparation system. The H 3 PO 4 acidification temperature was 72°C (for more details, see Spötl and Mattey 41 ). Isotope ratios are reported relative to the VPDB scale, and the 1 sigma precision is 0.06 and 0.08‰ for δ 13 C and δ 18 O, respectively.
Age-depth modelling. Age-depth modelling was performed using iscam 26 , an algorithm designed to align records from the same location aiming to construct a composite record. This method correlates dated proxy signals from several, contemporaneously growing stalagmites by varying their age-depth relation within allowed errors. It determines the most probable age-depth model and calculates the age uncertainty for the combined record. This increases the signal-to-noise ratio and minimizes the age uncertainties within the overlapping periods. While we do not want to go too much into detail about the algorithm, we want to discuss the most important features and parameters used. We applied a point wise linear interpolation for our age-depth modelling. We run 1,000,000 randomly selected iterations on the time series to provide the possibility to find the best option for the age-depth realization. The best option for an age-depth model is the iteration with the highest correlation coefficient between both time series. The significance of the highest correlation coefficient is usually tested by artificially constructed time series, which have similar characteristics than the measured isotope records (e.g., same dating information, proxy resolution, memory of the time series). For each stacking step, we run 2000 artificially constructed time series. To calculate the significance limits in a reasonable time we tested each of the 2000 artificially constructed time series in the same way as the measured time series, except that we used only 1000 randomly selected iterations due to computational power constraints. For each of the 2000 artificially constructed time series, the highest correlation, found within the 1000 iterations, is used to construct a probability density function, which defines significance limits for the correlation coefficient of the measured time series. All found correlation coefficients of the best-suited age-depth realizations for the individual stacking steps exceed the 99% confidence limits, except for the oldest part of M39_1_62 and M39_1_61B, which covers the time interval between~198 and~187 ka. For this period, the best correlation is not significant with respect to the applied testing method. Most likely the reason for that is the slow growth and reduced proxy resolution and number of U-Th measurements in this period.
The resulting age-depth models show a fast growth of all stalagmites betweeñ 168 and~159 ka and slower growth rates elsewhere (Fig. S3). The resulting isotope records plotted over age agree perfectly with each other (Fig. S4). Even in the time period (~198-187 ka), where the significance test of iscam failed to produce a significant correlation coefficient, the synchronization of events is given.
Conceptual climate model experiments. In this study we apply a previously established conceptual model 13 to compare the number of stadial-interstadial transitions and their length during the past two glacial cycles based on δ 18 O data from Swiss speleothems (this study) and ice core data from Greenland 14 (NGRIP). The variables and mechanisms that are accounted for in the model include the North Atlantic sea ice cover and sea surface temperatures, the strength of the AMOC as well as Greenland and Antarctic oxygen isotope ratios and are coupled by differential equations of first order. For details on the model we refer to Boers et al. 13 . Fig. 4 Comparison of climate characteristics from proxy data and stylized climate model. Sea level 34 and boreal summer insolation (a), number of warm years per 10 ka (b) and number of stadial-interstadial transitions per 10 ka (c) according to climate data as derived from the Melchsee-Frutt cave region and NGRIP (black) in comparison with results from the conceptual model (red with orange confidence interval). In the model, too cold conditions at NH summer insolation minima prevent the destabilization of sea ice and hence the triggering of stadial-interstadial oscillations. The simulations show excellent agreement with the proxy data in terms of the frequency of stadial-interstadial transitions and thus provides evidence for the proposed explanation for the quiet glacial episodes. Using the model with a constant threshold only does not reproduce the variability observed in the proxy data (Fig. S10).
In the original model, interstadials are triggered when North Atlantic Ocean surpasses a constant temperature threshold (10°C) above which the sea ice is rapidly vanishing. This part of the model was updated by applying a variable temperature threshold to account for potential influence of the NH summer insolation. The variable threshold was derived by normalising the northern Hemisphere summer insolation by its mean and standard deviation. This curve was added on the original constant threshold in a way, that for phases with originally elevated Northern Hemisphere summer insolation levels the threshold is reduced and for periods with an originally lower insolation the threshold is enlarged.
The two model versions were independently trained to reproduce the probability density function of the NGRIP δ 18 O time series for a period from 59 to 23 ky as well as to reproduce the average stadial durations in this interval. The best parameter set is estimated by analysing repeating simulations (100,000 iterations) for the training phase. The model is then run for the last two glacial-interglacial cycles. The stochastical and non-linear behaviour of the conceptual model is accounted for by repeating simulations (1000 iterations) and averaging the results. The results for the variable (Fig. 4) and constant (Fig. S10) model versions differ, with the variable threshold approach showing better agreement with the climate proxy derived characteristics of interstadial occurrence and length (compare Fig. 4 and S10).
Determination of interstadials. Interstadials during the penultimate glacial were identified from the composite record by interpolating the record on an annual scale and subsequently moving two neighbouring 300 year windows through the data. The difference of the averages of both windows provides a measure for the rapidness of changes in the δ 18 O record (Fig. S5). Minima in the difference highlight rapid increases in the O isotope composition and mark the start of an interstadial. The standard deviation of the difference of the two neighbouring running windows serves as a threshold (Fig. S5), which must be crossed to be identified as a stadial-interstadial transition. The standard deviation is 0.47 and all local minima below the negative value of this threshold are considered as significant (Fig. S5) and labelled in Figs. 2 and 3 of the main manuscript. The following maximum after each minimum, which represent a stadial-interstadial transition, highlights the return from interstadial to stadial climate conditions. For each detected interstadial, the beginning and end is defined as the first δ 18 O value that deviates from the typical stadial δ 18 O value range and the last value before the δ 18 O values return to the isotope range defined by the subsequent stadial phase.
We used a similar approach to detect the onset of interstadials in the δ 18 O records from China 32 and Peru 33 (Fig. 2). As both records show considerable orbital-scale variability, which hampers the detection of interstadial onsets, we used detrended and normalised versions of the records. For the Peru data set, we used a stacked record, which was detrended and normalized as was provided in the original publication 33 (Fig. S6). For the Chinese record 32 (Fig. S7) we normalised the time series and detrended the signal with the NH summer insolation at 60°N as was suggested by Barker et al. 16 . Then, we used the above-mentioned method to detect abrupt transitions in both normalised records, but applied 400 year running windows, as the noise in both records is larger than in our Melchsee-Frutt speleothem record. The corresponding standard deviation is again used as a threshold to identify abrupt transitions.

Data availability
All isotope data are publicly available on the NOAA paleoclimate data set webpage https://www.ncei.noaa.gov/access/paleo-search/study/38201 (including U-series data). Other data used for Figs. 1, 2, 3 are freely available at online data repositories. The data for Fig. 4