Connecting the Greenland ice-core and U/Th timescales via cosmogenic radionuclides: testing the synchroneity of Dansgaard–Oeschger events

Abstract

During the last glacial period Northern Hemisphere climate was characterized by extreme and abrupt climate changes, so-called Dansgaard–Oeschger (DO) events. Most clearly observed as temperature changes in Greenland ice-core records, their climatic imprint was geographically widespread. However, the temporal relation between DO events in Greenland and other regions is uncertain due to the chronological uncertainties of each archive, limiting our ability to test hypotheses of synchronous change. In contrast, the assumption of direct synchrony of climate changes forms the basis of many timescales. Here, we use cosmogenic radionuclides (${}^{\mathrm{10}}\mathrm{Be}$,${}^{\mathrm{36}}\mathrm{Cl}$, ${}^{\mathrm{14}}\mathrm{C}$) to link Greenland ice-core records to$\mathrm{U}/\mathrm{Th}$-dated speleothems, quantify offsets between the two timescales, and improve their absolute dating back to 45 000 years ago. This approach allows us to test the assumption that DO events occurred synchronously between Greenland ice-core and tropical speleothem records with unprecedented precision. We find that the onset of DO events occurs within synchronization uncertainties in all investigated records. Importantly, we demonstrate that local discrepancies remain in the temporal development of rapid climate change for specific events and speleothems. These may either be related to the location of proxy records relative to the shifting atmospheric fronts or to underestimated $\mathrm{U}/\mathrm{Th}$ dating uncertainties. Our study thus highlights the potential for misleading interpretations of the Earth system when applying the common practice of climate wiggle matching.

Introduction

Precise and accurate chronologies are critical for understanding past environmental and climatic changes. Global natural and anthropogenic archives can only be directly compared through the development of robust chronological frameworks, enabling studies of the spatiotemporal dynamics of past change. These findings are crucial for understanding the nature and cause of rapid climate changes in the past and hence characterizing the dynamics and feedbacks of past and projected future climate change (Thomas, 2016). However, the applicability, precision, and accuracy of the available dating methods pose strong constraints on our ability to infer leads and lags between climate records and, ultimately, mechanisms of change in the Earth system. Instead, the situation is often reversed: climate changes such as Dansgaard–Oeschger, or DO, events (Dansgaard et al., 1969, 1993) are typically assumed to occur synchronously across the Northern Hemisphere in different climate proxies from various regions and then used as chronological tie points. This so-called “climate wiggle matching” forms the chronological basis of a large part of paleoclimate records (e.g., Bard et al., 2013; Hughen et al., 2006; Henry et al., 2016; Turney et al., 2015), especially in the marine realm where other dating methods suffer from low precision and poorly constrained biases such as the marine radiocarbon reservoir age (Lougheed et al., 2013). Furthermore, it also plays a central role for one of the most widely used dating methods in paleosciences – the radiocarbon dating method. The current radiocarbon dating calibration curve (IntCal13, Reimer et al., 2013) is constructed from accurately dated tree-ring chronologies back to 13.9 ka BP (ka BP is kilo-years before present, which is 1950 CE). Beyond this time, which encompasses all DO events, about one-fourth of the data underlying IntCal13 obtain their absolute age from climate wiggle matching.

Climate wiggle matching has the obvious drawback that the leads and lags between different climate records cannot be studied once the records have been forced to align. The approach critically rests on the assumptions that (i) the climate change indeed occurred synchronously everywhere and that (ii) the (sometimes fundamentally different) proxies in question record the changes in a similar way and without intrinsic delays. These assumptions, however, can very rarely be rigorously tested, but when they are, ample evidence is revealed that questions their universal validity. Lane et al. (2013) showed that rapid climate change in the North Atlantic region may be time transgressive with regional leads and lags of the order of a century. Nakagawa et al. (2003) argued that the onset of Greenland Interstadial 1e (GI-1e; Rasmussen et al., 2014a) occurred multiple centuries after the associated climate shift in Japan, and subsequent revisions of the underlying timescales (Staff et al., 2013; Bronk Ramsey et al., 2012; Seierstad et al., 2014) did not resolve this conundrum. Buizert et al. (2015) inferred that the Southern Ocean response to DO events is delayed by about 200 years on average, while the atmosphere around Antarctica reacted instantaneously (Markle et al., 2016). Baumgartner et al. (2014) found asynchronies between ice-core proxies for local Greenland temperature ($\mathit{\delta }$${}^{\mathrm{15}}\mathrm{N}$) and the tropical and midlatitude hydrological cycle (${\mathrm{CH}}_{\mathrm{4}}$) during some DO events. They discussed the possibility that the climate changes in polar and low-latitude regions may indeed be synchronous, but that atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ concentrations rise with a delay during some DO events because of compensating changes in the source strengths of the Northern and Southern Hemisphere wetlands. Alternatively, their findings can be explained via a real delay between Greenland climate change and the activation of ${\mathrm{CH}}_{\mathrm{4}}$ source areas during certain DO events. Fleitmann et al. (2009) reported timing differences of DO events in Greenland ice cores and speleothems, albeit largely within dating uncertainties. However, they also found significant differences between speleothem records outside their chronological uncertainties. This is complemented by a recent study showing that the duration of a stadial–interstadial transition can differ by up to 300 years between different East Asian speleothems (Li et al., 2017), emphasizing the question of whether we should expect the onset, midpoint, or end point of DO events to occur simultaneously, as this choice will lead to different results when aligning the records.

In this paper, we attempt to provide improved constraints on the paradigm of climate synchroneity. We employ cosmogenic radionuclides as a climate-independent synchronization tool between the Greenland ice-core timescale (Andersen et al., 2006; Rasmussen et al., 2006; Seierstad et al., 2014; Svensson et al., 2006, 2008; Vinther et al., 2006) and the $\mathrm{U}/\mathrm{Th}$ timescale (Broecker, 1963; Edwards et al., 1987; Cheng et al., 2013a) and strongly reduce the absolute dating error of the Greenland ice cores back to 45 000 years BP. This allows us to compare the timing of DO-type variability seen in key paleoclimate records with unprecedented precision: the Greenland ice cores and $\mathrm{U}/\mathrm{Th}$-dated (sub)tropical speleothems.

Cosmogenic radionuclides as synchronization tools

Cosmogenic radionuclides (such as ${}^{\mathrm{14}}\mathrm{C}$, ${}^{\mathrm{10}}\mathrm{Be}$, and ${}^{\mathrm{36}}\mathrm{Cl}$) are produced in a nuclear cascade that is triggered when galactic cosmic rays (GCRs) collide with the Earth atmosphere's constituents (Lal and Peters, 1967). While the GCR flux outside the heliosphere can be assumed to be constant over the past million years (Vogt et al., 1990), the flux arriving at Earth is modulated by the strength of the helio-magnetic and geomagnetic fields (Masarik and Beer, 1999). This causes the production rates of cosmogenic radionuclides to be inversely related to changes in solar activity and/or the strength of the geomagnetic field. This modulation effect leaves a globally synchronous, externally forced signal in cosmogenic radionuclide records around the world. Hence, they can serve as a powerful synchronization tool for climate archives from different regions. The challenge lies in estimating potential non-production-related impacts on radionuclide concentrations in a given archive that may result from geochemical and meteorological processes.

After production, ${}^{\mathrm{14}}\mathrm{C}$ is oxidized to ${}^{\mathrm{14}}{\mathrm{CO}}_{\mathrm{2}}$ and enters the carbon cycle. Changing ${}^{\mathrm{14}}\mathrm{C}$ production rates thus alter the atmospheric ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio (expressed as per mille $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$, which is ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ corrected for fractionation and decay relative to a standard, denoted $\mathrm{\Delta }$ in Stuiver and Polach, 1977). Due to carbon cycle effects, these variations in $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ are dampened and delayed with respect to the causal production rate changes (Siegenthaler et al., 1980; Roth and Joos, 2013). In addition to variable production rates, changes in the exchange rates between the different carbon pools can alter $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. The world's oceans in particular have a significantly lower $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ than the contemporary atmosphere due to their long carbon residence time (Craig, 1957). Thus, variations in the ${}^{\mathrm{14}}\mathrm{C}$ exchange rates between the ocean and the atmosphere will alter atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ independent of production rate changes.

${}^{\mathrm{10}}\mathrm{Be}$ attaches to aerosols and is transported from the stratosphere to the troposphere within 1–2 years (Raisbeck et al., 1981), mainly via midlatitude tropopause breaks (Heikkilä et al., 2011). It has no active geochemical cycle, so its atmospheric concentration is a more direct recorder of production rate changes compared with $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. However, ${}^{\mathrm{10}}\mathrm{Be}$ transport and deposition in the troposphere are guided by local meteorology and thus susceptible to changes thereof (Heikkilä and Smith, 2013; Pedro et al., 2011). This can cause variations in ${}^{\mathrm{10}}\mathrm{Be}$ records that are not related to production rate changes. Furthermore, a so-called “polar bias” (i.e., an overrepresentation of polar as opposed to global production rate changes) has been proposed for ice-core records (Bard et al., 1997). This would lead to subdued geomagnetic and enhanced solar modulation of ice-core radionuclide records due to the geometry of the geomagnetic field. However, there is no consensus in different empirical studies and modeling experiments on whether this effect is present and the results may also vary regionally (Bard et al., 1997; Heikkilä et al., 2009a; Pedro et al., 2012; Adolphi and Muscheler, 2016; Muscheler and Heikkilä, 2011; Field et al., 2006).

The transport and deposition of ${}^{\mathrm{36}}\mathrm{Cl}$ in its aerosol phase are comparable to ${}^{\mathrm{10}}\mathrm{Be}$. However, in addition to an aerosol phase, ${}^{\mathrm{36}}\mathrm{Cl}$ also has a gaseous phase (H${}^{\mathrm{36}}\mathrm{Cl}$), which is likely dominant in the stratosphere (Zerle et al., 1997). In the troposphere, the partitioning between the aerosol and gas phase is not well understood. It may vary in space and time (Lukasczyk, 1994) and can change rapidly depending on pH (Watson et al., 1990). The gaseous H${}^{\mathrm{36}}\mathrm{Cl}$ phase can also be lost from acidic ice in low accumulation sites after deposition, which is less relevant for the high accumulation sites studied here (Delmas et al., 2004). In Greenland, similar to ${}^{\mathrm{10}}\mathrm{Be}$, the dominant deposition process of ${}^{\mathrm{36}}\mathrm{Cl}$ is wet deposition (Heikkilä et al., 2009b), which is supported by the overall similarity of ${}^{\mathrm{36}}\mathrm{Cl}$ and ${}^{\mathrm{10}}\mathrm{Be}$ variations recorded in ice cores (Wagner et al., 2001b; Muscheler et al., 2005).

As a result, all three radionuclides depend on the same production mechanism, which causes their production rates to covary globally. This signal can be exploited for global synchronization of paleorecords from natural archives. However, to identify these common changes, their different geochemistry needs to be accounted for. In the case of radiocarbon this is achieved through carbon cycle modeling to deconvolve the effects of the carbon cycle on the relation between ${}^{\mathrm{14}}\mathrm{C}$ production rates and $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ (Muscheler et al., 2004). For ${}^{\mathrm{10}}\mathrm{Be}$ and ${}^{\mathrm{36}}\mathrm{Cl}$, fluxes can be calculated from ice accumulation rates. This provides a first-order correction for changing paleoprecipitation rates on the ice sheet and their influence on the radionuclide concentrations. In reality, aerosol transport to the ice sheet is more complex and depends on changes in transport velocity, pathways, and scavenging effects en route (Schüpbach et al., 2018), which are difficult to constrain for ${}^{\mathrm{10}}\mathrm{Be}$ due to its stratospheric origin. Instead, comparisons of fluxes and concentrations to other climate proxies can inform us about potential climate influences on ${}^{\mathrm{10}}\mathrm{Be}$ and${}^{\mathrm{36}}\mathrm{Cl}$ transport and deposition (Adolphi and Muscheler, 2016). It is currently not possible to quantitatively correct either of the radionuclides for these non-production influences since neither past carbon cycle changes nor atmospheric circulation changes are sufficiently well known. However, the potential amplitude of non-production rate changes can be assessed through sensitivity experiments and added as an uncertainty for the production rate signal (Adolphi and Muscheler, 2016; Köhler et al., 2006).

The potential of this synchronization tool has been demonstrated multiple times to infer differences between the tree-ring and ice-core timescales (Adolphi and Muscheler, 2016; Muscheler et al., 2014a; Southon, 2002), test the accuracy of the radiocarbon calibration curve (Adolphi et al., 2017; Muscheler et al., 2008, 2014b), and synchronize ice cores from both hemispheres (Raisbeck et al., 2007, 2017).

Methods and data

Ice-core data

The ice-core ${}^{\mathrm{10}}\mathrm{Be}$ and ${}^{\mathrm{36}}\mathrm{Cl}$ data used in this study are shown in Fig. 1. We focus on records that have been robustly linked to the GICC05 timescale (Andersen et al., 2006; Rasmussen et al., 2006, 2008; Seierstad et al., 2014; Svensson et al., 2008). Hence, the majority of the data stems from the deep Greenland ice cores GRIP, GISP2, and NGRIP. In addition, we use Antarctic ${}^{\mathrm{10}}\mathrm{Be}$ fluxes from EDC, EDML, and Vostok that have been anchored to GICC05 by matching the solar variability present in all ${}^{\mathrm{10}}\mathrm{Be}$ records and volcanic tie points (Raisbeck et al., 2017).

Data used in this study. (a–g) Individual ice-core records of GRIP ${}^{\mathrm{10}}\mathrm{Be}$ (Baumgartner et al., 1997b; Muscheler et al., 2004; Wagner et al., 2001a; Yiou et al., 1997; Adolphi et al., 2014), GRIP ${}^{\mathrm{36}}\mathrm{Cl}$ (Baumgartner et al., 1997a, 1998; Wagner et al., 2001b; Wagner et al., 2000), GISP2 ${}^{\mathrm{10}}\mathrm{Be}$ (Finkel and Nishiizumi, 1997), and ${}^{\mathrm{10}}\mathrm{Be}$ from EDC, EDML, Vostok, and NGRIP (all Raisbeck et al., 2017). Each record represents deposition fluxes (green) and “climate-corrected” fluxes (purple; see text). In addition, each panel contains the stack of all ice-core records (black; see text). (h) ${}^{\mathrm{10}}\mathrm{Be}$ production rates modeled from two geomagnetic field intensity reconstructions: GLOPIS (green; Laj et al., 2004) and based on Black Sea sediments (purple; Nowaczyk et al., 2013) using the production rate model by Herbst et al. (2016). The ice-core radionuclide stack is shown in black. All records in (a–h) are shown on the GICC05 timescale (Seierstad et al., 2014) and normalized to (i.e., divided by) their mean. (i) Absolutely dated ${}^{\mathrm{14}}\mathrm{C}$ data from Lake Suigetsu (yellow; Bronk Ramsey et al., 2012), Hulu Cave (blue; Southon et al., 2012), Bahamas speleothems (purple; Hoffmann et al., 2010), and various tropical coral datasets (Bard et al., 1998; Cutler et al., 2004; Durand et al., 2013; Fairbanks et al., 2005; shown in light blue, olive, red, and green, respectively). The black lines encompass the $\pm \mathrm{1}\mathit{\sigma }$ uncertainties of IntCal13 (Reimer et al., 2013).

[Figure omitted. See PDF]

By calculating fluxes we make a first-order correction for the changing snow accumulation rates between stadials and interstadials and their influence on radionuclide concentrations (Wagner et al., 2001b; Johnsen et al., 1995; Rasmussen et al., 2013; Finkel and Nishiizumi, 1997). The accumulation rates for each ice core are based on their annual layer thickness – derived from their individual timescales – corrected for ice thinning. For the Greenland ice cores this thinning function is based on a 1-D ice flow model (Dansgaard and Johnsen, 1969; Johnsen et al., 2001, 1995; Seierstad et al., 2014). For the Antarctic ice cores we use the strain rate derived from the Bayesian ice-core dating effort AICC12 (Veres et al., 2013). These strain rates are inherently uncertain, and independently derived accumulation rate estimates differ by up to 10 %–20 % in the glacial (Gkinis et al., 2014; Rasmussen et al., 2013; Guillevic et al., 2013). However, these differences are largely systematic and change only on multi-millennial timescales. The shorter-term changes in accumulation rates are a more direct function of the timescale that determines the age–depth relationship, and thus annual layer thickness, and is very precise for increments of the core (Rasmussen et al., 2006). This is important to note, as we mainly exploit production rate changes on centennial-to-millennial timescales for synchronization.

To test for additional climate influences on ${}^{\mathrm{10}}\mathrm{Be}$ or ${}^{\mathrm{36}}\mathrm{Cl}$ deposition in the ice cores, we followed the approach by Adolphi and Muscheler (2016). For each ice core we calculated multiple linear regression models using $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ and snow accumulation rates as predictors for ${}^{\mathrm{10}}\mathrm{Be}$ (${}^{\mathrm{36}}\mathrm{Cl}$) fluxes and subtracted the obtained model from the ${}^{\mathrm{10}}\mathrm{Be}$ (${}^{\mathrm{36}}\mathrm{Cl}$) data. We denote the resulting record as the “climate-corrected flux” (Flux${}_{\mathrm{c}}$). This approach may correct climate effects on ${}^{\mathrm{10}}\mathrm{Be}$ (${}^{\mathrm{36}}\mathrm{Cl}$) deposition insufficiently, or it may overcorrect them, so it cannot be assumed per se that the resulting record is more reliable than the original fluxes. Nevertheless, it provides a first-order sensitivity test for the ice-core records with respect to climate-related transport and depositional effects on ${}^{\mathrm{10}}\mathrm{Be}$ (${}^{\mathrm{36}}\mathrm{Cl}$) fluxes.

To combine all ice-core records, we calculated their mean (denoted as “stack”, Fig. 1) using Monte Carlo bootstrapping (Efron, 1979). Using seven ice-core records in two versions (flux and flux${}_{\mathrm{c}}$) yields a total number of 14 samples. In each iteration, 14 samples are randomly drawn (with replacement; i.e., each record can be drawn multiple times), perturbed within measurement errors, and stacked. By repeating this procedure 1000 times we obtain an average relative standard deviation of 8 % between the derived stacks, which is comparable to the measurement uncertainty of individual measurements but larger than the expected error of the mean; this points to systematic differences between the records. For the period over which we have data from both hemispheres this standard deviation is only slightly higher (10 %). Even though this is only a relatively short period (see Fig. 1), it contains multiple DO events that are expressed differently in Northern and Southern Hemisphere climate. Thus, this agreement can serve as an indication that climate effects do not dominate the signal.

Radiocarbon data

For the purpose of this study we have to focus on radiocarbon records that are absolutely dated. Furthermore, the length and sampling resolution of the records need to be sufficient to resolve centennial-to-millennial production rate changes. The records that fulfill these criteria are shown in Fig. 1 and comprise ${}^{\mathrm{14}}\mathrm{C}$ data from various $\mathrm{U}/\mathrm{Th}$-dated coral records (Bard et al., 1998; Durand et al., 2013; Cutler et al., 2004; Fairbanks et al., 2005), as well as ${}^{\mathrm{14}}\mathrm{C}$ measured in two speleothems (Southon et al., 2012; Hoffmann et al., 2010). In addition, we use the ${}^{\mathrm{14}}\mathrm{C}$ record from Lake Suigetsu (Bronk Ramsey et al., 2012) since the $\mathrm{U}/\mathrm{Th}$-dated records do not directly reflect atmospheric ${}^{\mathrm{14}}\mathrm{C}$, only the ocean mixed layer (corals), and in the case of speleothems a mixture of atmospheric and soil ${\mathrm{CO}}_{\mathrm{2}}$ and carbonate bedrock from above the cave. The timescale of the Lake Suigetsu record is based on varve counting, corrected for long-term systematic errors by matching its ${}^{\mathrm{14}}\mathrm{C}$ record to the ${}^{\mathrm{14}}\mathrm{C}$ variations in speleothems (Bronk Ramsey et al., 2012). Hence, it is not truly independently dated. However, similar to ice-core layer counting, this varve count adds constraints, especially on centennial timescales, so $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations on these timescales should be relatively unaffected by this tuning to the speleothem ${}^{\mathrm{14}}\mathrm{C}$ data. Thus, even though the timescale may not be independent, this record can still be used to verify the existence of $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations in the atmosphere seen in the mixed layer records.

In addition, we use the available tree-ring records back to 14 000 cal BP (calibrated before present, 1950 CE) in the revised version by Hogg et al. (2016) (not shown in Fig. 1 for clarity).

Carbon cycle modeling

To be able to compare ice-core and radiocarbon records directly we have to account for the effects of the carbon cycle. Following earlier studies (Muscheler et al., 2004, 2008), we use a box-diffusion carbon cycle model (Siegenthaler et al., 1980) to model $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ from the ice-core radionuclide records. We assume that ice-core ${}^{\mathrm{10}}\mathrm{Be}$ (${}^{\mathrm{36}}\mathrm{Cl}$) variations are proportional to ${}^{\mathrm{14}}\mathrm{C}$ production rate changes (see also the following section) and model $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies from each realization of the ice-core stack, as well as the single ice-core records (Fig. 2). It can be seen that the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ records from the individual ice-core records differ in their long-term trends since the carbon cycle integrates over time so that relatively small but systematic differences in the radionuclide fluxes (possibly arising from uncertainties in the strain rates) have a significant effect on longer timescales. However, all records show the same overall evolution of $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. Furthermore, especially when subtracting the long-term trend and isolating variations on timescales shorter than 5000 years, the agreement is very high (on average within 15 ‰ at $\mathrm{1}\mathit{\sigma }$; Fig. 2b), which is the part of the signal that we will be exploiting in our synchronization effort.

Production rate ratio

Modeling $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ values from ${}^{\mathrm{10}}\mathrm{Be}$ measurements is based on the assumption that ${}^{\mathrm{10}}\mathrm{Be}$ and ${}^{\mathrm{14}}\mathrm{C}$ production rate changes are proportional to each other. However, different production rate models differ in their sensitivity of ${}^{\mathrm{14}}\mathrm{C}$ and ${}^{\mathrm{10}}\mathrm{Be}$ production rate changes to variations in the geomagnetic field (Cauquoin et al., 2014). For a given geomagnetic field change, the production rate model by Masarik and Beer (1999, 2009) yields 30 %–50 % lower ${}^{\mathrm{10}}\mathrm{Be}$ production rate changes than the calculations by Poluianov et al. (2016) and Herbst et al. (2016). For ${}^{\mathrm{14}}\mathrm{C}$, on the other hand, all models yield roughly similar amplitudes. This leads to differences in the ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ production rate ratio for a given change in the geomagnetic field. If Masarik and Beer (1999) are correct, the variations in ice-core ${}^{\mathrm{10}}\mathrm{Be}$ records have to be upscaled by 30 %–50 % to be proportional to ${}^{\mathrm{14}}\mathrm{C}$ production rate changes, while no such scaling is necessary when the other production rate models are used. In addition, the amplitudes in ${}^{\mathrm{14}}\mathrm{C}$ and ${}^{\mathrm{10}}\mathrm{Be}$ may differ due to the presence of polar bias (see Sect. 2). If this effect was present, then geomagnetic field changes should cause bigger variations in ${}^{\mathrm{14}}\mathrm{C}$ than ${}^{\mathrm{10}}\mathrm{Be}$.

Since the presence of a polar bias is debated and the physical reason for the differences between the production rate models is unresolved, we chose an empirical approach to scale the ice-core record appropriately.

Modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies from individual ice-core records (see legend; solid lines are based on radionuclide fluxes, while dashed lines are inferred from flux${}_{\mathrm{c}}$) and the realizations of the ice-core stack (black line shows the mean of all realizations, dark and light grey shading encompass 68.2 % and 95.4 % probability ranges). (a) The unfiltered model output. (b) The records after variations with frequencies < 1 $/$ 5000 a${}^{-\mathrm{1}}$ have been subtracted (FFT-based filter).

[Figure omitted. See PDF]

We use three geomagnetic field intensity reconstructions around the Laschamp geomagnetic field minimum (Laj et al., 2000, 2004; Nowaczyk et al., 2013) and calculate the resulting ${}^{\mathrm{10}}\mathrm{Be}$ production rate changes using the production rate models by Masarik and Beer (1999) and Herbst et al. (2016) (Fig. 3a–c). Subsequently, we scale the ice-core ${}^{\mathrm{10}}\mathrm{Be}$ record to minimize the root mean square error (RMSE) between ice-core and geomagnetic-field-based records (Fig. 3d). It can be seen that the RMSE reaches a minimum for a ${}^{\mathrm{10}}\mathrm{Be}$ scaling factor of $\sim \mathrm{1}$ (for Masarik and Beer, 1999) and $\sim \mathrm{1.3}$ (for Herbst et al., 2016). This represents a fortunate coincidence; irrespective of which production rate model is used, the amplitude of the ice-core ${}^{\mathrm{10}}\mathrm{Be}$ variations has to be increased by approximately 30 % to match ${}^{\mathrm{14}}\mathrm{C}$. If the production rate model by Masarik and Beer is used, then the amplitude of the ice-core ${}^{\mathrm{10}}\mathrm{Be}$ record is in agreement with geomagnetic field data, but due to the higher production sensitivity of ${}^{\mathrm{14}}\mathrm{C}$ (see above), ${}^{\mathrm{10}}\mathrm{Be}$ variations have to be increased by $\sim \mathrm{30}$ %. Similarly, if the production rate model by Herbst et al. is used, then the amplitude of the ice-core ${}^{\mathrm{10}}\mathrm{Be}$ record is 30 % smaller than implied by geomagnetic field data (possibly due to a polar bias), while the sensitivity of ${}^{\mathrm{14}}\mathrm{C}$ and ${}^{\mathrm{10}}\mathrm{Be}$ is the same. Again, the net effect is that the ${}^{\mathrm{10}}\mathrm{Be}$ variations have to be scaled up by 30 % for the comparison to ${}^{\mathrm{14}}\mathrm{C}$.

Comparison of ice-core-based and geomagnetic-field-based reconstructions of ${}^{\mathrm{10}}\mathrm{Be}$ production rates. (a–c) The ice-core stack (black) in comparison to ${}^{\mathrm{10}}\mathrm{Be}$ production rates based on geomagnetic field reconstructions and two different production rate models (Herbst et al., 2016, in pink and Masarik and Beer, 1999, in green). (a) The Black Sea geomagnetic field record (Nowaczyk et al., 2013), (b) the NAPIS geomagnetic field stack (Laj et al., 2000), and (c) the GLOPIS geomagnetic field stack (Laj et al., 2004). (d) The RMSE between the ice-core data and the geomagnetic-field-based records when variations in the ice-core record are scaled by different factors ($x$ axis). The colors correspond to the production rate models. The line styles indicate the geomagnetic field records (see legend) and the symbols denote the RMSE minima.

[Figure omitted. See PDF]

The state of the carbon cycle

As mentioned in Sect. 2, a quantification of transient carbon cycle changes and their influence on $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ is challenged by insufficient knowledge of inventories and processes. The contribution of single processes to $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes over the last glacial cycle is likely within 30 ‰ and, due to compensating effects, their combination is likely not bigger than 40 ‰ (Köhler et al., 2006). Here we use the Laschamp event to estimate the state of the ocean ventilation around 40 ka BP.

The datasets underlying IntCal13 all show an increase of about 320 ‰ in $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ into the Laschamp event (Fig. 4), albeit at different absolute levels (see Fig. 1). This is $\sim \mathrm{100}$ ‰ more than the compiled IntCal13 curve itself implies. This disagreement can be explained by differences in timing and absolute $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ between the different datasets leading to smoothing and dampening of $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations during the construction of IntCal13. Also, geomagnetic field changes yield a $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ change more in line with the individual ${}^{\mathrm{14}}\mathrm{C}$ datasets than with IntCal13, even when assuming a preindustrial carbon cycle.

To estimate the mean state of the carbon cycle during this period, we run our carbon cycle model with different (constant) values of ocean diffusivity. We find that modeled and measured $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ around the Laschamp event match best in amplitude when we run the model under conditions in which ocean ventilation is reduced to $\sim \mathrm{75}$ % of its preindustrial value (Fig. 4). This is in broad agreement with previous modeling experiments (Köhler et al., 2006; Roth and Joos, 2013) and proxy data (Henry et al., 2016).

In the following, we will use this estimate for the parameterization of our model. As mentioned above, a transient adjustment of carbon cycle parameters is uncertain and will hence not be attempted. Instead, we ascribe an associated uncertainty to the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ based on the carbon cycle sensitivity experiments by Köhler et al. (2006). Furthermore, it should be noted that by only using (filtered) $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies as synchronization targets, we (i) avoid systematic carbon cycle influences on $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ levels and (ii) minimize transient carbon-cycle-related changes in $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ (Adolphi and Muscheler, 2016).

The Laschamp event in measured and modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. (a–f) $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies in macrofossils from Lake Suigetsu (yellow; Bronk Ramsey et al., 2012), tropical corals (blue; Fairbanks et al., 2005), foraminifera from Cariaco Basin sediments (red; Hughen et al., 2006), foraminifera from Iberian Margin sediments (light blue; Bard et al., 2013), Bahamas speleothems (green; Hoffmann et al., 2010), and IntCal13 (black; Reimer et al., 2013). All data are shown as anomalies to their error-weighted mean prior to the Laschamp event, i.e., the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ increase. The dashed boxes encompass the time periods and $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ uncertainties (error of the error-weighted mean) used for the definition of the pre- and post-Laschamp event levels. (g, h) Modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ using the GLOPIS (Laj et al., 2004) geomagnetic field record and the ice-core stack as production rate inputs. The different colored lines reflect different carbon cycle scenarios (see legend; PI denotes preindustrial). The conversion of geomagnetic field intensity to ${}^{\mathrm{14}}\mathrm{C}$ production rate is based on the production rate model by Herbst et al. (2016). Note that the amplitude of the ${}^{\mathrm{10}}\mathrm{Be}$ variations has been increased by 30 % as discussed in Sect. 3.3.1.

[Figure omitted. See PDF]

Synchronization – effects of the carbon cycle and the archive

The synchronization method follows Adolphi and Muscheler (2016) and is outlined and tested in detail therein. In brief, sections of modeled (ice-core-based) $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies are compared to the measured $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. For our analysis we employ high-frequency changes in $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ since carbon cycle changes have only limited effects on atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ on shorter timescales (Adolphi and Muscheler, 2016). Similarly, as shown in Fig. 2, the agreement of the different ice-core records is better on shorter timescales. In this study, we employ two types of high-pass filtering: an FFT-based high-pass filter and simple linear detrending. The choice of filter is based on the data sampling resolution. For the highly resolved tree-ring data we use a 1000-year high-pass FFT filter, while the lower-resolved and more unevenly sampled coral, speleothem, and/or macrofossil data are filtered by linear detrending to avoid the interpolation to equidistant resolution required for FFT analysis. Similarly, the high sampling resolution of the tree-ring data allows us to compare the data in 2000-year windows, while we increase the window length to 4000 and 5000 years for the lower-resolved data prior to 14 ka BP. The exact frequencies and window lengths are also given in the Results section. Using the same statistics as for radiocarbon wiggle-match dating (Bronk Ramsey et al., 2001), we then infer a probability density function (PDF) for the timescale difference between the modeled and measured $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ records. For details on the statistics of this methodology we refer the reader to Adolphi and Muscheler (2016). Here we focus instead on additional uncertainties that arise when comparing modeled atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ to ${}^{\mathrm{14}}\mathrm{C}$ records from the ocean mixed layer (corals) or speleothems.

$\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations in the atmosphere are dampened and delayed compared to the causal production rate changes. Both factors, attenuation and delay, depend on the frequency of the production rate change (Roth and Joos, 2013; Siegenthaler et al., 1980). The dampening is largest at high frequencies and decreases with longer periods. On the other hand, the apparent peak-to-peak delay between sinusoidal production rate changes and the resulting $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ change increases with increasing wavelengths. Similar effects occur when comparing atmospheric and oceanic $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes to each other: the ocean reacts to atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes with a delayed and dampened response that is wavelength dependent. Hence, we need to take these factors into account when comparing a modeled atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ record to mixed layer marine records. However, the frequency dependence of the attenuation and delay makes it difficult to explicitly correct for this since atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes vary on different timescales simultaneously. Furthermore, the coral records vary in their sampling frequency and often it is not precisely known over how much time an individual ${}^{\mathrm{14}}\mathrm{C}$ sample integrates.

The delay between $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ in the atmosphere and ocean mixed layer. (a) Modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ from the ice-core stack around the Laschamp event. The modeled atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ is shown in black, while the ocean mixed layer is shown in grey. (b) The inferred delay from our synchronization method when comparing the atmospheric to the mixed layer signal for different low-pass filters of the mixed layer signal ($x$ axis).

[Figure omitted. See PDF]

Figure 5 shows a sensitivity test regarding these effects. We modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ from the ice-core stack around the Laschamp event and compared the atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ to the mixed layer $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ in the model. To simulate the effect of varying averaging effects of the coral samples, we low-pass filtered the mixed layer signal with increasing cutoff wavelengths. For each filter, we then inferred the apparent delay between the mixed layer (i.e., the “coral”) and the atmosphere. In doing so we infer that even though the signal is dominated by a long-lasting $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ increase, the inferred delay is small ($\sim \mathrm{30}$ years) as long as the coral samples do not integrate over long times. Only when assuming that each coral sample averages over more than 1000 years do we infer delays of about 120 years. Nevertheless, this experiment also shows that within reasonable bounds of averaging, the delay of the mixed layer to atmospheric signal is limited.

The speleothem $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ reacts differently than the ocean mixed layer. The so-called dead carbon fraction (DCF) of a speleothem consists of two main contributors: (i) respired soil organic matter that is older (in ${}^{\mathrm{14}}\mathrm{C}$ years) than the atmospheric ${}^{\mathrm{14}}\mathrm{C}$ signal and (ii) carbonate bedrock that contains no ${}^{\mathrm{14}}\mathrm{C}$. Applying the model of Genty and Massault (1999), we model speleothem $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ using different assumptions on the age of the respired soil organic matter and fraction of carbonate bedrock in drip-water ${\mathrm{CO}}_{\mathrm{2}}$. We do this for two examples: (i) a speleothem with an apparent DCF (i.e., offset from the atmosphere) of 5.8 % (resembling the Hulu Cave speleothem record by Southon et al., 2012) and (ii) a speleothem with an apparent DCF of 25.7 % (resembling the Bahamas speleothem by Hoffmann et al., 2010). By assuming different ages of the soil respired carbon ($\mathit{\tau }=\mathrm{10}$–400 years; see Fig. 6), we adjust the fraction of ${}^{\mathrm{14}}\mathrm{C}$-free ${\mathrm{CO}}_{\mathrm{2}}$ so that the apparent DCF for each speleothem is matched. The age of the soil respired carbon is defined following Genty and Massault (1999): if, for example, $\mathit{\tau }=\mathrm{100}$ years, then the activity of the soil respired CO${}_{\mathrm{2}}$ is the mean of the atmospheric activity over the past 100 years prior to sampling (also accounting for decay within these 100 years). For simplicity we assume a uniform age distribution for the soil respired carbon. Subsequently, we compare the modeled speleothem $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ to the original atmospheric input using our synchronization method and plot the inferred delay (Fig. 6b). From this experiment it can be seen that the controlling factor on the inferred delay is the age of the soil respired matter that acts as an integrator (low-pass filter) of the atmospheric ${}^{\mathrm{14}}\mathrm{C}$ signal. The fraction of ${}^{\mathrm{14}}\mathrm{C}$-free carbonate has no influence on the lag between $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes in the atmosphere and the speleothem, but only dampens the amplitude of the corresponding change. Realistic ages of soil respired carbon differ from region to region, but even though some slow-cycling fractions of soil organic matter may be up to several thousand years old (Trumbore, 2000), the major contributors to soil ${\mathrm{CO}}_{\mathrm{2}}$ are considerably younger and of the order of decades (Genty et al., 2001; Fohlmeister et al., 2011).

From these experiments we conclude that our systematic matching uncertainties to coral and speleothem records are probably below 100 years. We note that this uncertainty is asymmetric since the ocean–speleothem signal cannot lead the atmosphere, so the offset is unidirectional.

Effect of varying ages of soil respired ${\mathrm{CO}}_{\mathrm{2}}$ and fractions of ${\mathrm{CO}}_{\mathrm{2}}$ from ${}^{\mathrm{14}}\mathrm{C}$ dead carbonate on the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ in speleothems. (a) Atmospheric modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ from the ${}^{\mathrm{10}}\mathrm{Be}$ stack (black) and two modeled speleothem scenarios with a net DCF of 5.8 % (warm colors) and 25.7 % (cold colors). For each speleothem, a number of different ages for the respired soil organic matter have been assumed (see legend) and the input of ${}^{\mathrm{14}}\mathrm{C}$-free ${\mathrm{CO}}_{\mathrm{2}}$ from carbonate has been adjusted to obtain the correct apparent DCF value between 39 and 40.5 ka BP. (b) The inferred delay when we apply our synchronization method to match the atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ to the speleothem record.

[Figure omitted. See PDF]

Change-point detection window for each record. For each investigated climate event and record, the change-point detection algorithm has been applied between t1 and t2. The windows have been defined visually, ensuring a sufficient amount of data prior to and after the transition. For each record, only events that are well expressed in the climate proxy records at high resolution have been investigated. For the ice-core record, t1 and t2 typically encompass 500 years prior to and after the nominal transition ages by Rasmussen et al. (2014a). The exact values have been adjusted to exclude overlap with other transitions where necessary (Erhardt et al., 2018).

Event	GICC05	MCE	Hulu $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$		Sofular $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$		Sofular $\mathit{\delta }$${}^{\mathrm{13}}\mathrm{C}$		El Condor $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$		Diamante $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$
	(years BP)		t1	t2	t1	t2	t1	t2	t1	t2	t1	t2
Holocene	11 653	99	12 453	10 503	12 703	10 703	12 703	11 003	12 453	11 203	13 403	11203
GI-1e	14 642	186	15 442	13 942	15 442	13 942	15 442	13 942	15 442	14 192	16 392	14192
GS-3 dust peak 1	24 130	645	25 380	24 080	–	–	–	–	–	–	25 780	24 630
GI-3	27 730	832	28 580	27 680	28 780	27 880	28 780	27 780	–	–	29 030	28 080
GI-4	28 850	898	30 100	28 900	30 150	29 400	30 150	29 200	29 900	29 000	30 100	29 100
GI-5.1	30 790	1008	31 540	30 790	–	–	–	–	31 590	30 740	32 040	30 840
GI-5.2	32 450	1132	33 300	32 200	33 100	32 400	33 300	32 200	33 250	32 000	33 050	32 450
GI-6	33 690	1195	34 590	33 640	34 740	33 690	34 990	33 540	34 240	33 490	–	–
GI-7c	35 430	1321	36 680	34 980	36 380	35 480	36 380	35 230	36 230	34 880	36 480	34 980
GI-8c	38 170	1449	39 420	37 420	39 420	37 220	39 120	37 220	39 220	37 220	–	–
GI-9	40 110	1580	40 860	40 060	40 960	39 960	41 160	39 960	–	–	–	–
GI-10	41 410	1633	42 110	41 060	42 460	41 590	42 460	41 460	42 210	40 960	–	–
GI-11	43 290	1736	44 240	42 940	44 840	43 540	–	–	44 040	42 440	–	–

Synchronization of GICC05 to tree-ring and Hulu Cave records during the last deglaciation. (a) Ice-core-based modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies on the original GICC05 timescale (thin black line, light grey shading encompasses the $\pm \mathrm{10}$ ‰ uncertainty ($\pm \mathrm{1}\mathit{\sigma }$) of the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ based on the carbon cycle sensitivity experiments by Adolphi and Muscheler, 2016) and synchronized timescale (bold grey line). Tree-ring data underlying IntCal13 are shown in pink. Revised Northern Hemisphere tree-ring data according to Hogg et al. (2016) are shown in orange (preboreal pines), dark blue (late glacial pine), and light blue (Younger Dryas B chronology). New kauri $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ data by Hogg et al. (2016) are shown in purple (FIN11) and green (Towai). Hulu Cave H82 $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ data are shown as white squares. All symbols are shown with $\pm \mathrm{1}\mathit{\sigma }$ error bars. All data are FFT filtered to isolate $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations on timescales < 1000 years. (b–d) Inferred probability distributions of timescale differences between GICC05 and tree rings (orange) and Hulu Cave (black). The symbols and error bars denote means and 68.2 % and 95.4 % confidence intervals of the inferred timescale difference. Each of the lower panels refers to a 2000-year subsection of the data indicated at the top of each panel.

[Figure omitted. See PDF]

Change-point detection in climate records

To test the synchroneity of rapid climate changes, we compare the timing of DO events seen in Greenland ice cores (Andersen et al., 2004) to a number of well-known $\mathrm{U}/\mathrm{Th}$-dated speleothems that show DO-type variability from Hulu Cave (Cheng et al., 2016), Sofular Cave (Fleitmann et al., 2009), El Condor, and Cueva del Diamante (both Cheng et al., 2013b).

We use a probabilistic model to detect the onset, midpoint, and end of the rapid climate transitions in each individual record. The employed model describes the abrupt changes as a linear transition between two constant states. Any variability due to the long-term fluctuations of the climate records around the transition model is described by an AR(1) process that is estimated in conjunction with the transition model. The model is independently fitted to windows of data on their individual timescales (Table 1 and Fig. 13) around the rapid transitions. Inference was performed using Markov chain Monte Carl sampling (MCMC) to obtain PDFs of the timing of the onset, the length, and the amplitude of each transition in each record. Using these PDFs we can calculate delays of the onset, midpoint, and end of the climate transitions between different records, propagating the respective uncertainties of the parameters. For each record, only events that are well expressed and measured in high resolution have been fitted. The approach and inference procedure are described in more detail in Erhardt et al. (2018).

Synchronization results between 18 000 and 25 000 years BP. (a) The thin black line shows the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ curve based on the ice-core stack on its original timescale. The bold black line and grey shading show the synchronized ice-core record including assumed $\pm \mathrm{1}\mathit{\sigma }$ uncertainties of $\pm \mathrm{30}$ ‰. The different colored symbols indicate various ${}^{\mathrm{14}}\mathrm{C}$ datasets underlying IntCal13, which is shown as the green envelope. Lower panels: each panel shows PDFs of the inferred timescale difference between the ice-core stack and IntCal13 (green), a combination of all $\mathrm{U}/\mathrm{Th}$-dated records (speleothems–corals, pink), the H82 speleothem (blue), and Lake Suigetsu (yellow). Symbols of similar color show the inferred mean and 68.2 % and 95.4 % confidence intervals. Color-coded text indicates ${\mathit{\chi }}^{\mathrm{2}}$ probabilities for the goodness of fit between modeled and measured $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ curves after synchronization. Small (e.g., < 0.1) values would indicate significant disagreement. Note that all ${\mathit{\chi }}^{\mathrm{2}}$ probabilities are relatively high, indicating that our uncertainty estimate for the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ is very conservative. Each of the lower panels refers to a specific subsection of the data indicated at the top of each panel.

[Figure omitted. See PDF]

Timescale differences between GICC05 and the $\mathrm{U}/\mathrm{Th}$ timescale

In the following sections we will show the synchronization results for different time windows. We focus our analysis on three distinct windows: 10–14, 18–25, and 39–45 ka BP. The youngest window is defined by the presence of high-resolution tree-ring data for ${}^{\mathrm{14}}\mathrm{C}$ back to 14 ka BP. Going further back in time it becomes increasingly challenging to unequivocally identify common structures in the various $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ records that are suitable for synchronization because the resolution of the individual records decreases back in time, while their differences to each other grow steadily (see Fig. 1i). Hence, we focus on the well-known Laschamp event around 41 ka BP and the period between 18 and 25 ka BP, i.e., preceding the major carbon cycle changes associated with the deglaciation. We omit the period between 25 and 39 ka BP. As discussed in Reimer et al. (2013) and seen in Fig. 1i there is substantial disagreement between the datasets underlying IntCal13 at that time that are impossible to reconcile within their respective age and/or ${}^{\mathrm{14}}\mathrm{C}$ uncertainties. Hence, any structure in the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ records may also be unreliable and thus lead to erroneous synchronization results.

10 000–14 000 years BP

In the 10–14 ka BP interval, we synchronize the ice-core stack to high-resolution tree-ring and speleothem $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ data (Fig. 7). The high sampling resolution of the ${}^{\mathrm{14}}\mathrm{C}$ records allows us to focus on centennial-to-millennial $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes (< 1000 years) for which carbon cycle influences on $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ can be expected to be small (Adolphi and Muscheler, 2016). In concordance with earlier studies (Muscheler et al., 2014a) we find that GICC05 is $\sim \mathrm{65}$ years older than the tree-ring timescale at the onset of the Holocene, but that this offset vanishes over the course of the Younger Dryas interval.

While Muscheler et al. (2014a) argued that this changing offset may be in part due to errors in the timescale of the floating late glacial pines, we can now support this change in the timescale difference through the $\mathrm{U}/\mathrm{Th}$-dated speleothems: the synchronization of the ice-core stack to the H82 speleothem from Hulu Cave (Southon et al., 2012) leads to fully consistent results as inferred from the tree rings. This indicates that the most likely explanation is an ice-core layer counting bias, i.e., that the GICC05 timescale suggests too-old ages at the onset of the Holocene, but is accurate within a few decades during GI-1.

Interestingly, we do not observe any significant differences between the results stemming from tree rings and the speleothem records. As shown in Sect. 3.4, we could expect a delay in the speleothem $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ compared to the atmosphere if the respired soil organic carbon contribution to the soil ${\mathrm{CO}}_{\mathrm{2}}$ was very old. This would result in GICC05 appearing older in comparison to the speleothem than relative to the tree rings. The lack of this delay implies that the majority of the respired soil organic carbon at Hulu Cave must be younger than $\sim \mathrm{100}$ years (see Fig.  6). This is supported by the fact that the centennial $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations in the tree-ring and speleothem data have the same amplitude (Fig.  7). If old organic carbon significantly contributed to the soil ${\mathrm{CO}}_{\mathrm{2}}$, we would instead expect to see a stronger smoothing of short-term $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variations.

18 000–25 000 years BP

Due to the irregular and lower sampling resolution of the ${}^{\mathrm{14}}\mathrm{C}$ records beyond 15 000 cal BP, we chose to linearly detrend each dataset (instead of band-pass filtering) to remove offsets between the different ${}^{\mathrm{14}}\mathrm{C}$ datasets (see Fig. 1i) and highlight common variability. Furthermore, we have to increase the length of the comparison data windows to 4000 years to ensure sufficient structure in the ${}^{\mathrm{14}}\mathrm{C}$ sequences entering the comparison. Each window is detrended separately in the analysis to isolate short-term $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variability. We note, however, that detrending each ${}^{\mathrm{14}}\mathrm{C}$ dataset over the entire timeframe (18–25 ka BP) instead does not alter the results significantly. Compared to the high-frequency $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes studied between 10 and 14 ka BP, the longer-term variations used for synchronization here may have been increasingly affected by carbon cycle changes. To account for this, we increase the uncertainty estimate of the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes to $\pm \mathrm{30}$ ‰ ($\pm \mathrm{1}\mathit{\sigma }$), which is sufficiently large to account for estimated carbon-cycle-driven $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes from modeling experiments during the entire glacial (Köhler et al., 2006). We note that this is a conservative estimate, given that during this period neither modeling (Köhler et al., 2006; Muscheler et al., 2004) nor data (Eggleston et al., 2016) suggest large carbon cycle changes.

It can be seen in Fig. 8 that it is challenging to infer robust covariability in multiple ${}^{\mathrm{14}}\mathrm{C}$ records. However, the millennial evolution of $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ does show common changes in the 18–25 ka BP interval. Synchronizing the ice-core stack to data from (i) Hulu Cave H82 speleothem, (ii) Lake Suigetsu macrofossils, (iii) the IntCal13 stack, or (iv) a combination of all $\mathrm{U}/\mathrm{Th}$-dated records (speleothems–corals) leads to consistent results within uncertainties for each choice of time windows: all records imply that GICC05 shows younger ages compared to the ${}^{\mathrm{14}}\mathrm{C}$ records around this time.

Close-up of measured and modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies between 21 and 23 ka BP. The thin grey line shows modeled atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ from the ice-core stack on the GICC05 timescale. The bold black and dashed red lines show the modeled atmospheric and ocean mixed layer $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ curves after synchronization to the ${}^{\mathrm{14}}\mathrm{C}$ records (yellow: Lake Suigetsu; blue: Hulu Cave; purple and white: corals. The inset panel shows the PDF of the inferred timescale difference between GICC05 and the combination of all ${}^{\mathrm{14}}\mathrm{C}$ records. The black line is based on using only the modeled atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. The red dashed line is based on comparing coral and speleothem data to the modeled mixed layer $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ and Lake Suigetsu data to modeled atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. The green line shows modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ based on geomagnetic field changes.

[Figure omitted. See PDF]

The most significant structure that is present in all measured and modeled ${}^{\mathrm{14}}\mathrm{C}$ records during this time is the centennial $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ increase around 22.1 ka BP (see Fig. 9). Comparing the ice-core stack to $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ between 21 and 23 ka BP indicates an offset of $\sim \mathrm{550}$ years between GICC05 and the $\mathrm{U}/\mathrm{Th}$ timescale around this time (GICC05 being younger). To account for the potential delay of coral and speleothem $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ compared to the atmosphere, we also modeled the mixed layer $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ signal from the ice-core stack and synchronized this signal to the measured ${}^{\mathrm{14}}\mathrm{C}$ data (Fig. 9). As discussed in Sect. 3.4, we find very little difference in the inferred timing since the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ variation is relatively rapid (centuries). Comparing the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies to geomagnetic field data shows that a small part of the longer-term development of this structure is probably driven by geomagnetic field changes. The amplitude ($\sim \mathrm{50}$ ‰) and short duration (centuries) of the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ increase, however, suggest that this change is mainly driven by a series of strong solar minima comparable to the grand solar minimum period around the onset of the Younger Dryas (Muscheler et al., 2008). We used this tie point (Fig. 9) in the final synchronization as it is the best-defined feature in this time interval and consistent within error with the estimates shown in Fig. 8.

39 000–45 000 years BP

Our oldest tie point is the previously discussed Laschamp event around 41 ka BP. The only independently and absolutely dated ${}^{\mathrm{14}}\mathrm{C}$ record around this time that has a sufficient sampling resolution is the Bahamas speleothem by Hoffmann et al. (2010). While offset in absolute $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ (see Fig. 1), the $\mathrm{U}/\mathrm{Th}$-dated coral data support the amplitude and timing of the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ increase seen in the speleothem even though precise synchronization is hampered by the low sampling resolution of the corals. The Lake Suigetsu record is characterized by large uncertainties and scatter around this time. As discussed in Sect. 3.3.2, IntCal13 is smoothed around Laschamp, having a smaller amplitude and a less sharp rise in $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$. For this tie point, we merely remove the error-weighted mean between 39 and 45 ka BP from each dataset, since detrending would remove the largest part of the signal. Hence, there are large $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ modeling uncertainties associated with unknown carbon cycle changes, and we assume a Gaussian $\pm \mathrm{1}\mathit{\sigma }$ error of 50 ‰, which we consider conservative since sensitivity experiments imply that the impact of carbon cycle changes on $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ was likely below 40 ‰ (Köhler et al., 2006).

Synchronization of ${}^{\mathrm{10}}\mathrm{Be}$ and ${}^{\mathrm{14}}\mathrm{C}$ around the Laschamp event. The black lines encompass the modeled $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies ($\pm \mathrm{1}\mathit{\sigma }$) from the ice-core data shifted by $+\mathrm{252}$ years (68.2 % confidence interval $=$ $-\mathrm{103}$ to 477 years) according to their best fit to the speleothem ${}^{\mathrm{14}}\mathrm{C}$ data. The green patch shows the $\pm \mathrm{1}\mathit{\sigma }$ envelope of IntCal13. The blue and purple symbols show $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ from Bahamas speleothem and corals, respectively. The yellow symbols show $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ anomalies based on Lake Suigetsu macrofossils. All datasets have been centered to 0 ‰ by subtracting the error-weighted mean of each dataset. The inset shows the PDF of the inferred age differences between the ice-core data and IntCal13 (green), Lake Suigetsu (yellow), and the Bahamas speleothem (blue). The dashed blue line corresponds to age differences from the modeled mixed layer $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ and the Bahamas speleothem.

[Figure omitted. See PDF]

Synchronizing the ice-core stack to the speleothem, Lake Suigetsu, and IntCal13 data yields significantly different results. We infer that GICC05 produces ages about 250 years younger than the $\mathrm{U}/\mathrm{Th}$-dated speleothem data (Fig. 10). The IntCal13 record, however, implies a larger difference of $\sim \mathrm{1000}$ years. Using Lake Suigetsu data, on the other hand, leads to multiple probability peaks, two of which are in agreement with the speleothem and one with the IntCal13 record. The large scatter of the Lake Suigetsu data, however, leads to poor statistics (low ${\mathit{\chi }}^{\mathrm{2}}$ probabilities). Furthermore, the Lake Suigetsu timescale is only constrained by varve counting back to 39 ka BP and based on extrapolation for older sections (Bronk Ramsey et al., 2012); it hence provides less precise constraints on the timing of the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ increase.

Comparison of the ice-core stack (blue) to Ar–Ar dates of the Laschamp excursion (yellow: Singer et al., 2009; pink: Laj et al., 2014) and relative geomagnetic field intensity (black; NRM/ARM; reversed $y$ axis) from a $\mathrm{U}/\mathrm{Th}$-dated speleothem (Lascu et al., 2016). The individual speleothem $\mathrm{U}/\mathrm{Th}$ dates are shown on the bottom of the figure with their $\pm \mathrm{2}\mathit{\sigma }$ uncertainties. Each panel shows a different shift of GICC05 according to the results from Fig. 10.

[Figure omitted. See PDF]

To test which of these links is the most likely we turn to independent radiometric ages of the Laschamp excursion. Pooled Ar–Ar, K–Ar, and $\mathrm{U}/\mathrm{Th}$ ages on lava flows place the period of (nearly) reversed field direction at $\mathrm{40}\mathrm{700}\pm \mathrm{950}$ years BP (Singer et al., 2009) or $\mathrm{41}\mathrm{300}\pm \mathrm{600}$ years BP (Laj et al., 2014). In addition, a North American speleothem provides a $\mathrm{U}/\mathrm{Th}$-dated transient evolution of the geomagnetic field (Lascu et al., 2016), with the lowest intensities occurring at $\mathrm{41}\mathrm{100}\pm \mathrm{350}$ years BP. Comparing the ice-core ${}^{\mathrm{10}}\mathrm{Be}$ stack to these data clearly shows that all of these records rule out the $+\mathrm{1000}$ year time shift implied by IntCal13, as it would induce a significant disagreement between radiometrically dated magnetic field records and the dating of the ${}^{\mathrm{10}}\mathrm{Be}$ peak in the ice cores (Fig. 11). We hence argue that the 252-year offset inferred from the comparison to the Bahamas speleothem is the most likely estimate of the timescale difference between GICC05 and the $\mathrm{U}/\mathrm{Th}$ timescale around this time. Similar to before, assuming that the speleothem represents a mixed layer signal instead of direct atmospheric $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ does not significantly affect the inferred timescale differences (see Fig. 10 inset, blue dashed line).

Transfer function between the $\mathrm{U}/\mathrm{Th}$ timescale and GICC05. The transfer function is shown in black with dark and light grey shading encompassing its 68.2 % and 95.4 % confidence intervals. The black dots with error bars show the match points used between ${}^{\mathrm{14}}\mathrm{C}$ and ${}^{\mathrm{10}}\mathrm{Be}$. The red star shows the difference between ages of a glacial kauri tree ${}^{\mathrm{14}}\mathrm{C}$ sequence on Lake Suigetsu ${}^{\mathrm{14}}\mathrm{C}$ and GRIP ${}^{\mathrm{10}}\mathrm{Be}$ (Turney et al., 2016). The blue open square shows the age difference between the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age of the Campanian Ignimbrite (Giaccio et al., 2017) and a tentatively associated spike in the GISP2 ${\mathrm{SO}}_{\mathrm{4}}$ record (Fedele et al., 2007) on the GICC05 timescale (Seierstad et al., 2014).

[Figure omitted. See PDF]

Transfer function

To construct a continuous transfer function between GICC05 and the $\mathrm{U}/\mathrm{Th}$ timescale we apply a Monte Carlo approach. Each iteration consists of (i) randomly sampling the PDFs at each tie point and (ii) interpolating between the tie points using an AR process that is constrained by the GICC05 maximum counting error (MCE). We use the tie points shown in Figs. 7, 9, and 10, i.e., three tie points between ice cores and tree rings during the deglaciation, one tie point between ice cores and the combination of corals, speleothems, and Lake Suigetsu during the LGM, and one tie point between ice cores and the Bahamas speleothem around the Laschamp event. For the interpolation, we use the time derivative of the MCE (i.e., its growth rate) as an incremental error estimate. During periods when the growth rate is > 0, GICC05 may be stretched (compressed), while a growth rate of 0 does not allow this, independent of what the absolute MCE is at that time. By multiplying this growth rate with a random realization of an AR process ($\mathrm{\Phi }=\mathrm{0.9}$, $\mathit{\sigma }=\mathrm{1}$), we simulate how much of that uncertainty has been realized in each iteration of the Monte Carlo simulation. Subsequently integrating over the resulting time series of simulated miscounts, we again obtain an absolute error estimate, i.e., one possible realization of the MCE. The parameters for the AR process were chosen so that the simulated realization of the MCE explores the whole absolute counting error space, without frequently exceeding the permitted growth rate of the MCE. A larger $\mathrm{\Phi }$ would increase interpolation uncertainty, but also frequently violate the constraints of the layer count. A smaller $\mathrm{\Phi }$, on the other hand, would decrease the uncertainty due to shorter decorrelation length (see also the discussion in Rasmussen et al., 2006). In each iteration, this realization is then anchored at the sampled tie points (step i) by linearly correcting the offset between the sampled tie points and the simulated counting error. Hence, this procedure provides us with a correlated interpolation uncertainty over time, taking into account some of the constraints provided by the ice-core timescale itself, but giving priority to our inferred tie points. We note that this treatment of the MCE as an AR process leads to larger interpolation errors compared to assuming a white noise model, which would lead to very small uncertainties that average out over a long time (see also discussion in Rasmussen et al., 2006). Furthermore, we treat the MCE as $\pm \mathrm{1}\mathit{\sigma }$ instead of $\pm \mathrm{2}\mathit{\sigma }$ as proposed by Andersen et al. (2006), which additionally increases our interpolation error. We stress that this procedure does not aim to provide a realistic model of the ice-core layer counting process and its uncertainty, which is clearly more complex (see Andersen et al., 2006; Rasmussen et al., 2006), nor should it be interpreted such that the MCE was a $\mathrm{1}\mathit{\sigma }$ uncertainty. However, our approach allows us to infer a conservative estimate of the interpolation uncertainty, while at the same time it takes advantage of the fact that GICC05 is a layer-counted timescale and hence cannot be stretched or compressed outside realistic bounds. This procedure was repeated 300 000 times, which was found sufficient to obtain a stationary solution leading to 300 000 possible timescale transfer functions.

Figure 12 shows the resulting mean transfer function along with its confidence intervals. First, it can be seen that all tie points fall into the uncertainty envelope of GICC05. The implied change in the timescale difference between the youngest two tie points (i.e., over the course of GS-1) and between 13 000 and 22 000 years BP is slightly larger than allowed by the MCE, although the latter is consistent within the uncertainties of the tie point at 22 000 years BP. We can see that the use of the MCE to determine the interpolation error leads to small uncertainties wherever the change in the timescale difference is large (e.g., over the 13 000–22 000 years BP interval): stretching GICC05 by as much as the counting error allows, requires that every uncertain layer has in fact been a real annual layer, leaving little room for additional error. Between 22 000 and 42 000 years BP, the interpolation uncertainties are determined by the MCE and thus grow or shrink at a rate determined by the MCE.

Our results are in very good agreement with the results by Turney et al. (2016) around Heinrich 3. In this study, a kauri-tree ${}^{\mathrm{14}}\mathrm{C}$ sequence was calibrated onto Lake Suigetsu ${}^{\mathrm{14}}\mathrm{C}$ and also matched on GICC05 via ${}^{\mathrm{10}}\mathrm{Be}$. The difference of the inferred ages (i.e., kauri on Suigetsu vs. kauri on GICC05) matches with our proposed transfer function (red star in Fig. 12).

Figure 12 also shows the inferred offset between the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age of the Campanian Ignimbrite (Giaccio et al., 2017) and a tentatively attributed ${\mathrm{SO}}_{\mathrm{4}}$ spike in the GISP2 ice core (Fedele et al., 2007). Even though it obviously requires a well-characterized tephra find in the ice cores to ensure that the ${\mathrm{SO}}_{\mathrm{4}}$ peak is indeed associated with the Campanian Ignimbrite, at least from a chronological point of view, our transfer function does not preclude this link. However, no matching shards were identified in this period (Bourne et al., 2013).

Timing of abrupt climate changes in different climate records. The climate archive and proxy are indicated in each panel. The black lines show the mean of the fitted ramps and their 95 % confidence intervals (dashed). The dots mark the midpoint of the mean transition. The $\mathrm{U}/\mathrm{Th}$ dates and their $\pm \mathrm{1}\mathit{\sigma }$ uncertainties of each climate record are shown at the bottom of the figure in color coding corresponding to the respective climate record. Each time series is shown on its original timescale not applying any synchronization.

[Figure omitted. See PDF]

The timing of DO events

To investigate the synchroneity of climate changes recorded in different parts of the globe, we compare ice-core data to a selection of well-dated speleothem records. The well-known Hulu–Dongge Cave records have become iconic blueprints for intensity changes in the East Asian summer monsoon (EASM) anchored on a precise $\mathrm{U}/\mathrm{Th}$ timescale (Cheng et al., 2016; Dykoski et al., 2005; Wang et al., 2001). The speleothem records from Cueva del Diamante and El Condor reflect changes in precipitation amount over eastern Amazonia associated with the South American monsoon (Cheng et al., 2013b). The speleothem records from Sofular Cave, Turkey, are not straightforward in their mechanistic interpretation, but likely reflect a mix of temperature and seasonality of precipitation ($\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$), type and density of vegetation, soil microbial activity ($\mathit{\delta }$${}^{\mathrm{13}}\mathrm{C}$), and hence effective moisture and temperature (Fleitmann et al., 2009). While this list of speleothem data can certainly be expanded in future studies, we chose these four speleothem records from three different regions that are all well dated and sensitive to the position of the ITCZ and compare it to the ice-core records. We used the NGRIP Ca record (Bigler, 2004), which shows the largest signal-to-noise ratio across DO events (compared to, e.g., $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$), making their identification more precise. In addition, the Ca aerosols originate from Asian dust sources (Svensson et al., 2000) and are thus more directly related to Asian hydroclimate (Schüpbach et al., 2018), making them potentially more comparable to, for example, the Hulu Cave record. Potential phasing differences between different climate proxies in the ice core are small compared to our synchronization uncertainties (Steffensen et al., 2008).

Timing differences of the onset (a, b), midpoint (c, d), and end (e, f) of rapid climate changes in NGRIP and speleothems (colored PDFs; see legend) and the timescale transfer function inferred from radionuclide matching (black line and grey shadings as in Fig. 12). (a, c, e) The PDFs of timing differences including only uncertainties from the determination of the change points in the climate records; panels (b, d, f) also include the speleothem dating uncertainties.

[Figure omitted. See PDF]

Figure 13 shows the ice-core and speleothem climate records on their original individual timescales, along with the fitted ramps to the rapid climate changes. Note that we could not fit each climate event for every record, since the method requires a minimum number of data points defining the levels before and after each transition to produce reliable estimates. Already visually, a lag of climate changes in Greenland compared to the speleothem records can be consistently identified between 20 and 35 ka BP when all records are on their original timescales. Combining the PDFs of the detected change points in Greenland and the speleothems allows us to infer a probability estimate of the timing difference between climate events in Greenland and speleothems. These differences are shown in Fig. 14 along with our transfer function based on the matching of radionuclides from Fig. 12. This comparison shows that the differences in the timing of starting point, midpoint, and end point of DO events in speleothems and ice cores largely fall within the uncertainties of our radionuclide-based timescale transfer function. Thus, rapid climate changes occur synchronously in Greenland and the (sub)tropics. Notable exceptions are (i) the transition from GS-1 to the Holocene around 11.6 ka BP, (ii) Heinrich event 2 at 24 ka BP, and (iii) DO 11 around 43 ka BP. However, there is large scatter among the different speleothem-based estimates at these events, indicating that these events are asynchronous in the different speleothem records on their respective timescales. Consequently, some of these records also imply asynchronous climate shifts with Greenland ice cores. This may either be interpreted as an indication of time-transgressive climate changes or as a bias in individual speleothems either in how climate is recorded in the speleothem or their dating (for example, through detrital thorium).

Average lead–lag between the onset of DO events in the speleothems and NGRIP. Each panel (color) shows the PDF for the lead–lag of the onset in the speleothem compared to NGRIP, averaged over all investigated DO events (i.e., excluding the GS-3 dust peak, H2). Panel (f) shows the PDF of the average of all DO events and speleothems. The dark and light shading of the PDF in each panel indicates 68.2 % and 95.4 % intervals.

[Figure omitted. See PDF]

Averaging over all DO events, we can estimate an overall probability of leads and lags. By using the individual realizations of the radionuclide-based transfer function (see Sect. 4.4) we take into account that the uncertainties of the transfer function are strongly autocorrelated. For each realization, we randomly sample the PDFs for the onset of the DO events for the ice-core and speleothem records (see Sect. 3.5), perturb the speleothem-based estimates within their $\mathrm{U}/\mathrm{Th}$ dating errors, determine the lead or lag between the DO onset in ice-core and speleothem records, and correct it for the expected lag from the realization of our transfer function. By averaging over all DO events we thus obtain a mean lag for each realization and speleothem. In addition, we combine the different speleothem-based estimates of each realization by averaging over their mean lags to obtain an overall (speleothem and DO event) mean lag. Converting the obtained lags from each realization into histograms we estimate the PDFs of average lags between ice-core and speleothem records.

Our lag estimates critically depend on our ability to fill the gaps between the widely spaced tie points, and thus on our assumptions about the ice-core layer counting uncertainty, and how well our AR(1) process model can capture these (Sect. 4). However, we note that by treating the MCE as a highly correlated $\mathrm{1}\mathit{\sigma }$ (instead of 2$\mathit{\sigma }$) uncertainty, our error estimate can be regarded as very conservative since it allows for large systematic drifts in each realization of the transfer function that will result in large errors of the mean.

The resulting PDFs of the lag between speleothems and ice cores are shown in Fig. 15. The uncertainties are mainly determined by our synchronization uncertainty. Thus, the uncertainty is only marginally reduced when averaging over all speleothems (Fig. 15, bottom). Because each realization of the transfer function varies smoothly, the offset between speleothem and ice-core records will be systematic for all speleothems in each realization and is thus only marginally reduced by averaging.

We find that all speleothem records except Cueva del Diamante (Cheng et al., 2013b) indicate synchroneity with NGRIP within $\mathrm{1}\mathit{\sigma }$ and that the delay obtained for Cueva del Diamante falls within $\mathrm{2}\mathit{\sigma }$. We note that the speleothem data from El Condor (Cheng et al., 2013b) from the same region as Cueva del Diamante do not indicate a significant lag to Greenland. Overall, our analysis cannot reject the null hypothesis of synchronous DO events in Greenland ice cores and (sub)tropical speleothems (lag: $\mathit{\mu }\pm \mathrm{1}\mathit{\sigma }=\mathrm{26}\pm \mathrm{189}$ years).

Climate changes around H2. (a–h) NGRIP Ca (Bigler, 2004) on the synchronized timescale (Fig. 14), Hulu Cave $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Cheng et al., 2016), El Condor $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Cheng et al., 2013b), Cueva del Diamante $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Cheng et al., 2013b), (e) Pacupahuain $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Kanner et al., 2012), Paixao $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Stríkis et al., 2018), Lapa Sem Fim $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Stríkis et al., 2018), and Jaragua Cave $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Novello et al., 2017). The arrows on the right-hand side of each axis point in the direction of the signature of increased precipitation on $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ through the amount effect (Dansgaard, 1964). The light grey box marks H2. The dark grey box highlights the preceding $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ anomaly in El Condor and Diamante caves. (i, j) Precipitation (color) and wind (arrows) response to freshwater forcing in the CCSM3 climate model (freshwater-only experiment of TraCE21k, all other forcings are held at 19 ka BP conditions; He, 2011). The red (blue) line depicts the latitude of the precipitation maximum during strong (weak) AMOC states. Only wind anomalies > 1 m s${}^{-\mathrm{1}}$ are plotted. The cave sites are indicated as dots. Panel (i) shows the winter (December–February) response, while panel (j) shows the summer (June–August) response. Anomalies are plotted as weak–strong AMOC mode.

[Figure omitted. See PDF]

Discussion

Our proposed transfer function quantifies the long-term differences between the Greenland ice-core and $\mathrm{U}/\mathrm{Th}$ timescale and allows for their synchronization. Even though based on only a few tie points, this can be used to evaluate the absolute dating accuracy of Greenland ice-core records during the past 45 ka BP, while maintaining the strength of their precise relative dating. In combination with similar work done for the Holocene (Adolphi and Muscheler, 2016; Muscheler et al., 2014a), the picture emerges that the GICC05 counting error may be systematic: when accumulation and data resolution are high (e.g., in parts of the Holocene), too many annual layers have been counted, whereas during periods of low accumulation (e.g., GS-1 and GS-2) there is a tendency to identify too few annual layers. In principle, this is well captured by the GICC05 uncertainty estimate as the derivative of our transfer function is (within error) consistent with the increase in the counting error. However, our results caution against the use of the GICC05 counting error as a $\mathrm{2}\mathit{\sigma }$ uncertainty as is often done in the literature. Originally, Andersen et al. (2006) pointed out that the MCE is not a true $\mathit{\sigma }$ uncertainty but proposed that a Gaussian distribution with $\mathrm{2}\mathit{\sigma }=$ MCE could serve as a pragmatic approximation. In combination with results from the Holocene (Adolphi and Muscheler, 2016) our study implies that the counting error can be strongly correlated over extended periods of time. This is in line with the discussion in Rasmussen et al. (2006), who point out that the main contribution to a potential bias in the layer count is the definition of how an annual layer is manifested in the proxy data. The data resolution and the manifestation of annual layers change between different climate states (Rasmussen et al., 2006), likely due to changes in aerosol transport and deposition resulting from variations in the atmospheric circulation and seasonality of precipitation (Merz et al., 2013; Werner et al., 2001). According to our analysis, the largest relative (i.e., year year${}^{-\mathrm{1}}$) change in the difference between GICC05 and the $\mathrm{U}/\mathrm{Th}$ and tree-ring timescale occurs over GS-1 (11 653–12 846 years BP) and GS-2 (14 652–23 290 years BP). Both of these periods have likely been characterized by an increased relative contribution of summer precipitation to the annual ice layer (Werner et al., 2000; Denton et al., 2005), and the annual layers in the ice core have been identified in a similar way in both intervals (Rasmussen et al., 2006). In the 11–13 ka BP interval, the offset between GICC05 and the tree-ring timescale changes from $-\mathrm{60}$ (95.4 % range: $-\mathrm{77}$ to $-\mathrm{42}$) years to zero (95.4 % range: $-\mathrm{12}$ to $+\mathrm{21}$) years. During the same interval, the GICC05 maximum counting error grows by 46 years. Although consistent within the absolute error margins, this stretch of GICC05 over GS-1 thus slightly exceeds the range allowed by the GICC05 counting error. Muscheler et al. (2014a) discussed the fact that this stretch may be partly explained by errors in the placement of the oldest part of the tree-ring chronologies. However, here, we use a revised late glacial tree-ring dataset in which the different chronologies are connected much more robustly (Hogg et al., 2016). Furthermore, our analysis of the fully independent Hulu Cave ${}^{\mathrm{14}}\mathrm{C}$ data yields similar results (Fig. 7). Hence, our analyses clearly show that the GS-1 interval is about 60 years too short in the GICC05 timescale.

Between 15 and 22 ka BP, our analysis yields a change in the GICC05 offset from $+\mathrm{118}$ (95.4 % range: 2–220) years to $+\mathrm{549}$ (95.4 % range: 207–670) years, while the GICC05 counting error grows by 335 years. Thus, again, our transfer function changes a little faster than the maximum counting error allows during this interval. We note that our ${}^{\mathrm{14}}\mathrm{C}$–${}^{\mathrm{10}}\mathrm{Be}$ match point around 22 000 years BP has a relatively low signal-to-noise ratio in the ${}^{\mathrm{14}}\mathrm{C}$ data (see Figs. 8–9) and should thus be regarded as tentative. However, as shown in Fig. 8 our results are generally robust against different choices of subsets of the ${}^{\mathrm{14}}\mathrm{C}$ data and time windows. Nevertheless, it can also be seen that the estimates of the most likely age difference (i.e., the peak of the PDFs) differ slightly for different choices of ${}^{\mathrm{14}}\mathrm{C}$ data. Hulu Cave yields a most likely offset of $\sim \mathrm{325}$ years, while Suigetsu implies a bigger age difference of $\sim \mathrm{550}$ years that coincides with a secondary probability peak in the Hulu Cave PDF. We note that assuming increased amounts of old soil organic carbon contributing to speleothem formation would lead to an even stronger difference between these estimates (see Sect. 3.4). Hence, we propose an age difference of $+\mathrm{549}$ (95.4 % range: 207–670) years based on the combination of all data (Fig. 9) that is consistent within error with the estimates based on the single datasets shown in Fig. 8, but stress that this tie point should be reevaluated as new suitable ${}^{\mathrm{14}}\mathrm{C}$ data become available in the future.

Assuming that the $\mathrm{U}/\mathrm{Th}$ dates are absolute, our transfer function can be used to account for the bias in the GICC05 timescale and thus facilitate comparisons of ice-core records to other absolutely dated archives. However, we note that our synchronization does not necessarily lead to consistent timescales with radiocarbon-dated records. As discussed in Sect. 3.3.2 (Fig. 4) and Sect. 4.3 (Figs. 10 and 11), discrepancies in the datasets underlying IntCal13 can lead to erroneous structures in the calibration curve. The reduced amplitude of the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ change around the Laschamp geomagnetic field minimum in IntCal13 compared to its underlying data implies that IntCal $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ must be offset prior to and/or after the Laschamp event. This underlines the challenges in radiocarbon calibration around this time pointed out by Muscheler et al. (2014b). More recently, Giaccio et al. (2017) also pointed out that paired ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ and ${}^{\mathrm{14}}\mathrm{C}$ dating of the Campanian Ignimbrite around 40 ka BP yields inconsistent ages when the ${}^{\mathrm{14}}\mathrm{C}$ age is calibrated with IntCal13. Since IntCal13 in principle should be tied to the $\mathrm{U}/\mathrm{Th}$ age scale for sections older than 13.9 ka BP, this implies either an inconsistency between $\mathrm{Ar}/\mathrm{Ar}$ and $\mathrm{U}/\mathrm{Th}$ dating or in the reconstructed ${}^{\mathrm{14}}\mathrm{C}$ levels of the calibration curve. The latter would be congruent with the conclusions by Muscheler et al. (2014b). If the problem was indeed the IntCal ${}^{\mathrm{14}}\mathrm{C}$ reconstruction, a synchronization of ice-core ${}^{\mathrm{10}}\mathrm{Be}$ to IntCal ${}^{\mathrm{14}}\mathrm{C}$ would not resolve this bias, since the problem would not be one of chronology, but of ${}^{\mathrm{14}}\mathrm{C}$ measurement and/or archive.

Our analysis provides the first rigorous test of whether DO events recorded in speleothems and ice cores occur synchronously. We reject the hypothesis of leads and lags larger than 189 years at the one sigma level, consistent with the findings of Baumgartner et al. (2014). Since we compare to speleothem records from different regions, this also highlights that the ITCZ likely migrated synchronously (within uncertainties) over the different ocean basins and continents during the onset of DO events (Schneider et al., 2014). However, there are also differences between the different speleothem records, which could be due to limitations in their dating or related to how directly individual archives record the rapid climate changes. The most notable examples are the onset of the Holocene and GI-11, which appear to occur asynchronously in the different speleothems (see Figs. 13 and 14). Another example is the younger GS-3 dust peak in the Greenland ice cores that appears to coincide with the East Asian summer monsoon decline seen in Hulu Cave, but postdates the precipitation increase seen in El Condor and Diamante. This change in the speleothems is typically attributed to the southward shift of the ITCZ as a response to Heinrich event 2.

Figure 16 shows the period around H2. First, we note that the younger of the two GS-3 dust peaks in the Greenland ice cores (Rasmussen et al., 2014a) occurs coevally (within chronological uncertainty) with the ITCZ movement recorded by the speleothems. At this time, the East Asian summer monsoon is strongly reduced as implied by decreased Hulu Cave $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ (Cheng et al., 2016). Coevally, precipitation increases in the South American summer monsoon region (Novello et al., 2017; Stríkis et al., 2018). The records thus exhibit more pronounced stadial conditions than normally seen during (non-Heinrich) DO events. However, taken at face value, the precipitation increase at El Condor and Cueva del Diamante, the two northernmost sites shown in Fig. 16 (Cheng et al., 2013b), significantly predates the event seen in Greenland and Hulu Cave. It also predates the more southern South American sites Lapa Sem Fim (Stríkis et al., 2018) and Jaragua (Novello et al., 2017) by more than 500 years. This could either point to errors in the dating of the El Condor and Diamante speleothems or be related to their latitudinal position. A freshwater-only experiment (all other boundary conditions held constant at 19 ka BP levels) with the Community Climate System Model 3 (TraCE-MWF; He, 2011) shows that during a weak AMOC state, reduced advection of moisture from the tropical Atlantic leads to lower precipitation north of the ITCZ, while the ITCZ position over South America itself changes very little (Fig. 16). El Condor and Cueva del Diamante are both located very close to the LGM position of the ITCZ. It is hence possible that when Northern Hemisphere summer insolation reached its lowest values over the past 50 ka BP around H2, the ITCZ migrated to a position south of El Condor and Cueva del Diamante and during its transition caused the reconstructed precipitation change. As a result, the precipitation response to freshwater forcing would change sign at these cave sites. The sites located slightly further south only show a weak (Pacupahuain) or no (Paixao) response during this period, but are both characterized by increased variability. The two southernmost sites, on the other hand (Jaragua and Lapa Sem Fim), remain south of the ITCZ throughout and hence show a clear increase in precipitation coeval with the signal in Greenland and Hulu Cave. In this context, the precipitation increase in El Condor and Cueva del Diamante around 25 ka BP (i.e., prior to H2) may signify when the ITCZ transitions over the sites. The subsequent reduction in AMOC strength during H2 then leads to a decrease in precipitation in northwest South America, an increase further south, and little change in between. Tentative support for this can be drawn from the response of the El Condor and Cueva del Diamante speleothems to GI-2.2 and GI-2.1 for which, albeit weakly, the $\mathit{\delta }$${}^{\mathrm{18}}\mathrm{O}$ records imply an increase in precipitation during GI-2 opposite to their response to DO events during MIS-3 (Figs. 13, 16). Thus, this analysis indicates that the seemingly asynchronous response to climate change in different proxy records may indeed only reflect site-specific changes in the proxy response. Alternatively, we cannot rule out undetected issues with the $\mathrm{U}/\mathrm{Th}$ ages of these speleothems. A detailed analysis of this observation feature is beyond the scope of this paper, but in the context of a timescale perspective, which is the focus of this work, it highlights the caveats of climate wiggle matching between single records, even if the mechanistic link between regional climate changes may be relatively well understood.

Conclusion

We present the first radionuclide-based comparison between the Greenland Ice Core Chronology 2005 (GICC05) and the $\mathrm{U}/\mathrm{Th}$ timescale. We find that GICC05 is accurate within its stated absolute uncertainties, but also that the maximum counting error of the GICC05 may be at the limit to capture the total uncertainty accumulated within certain climatic periods. Our analysis indicates that the relationship between GICC05 and the $\mathrm{U}/\mathrm{Th}$ timescale over the last 45 ka drifts over time and reaches its maximum offset around 22 ka BP. We propose a transfer function that quantifies this drift and facilitates analysis of ice-core and $\mathrm{U}/\mathrm{Th}$ records, such as speleothems, on a common timescale. Thus, this transfer function allows for the further integration of key timescales in paleosciences and contributes to the INTIMATE (INTegration of Ice-core, MArine, and TEerrestrial records) initiative (Bjorck et al., 1996; Rasmussen et al., 2014b; Bronk Ramsey et al., 2014). Provided that $\mathrm{U}/\mathrm{Th}$ ages are regarded as accurate, the transfer function strongly reduces the absolute dating uncertainty of Greenland ice cores back to 45 ka BP. We reject the hypothesis of leads or lags larger than 189 years between Greenland, East Asia, and South America at the one sigma level. We show that the southward ITCZ shift around 24.5 ka BP seen in speleothems, typically associated with H2, coincides with the younger GS-3 dust peak recorded in Greenland ice cores. However, we also highlight inconsistencies between speleothem records around the onset of the Holocene, late GS-3, and GI-11 and thus caveats to the commonly applied practice of climate wiggle matching.

By comparing various ${}^{\mathrm{14}}\mathrm{C}$ records underlying IntCal13 as well as ice-core ${}^{\mathrm{10}}\mathrm{Be}$ data and geomagnetic field records, we infer that the current radiocarbon calibration curve underestimates the amplitude and rapidity of the $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ change around the Laschamp event at 41 ka BP. This adds to previous studies (Giaccio et al., 2017; Muscheler et al., 2014b) highlighting the fact that there are likely systematic errors in IntCal13 that will directly translate into errors of radiocarbon-based chronologies around that time. The combination of several internally inconsistent datasets in IntCal13 can lead to erroneous timing and amplitude of $\mathrm{\Delta }$${}^{\mathrm{14}}\mathrm{C}$ changes. Hence, great care has to be taken when attempting to use sections older than 30 ka BP of IntCal13 directly for studies of ${}^{\mathrm{14}}\mathrm{C}$ production rates and/or carbon cycle changes.

The transfer function shown in Fig. 12 will be made available as a Supplement to this paper and on NOAA.

The supplement related to this article is available online at: https://doi.org/10.5194/cp-14-1755-2018-supplement.

FA designed and carried out analyses and wrote the paper in correspondence with CBR and RM. TE designed and applied break-point detection analysis and wrote the corresponding “Methods and data” section (Sect. 3.5). FA, RM, CBR, SOR, CT, and AC initiated the project. RLE and HC provided speleothem data. AS and SOR provided insights into the ice-core chronology. HF and TE gave insights into aerosol transport and deposition. All authors discussed and commented on the paper.

The authors declare that they have no conflict of interest.

Acknowledgements

FA was supported through a grant by the Swedish Research Council (Vetenskapsrådet no. 2016-00218). CBR was partially supported through the UK Natural Environment Research Council (NERC) Radiocarbon Facility (NRCF010002). TE and HF acknowledge the long-term support of ice-core research at the University of Bern by the Swiss National Science Foundation (SNSF) and the Oeschger Center for Climate Change Research. SOR gratefully acknowledges support from the Carlsberg Foundation to the project ChronoClimate. This work was partially supported by the Swedish Research Council (grant DNR2013-8421 to RM), the NSF 1702816 to RLE, and the Australian Research Council DP170104665 to CT and AC. We gratefully acknowledge the financial support of the University of Adelaide Environment Institute for the initial Marble Hill meeting that initiated this work. Edited by: Denis-Didier Rousseau Reviewed by: Paula Reimer, Niklas Boers, Frédéric Parrenin, and Jeff Severinghaus