Introduction

Climate warming leads to environmental changes with consequent feedbacks on climate1,2. Feedbacks involving the biosphere are generally positive owing to the nonlinear stimulation of all biological processes by increasing temperature1,3. However, the magnitude of biosphere feedbacks on centennial timescales relevant to current global warming is poorly known3,4,5,6. Estimates of the strength of individual feedbacks based on modern observations (e.g. ref. 7) are hampered by the short length of the available records and uncertainties due to the influence of anthropogenic land-use change in recent decades. Earth System Models have been used to estimate the feedback strength8,9,10,11, but many biosphere processes are either not included or are poorly represented in the current generations of models12. Indeed, even when biosphere feedbacks are included, these modules are often not used in future projections or in simulations of the past.

Dansgaard-Oeschger (D-O) events are rapid climate fluctuations that occurred about 25 times during the last glacial period (ca 115 to 11.7 ka). They are characterised by a rapid warming over a few decades followed by a slower cooling over centuries to millennia13,14, with individual events registering warming of between 5 and 16 °C in Greenland15. This pattern is generally thought to reflect changes in the strength of the Atlantic Meridional Overturning Circulation (AMOC), whereby there is less poleward ocean heat transport when the AMOC is weak leading to cooling conditions around Greenland and vice versa16,17. The rapid warming events correspond to recovery of the AMOC. The cause of these events is still under debate and several mechanisms have been invoked, including ice-sheet instability18, sea-ice fluctuations linked to ice-shelf growth and decay19,20, sea-ice variability21,22, shifts in atmospheric circulation23,24 or in tropical climate modes24,25. The imprint of the D-O events is, nonetheless, reflected in large and globally synchronous changes in regional climates26,27,28 transmitted through the atmospheric circulation everywhere except Antarctica and surrounding regions, where the signal is dominated by a slower oceanic response to changes in the north29.

Ice-core records indicate that all of the D-O warmings were characterised by increased atmospheric CO2, CH4 and N2O concentrations30,31,32, showing that these events had an impact on global biogeochemical cycles4. While it has been suggested that the reinvigoration of the AMOC during D-O warming events could itself result in the physical release of CO2 to the atmosphere, diagnoses using a simple box model indicate the observed centennial-scale CO2 change is largely a result of carbon release due to the warming30. D-O events provide an opportunity to quantify the warming-induced greenhouse-gas feedbacks to climate on a centennial timescale relevant to contemporary climate change. Here, we exploit this opportunity to provide new estimates for CO2, CH4 and N2O climate feedbacks.

Feedback estimates from the Dansgaard-Oeschger events

The concept of feedback has been discussed in many previous studies, although terminologies differ2,3 (see Methods for quantitative explanations). To estimate feedback strengths in terms of the associated change in radiative forcing (W m−2) per degree (K) of global mean temperature change, we (a) identified the concentration changes in greenhouse gases from ice-core records across D-O events and converted them to radiative forcing; (b) used LOVECLIM model outputs to obtain the global mean temperature change during D-O events between 50 and 30 ka; and (c) combined both to derive feedback strengths, on the assumption that, on this timescale, the increase in global mean temperature leads to the increase in greenhouse gases.

Ice-core records of the concentration of CO2 (ref. 30), CH4 (ref. 31) and N2O (ref. 32) during the period between 50 and 30 ka (Fig. 1, Supplementary Table 1) were converted to a common timescale (AICC2012) based on the age-depth relationships for each chronology33. We estimated the change in CO2, CH4 and N2O concentration associated with the warming phase of each D-O event (Supplementary Fig. 1.1–1.8), using the dating of the beginning of these events from ref. 14, which was also converted to the AICC2012 timescale. The concentration changes of the three greenhouse gases were converted to radiative forcing using equations given in ref. 34, as adopted by IPCC WG1 AR6 (ref. 35), with concentration measurement uncertainties propagated into the corresponding radiative forcing uncertainties.

Fig. 1: Ice-core records between 50 and 30 ka.
figure 1

Changes in (a) CO2, (b) CH4, (c) N2O, (d) Greenland temperature, and (e) δD excess. The vertical lines show the official start dates of the numbered D-O warming events. All data are on the AICC2012 timescale (BP 1950).

Full size image

There are too few quantitative reconstructions of temperature changes, especially over land, to be able to make reliable estimates of changes in global mean temperature during the D-O events. We therefore use model-based estimates of the change in global mean temperature. The LOVECLIM model provides a global simulation of temperature changes during the interval 50–30 ka (ref. 36) in response to realistic time-varying changes in orbital parameters, atmospheric trace gas concentrations and ice-sheet configuration, and by adding meltwater pulses at the correct times required to trigger each D-O event. Evaluation of the experiments against individual records36,37 as well as comparison with the global compilation of palaeoclimate data in ref. 38 shows that it simulates the pattern of regional changes during individual D-O events during Marine Isotope Stage 3 well (Supplementary Fig. 2.1 to 2.8). We derived global mean temperature change by area-weighted averaging of the 64 × 32 grid cells, using the cosine of latitude as a weight (Fig. 2). The change in global mean temperature was identified in the same way as greenhouse gases (Supplementary Fig. 1.1–1.8).

Fig. 2: Global mean temperature anomaly to 30 ka.
figure 2

The data were obtained from LOVECLIM simulations and binned in 25 years. The global mean temperature was area-weighted, using the cosine of latitude as a weight for each grid. The age is at absolute timescale. The vertical lines show the official start dates of the numbered D-O warming events on AICC2012 timescale (BP 1950).

Full size image

The D-O events are not characterised by the ubiquitous warming of recent decades39 since, although most of the land was warming, the ocean warmed in the northern hemisphere and cooled in the southern hemisphere (Supplementary Fig. 2.1 to 2.8). Nevertheless, overall both ocean and land temperatures increased on average (Supplementary Fig. 2.9) and the land/ocean warming ratio was 1.48 ± 0.08 (95 % CI), comparable to present-day warming40,41. The amplitude and rate of global mean temperature increase (Supplementary Table 2) were also comparable to those of present day, which is 0.95–1.20 K increase by the decade 2011 ~2020 compared to pre-industrial times (1850–1900) with a rate of 0.0068–0.0085 K/year39. These similarities mean that D-O events usefully constrain present-day greenhouse-gas climate feedbacks.

The value of feedback strength (in unit of W m−2 K−1) is the ratio of the radiative forcing brought about by the increases in CO2, CH4 and N2O to the increase in global mean temperature during D-O events (Fig. 3). A maximum likelihood method42 is used to derive this ratio because it considers uncertainty of both the x- and y-variables (i.e. the driver and the response), in contrast with ordinary least squares regression which assigns uncertainty only to the y-variable. Based on the 8 D-O events that occurred between 50 and 30 ka, we estimated a feedback strength of 0.155 ± 0.035 W m−2 K−1 for CO2, 0.114 ± 0.013 W m−2 K−1 for CH4 and 0.106 ± 0.026 W m−2 K−1 for N2O (Table 1).

Fig. 3: Maximum likelihood estimation of feedback strengths.
figure 3

The figure shows the relationship between the increase in global mean temperature and radiative forcing induced by changes in (a) CO2, (b) CH4, (c) N2O concentrations and (d) combined radiative forcing of CO2, CH4 and N2O. Each D-O event is numbered; the horizontal and vertical lines show the 95 % confidence intervals. The measurements of CH4 concentration are very accurate so the vertical lines are too small to be observable on these plots. The solid line shows the maximum likelihood estimation of the ratio of radiative forcing to global mean temperature increase, the dashed lines show the 95 % confidence intervals of the ratio.

Full size image
Table 1 Feedbacks estimated from D-O events, with 95% confidence intervals.
Full size table

For comparison purposes, we adopt the dimensionless quantity ‘gain’, a measure of the extent to which the change in global mean temperature would be reduced (if gain is positive) or increased (if gain is negative) in the absence of the feedback. Gains are estimated by multiplying the feedback strengths (W m−2 K−1) by the climate sensitivity parameter (K W−1 m2). Climate sensitivity is defined as the equilibrium temperature increase of the Earth’s surface due to a radiative forcing (3.7 W m−2) equal to doubling atmospheric CO2 concentration compared to the pre-industrial level, after all the fast physical climate feedbacks (but not ice sheets and greenhouse-gas concentrations) are taken into account. Recent estimates of this equilibrium climate sensitivity (ECS), using different lines of evidence, yield ranges of 2.2–3.4 K (66% CI)43, 2.3–4.7 K (95% CI)44 and 2.4–4.5 K (95% CI)45. Assuming these estimates are independent, we derive an ECS of 3.23 ± 0.66 K (95% CI) and thus a climate sensitivity parameter (λ0) of 0.87 ± 0.18 K W−1 m2 (95% CI), which yields gains of 0.135 ± 0.041, 0.100 ± 0.023, 0.093 ± 0.029 for CO2, CH4 and N2O, respectively (Table 1). A larger number of estimates of climate sensitivity is given in IPCC WG1 AR6 (ref. 35). Assuming that the combined likely range of these estimates can be treated as equivalent to a 95% confidence interval, we derive an ECS of 3.0 ± 0.5 K and a value of λ0 of 0.81 ± 0.14 K W−1 m2, which yields gains of 0.125 ± 0.035, 0.093 ± 0.019, 0.086 ± 0.025 for CO2, CH4 and N2O, respectively (Table 1).

Comparison with previous estimates

Model-based feedback estimates have been derived from simulations of the response to anthropogenic emissions, and separate the carbon-concentration feedback and the carbon-climate feedback8. Changes in the atmospheric carbon concentration caused by emissions are buffered by the land and ocean uptake through the carbon-concentration feedback (a negative feedback); the amount of carbon these sinks can absorb is reduced by the carbon-climate feedback (a positive feedback)8. In the present-day context, anthropogenic CO2 emissions are the main driver of changes in the carbon cycle and warming is the response of the emissions; in the D-O context, warming is the main driver and changes in the atmospheric CO2 is the response. The feedback we quantified using D-O warmings equals the carbon-climate feedback defined in ref. 8. See Supplementary Notes for a more detailed explanation.

Model estimates of the carbon-climate feedback based on simulations from the Coupled Climate Carbon Cycle Model Intercomparison (C4MIP)8 show considerable variability, with estimates of the gain ranging from 0.04 to 0.31 (Fig. 4; see Supplementary Table 3 and Supplementary Fig. 3 for comparison of feedback strengths). The range is somewhat reduced in models from the fifth and sixth phases of the Coupled Model Intercomparison Project (CMIP510, CMIP646): 0.03 to 0.18 in CMIP5 and −0.002 to 0.18 in CMIP6. Our estimate of the CO2 gain derived from the D-O warming events is not consistent with high-end estimates from C4MIP, nor with low-end estimates from C4MIP, CMIP5 and CMIP6.

Fig. 4: Comparison of gains.
figure 4

Gains of this paper and previous estimates for (a) CO2, (b) CH4 and (c) N2O. Horizontal lines show the 95 % confidence intervals on each estimate. The shaded bars show the estimates from this paper, using 3.23 ± 0.66 K as the ECS. All the models estimate feedbacks at 2100. Stocker et al.11 only simulates the land climate feedbacks. The recalculated estimates for the Little Ice Age are based on the full 7000-member ensemble of global mean temperature reconstructions provided by the PAGES2k Consortium, using the 95% range to approximate 95% confidence intervals.

Full size image

There are relatively few model-based estimates of the feedbacks associated with either CH4 or N2O (Fig. 4). IPCC AR6 (ref. 47) estimated the CH4-climate feedback due to the effect of temperature on methanogenesis in wetlands as 0.03 ± 0.01 W m–2 K–1 (1 standard deviation, based on limited evidence) and an additional, highly uncertain feedback of 0.01 (0.003 to 0.04, 5th to 95th percentile range, also based on limited evidence) W m–2 K–1 due to permafrost thaw. Our results suggest that the CH4-climate feedback is larger than that assessed by AR6. Xu-Ri et al.9 simulated terrestrial N2O feedback estimate to be 0.11 W m−2 K−1, which corresponds to a gain of 0.10 ± 0.02. This is within the range estimated from the D-O warming events. Stocker et al.11 estimated the terrestrial feedbacks associated with CO2, CH4 and N2O to be 0.079, 0.011 and 0.023 W m–2 K–1 using the LPX-Bern vegetation model. IPCC AR6 (ref. 47) estimated the land N2O feedback as 0.02 ± 0.01 W m–2 K–1 (with low confidence) and the oceanic N2O feedback as −0.008 ± 0.002 W m–2 K–1 (based on limited evidence). Thus, AR6 indicates that the total N2O feedback is positive and dominated by the land, while the ocean feedback is smaller and of opposite sign. The combined (land plus ocean) feedback strength for N2O according to AR6 ((0.02 − 0.008) ± √(0.012 + 0.0022) = 0.012 ± 0.010 W m–2 K–1) however is considerably smaller than the value indicated by the D-O records.

Modern observations have been used to constrain model-based estimates of biosphere feedbacks. Gedney et al.7 used multi-site flux measurements as a constraint on simulated wetland CH4 emissions to obtain feedback estimates in the range of 0.01 to 0.11 W m−2 K−1 (Fig. 4). Other studies have used the emergent constraint approach to estimate the sensitivity of tropical land carbon storage to warming48,49, but only address part of the CO2 feedback and cannot be used to derive estimates of the gain. This lack of strong observational constraints has motivated the use of past climate changes to estimate greenhouse-gas feedbacks to climate50,51,52.

Previous attempts to quantify greenhouse-gas feedbacks using past climate changes have focused on the volcanically forced cooling during the Little Ice Age (LIA: 1500-1750 CE) which was associated with a decrease in CO2 of ca 8 ppm53. However, these estimates vary considerably and have high uncertainties (Fig. 4), in part associated with the temperature reconstruction used and in part due to differences in methodology. Scheffer et al.51 used alternative reconstructions of the LIA temperature change, derived from Mann and Jones54 and Moberg et al.55, and obtained estimates of 1/(1 − gCO2) of 1.28–2.93 and 1.07–1.25 (corresponding to a gain of 0.22–0.66 and 0.07–0.20, respectively). However, Cox and Jones52 obtained a much larger estimate of 40 ± 20 ppm CO2 per K, corresponding to a gain of 0.46 ± 0.25, using the Moberg et al. reconstruction55. Our recalculation of the CO2 feedback using the full 7000-member ensemble of temperature reconstructions provided by the PAGES2k Consortium56 produced a lower estimate of the gain than either Scheffer et al.51 using the Mann and Jones reconstruction54 or Cox and Jones52, but still with very large uncertainties that encompass almost all of the previous LIA estimates (Fig. 4). This uncertainty is also seen in recalculations of the gain associated with changes in CH4 and N2O over the LIA, suggesting that the LIA does not provide a sufficiently strong constraint to provide reliable estimates. In contrast, the D-O warming events provide a strong constraint because the temperature changes, and the responses, are relatively large. Furthermore, replication over 8 events considerably reduces uncertainty compared to using a single event such as the LIA.

We rely on the LOVECLIM model to derive estimates of global temperature because there are insufficient observationally based, quantitative reconstructions to estimate these reliably. Although a number of modelling groups have made simulations to mimic D-O events during the glacial by adding freshwater forcing57,58,59,60,61, none of these have used realistic forcings for individual D-O events. Comparison of the spatial patterns of the LOVECLIM simulated temperature changes for individual D-O events with records from the Voelker data compilation38 (Supplementary Fig. 2.1–2.8) indicate that there is good qualitative agreement in the sign of the change, with >75% of the grid cells being correctly predicted (Supplementary Table 4). Although LOVECLIM is a low-resolution model and the simulations were made with fixed cloud cover, neither of these constraints should have a major impact on the estimates of global temperature62. Furthermore, analyses based on estimating global temperature from observed temperature changes in Greenland over the interval between 80 and 20 ka using the relationship between simulated Greenland and global temperature obtained from the LOVECLIM simulations (see Supplementary Notes) produce comparable estimates of feedback strength. Thus, although the use of model outputs is a potential source of additional uncertainty, in the absence of a compelling alternative this approach provides a way to estimate greenhouse-gas climate feedbacks on centennial scales.

We have assumed that there is a strong relationship between global temperature changes and greenhouse-gas emissions during D-O warming events in order to estimate the climate feedback. Some of the changes in emissions may reflect change in hydroclimate, particularly in tropical regions27, but we presume that such changes are also conditioned by changes in temperature and thus will be reflected in the global temperature record. Similarly, changes in the balance of marine versus terrestrial sources of greenhouse-gas emissions, particularly CO2, are influenced by the changes in global temperature. There is currently insufficient information to disentangle the regional sources of greenhouse-gas emissions during the D-O events. However, the global feedback estimates obtained from analysis of the D-O events indicates that these feedbacks are non-negligible and poorly represented in current models.

Conclusions

We have used D-O cycles to estimate the climate feedbacks associated with CO2, CH4 and N2O. Based on the ECS estimate of 3.23 ± 0.66 K (3.0 ± 0.5 K), these feedbacks would amplify the equilibrium global mean temperature increase by 10–21% (10–19%), 8–14% (8–13%) and 7–14% (6–13%), respectively. The combined feedback from changes in all three greenhouse gases is 32–66% (31–56%). This finding indicates the importance of including greenhouse-gas feedbacks explicitly in climate change predictions. Although broadly compatible with previous estimates, our estimates have smaller uncertainties than previous observationally constrained estimates and provide a stronger constraint on model-based estimates.

Methods

Quantitative explanation of feedback terms

The concept of feedback has been explained quantitatively in many previous studies, although terminologies differ2,3,63,64. Briefly, a perturbation to the energy balance of a system, termed radiative forcing, pushes the system to a new equilibrium state with a change in temperature2,3,63. A reference system without feedbacks, gives a temperature increase (ΔT0) with a radiative forcing (ΔR0) when it reaches equilibrium; the ratio of ΔT0 to ΔR0, denoted λ0, is the climate sensitivity parameter of this reference system2,3,63.

$$\begin{array}{c}\triangle {T}_{0}\,=\,{\lambda }_{0}\triangle {R}_{0}\end{array}$$
(1.1)

Feedbacks results in additional radiative forcing. The temperature increase at equilibrium with feedback (ΔT) is thus:

$$\begin{array}{c}\triangle T\,=\,{\lambda }_{0}\left(\triangle {R}_{0}\,+\,\triangle {R}_{1}\,+\,\triangle {R}_{2}\,+\ldots\,+\,\triangle {R}_{n}\right)\end{array}$$
(1.2)

Assuming ΔR1, ΔR2, …, ΔRn proportional to ΔT with parameters c1, c2, …, cn (refs. 2,3) gives:

$$\begin{array}{c}\triangle T\,=\,{\lambda }_{0}\left(\triangle {R}_{0}\,+\,{c}_{1}\triangle T\,+\,{c}_{2}\triangle T\,+\,\ldots \,+\,{c}_{n}\triangle T\right)\end{array}$$
(1.3)

Combining Eqs. 1.1, 1.3 gives:

$$\begin{array}{c}\triangle T\,=\,\triangle {T}_{0}\,+\,{{\lambda }_{0}c}_{1}\triangle T\,+\,{{\lambda }_{0}c}_{2}\triangle T\,+\,\ldots \,+\,{\lambda }_{0}{c}_{n}\triangle T\end{array}$$
(1.4)

The terms c1, c2, …, cn express feedbacks in radiative forcing per degree of temperature increase (W m−2 K−1). This metric can be converted to a dimensionless measure called gain (g1, g2, …, gn), by multiplying by λ0 (refs. 2,3):

$$\begin{array}{c}\triangle T\,=\,\triangle {T}_{0}\,+\,{g}_{1}\triangle T\,+\,{g}_{2}\triangle T\,+\,\ldots \,+\,{g}_{n}\triangle T\end{array}$$
(1.5)

The relationship between the equilibrium temperature increase with and without feedback is thus:

$$\begin{array}{c}\triangle T\,=\,\frac{\triangle {T}_{0}}{\left(1\,-\,{g}_{1}\,-\,{g}_{2}\,-\,\ldots \,-\,{g}_{n}\right)}\,=\,\frac{\triangle {T}_{0}}{\left(1\,-\,{g}_{{total}}\right)}\end{array}$$
(1.6)

Equation 1.6 shows that: (a) a positive gain amplifies ΔT0 and a negative gain diminishes ΔT0; (b) the gain shows by what fraction ∆T0 is less than ∆T; a gain of 0.2, for example, means that ∆T0 is 20% less than ∆T, which means that ∆T is 25% more than ∆T0; (c) independent gains sum to gtotal, but their impacts on amplifying ΔT0 are not additive; two gains of 0.2, for example, combine to make ∆T 67% more than ∆T0 (refs. 2,3).

Ice-core data

We used the ice-core records of atmospheric CO2, CH4 and N2O concentrations detailed in Supplementary Table 1. The age models were converted to the Antarctic Ice-Core Chronology 2012 (AICC2012) timescale33 prior to analysis. The average resolution of the records on the AICC2012 timescale is 134 years for CO2, 18 years for CH4, and 59 years for N2O over the period 50–30 ka. Greenland temperatures were taken from the NGRIP ice core65, and again the original chronology was converted to the AICC2012 chronology before analysis. The average resolution for Greenland temperature is 19 years. We compare this with the δD excess record from the EPICA Dome C (EDC) ice core66, which has an average resolution of 49 years. Strictly speaking δD excess is interpreted as temperature changes in the source area rather than global temperature67.Nevertheless, it does clearly show the temperature changes associated with the D-O events.

The conversion to the AICC2012 chronology introduces additional uncertainties to those inherent in the original ice-core age models, particularly for the earlier part of the record33. Nevertheless, these are unlikely to have a remarkable effect given the method of determining the minimum and maximum response and thus of estimating the amplitude of change.

LOVECLIM temperature simulations

We used a transient simulation of the interval 50–30 ka performed with the LOVECLIM model36 to obtain an estimate of global mean temperature during the D-O events. LOVECLIM is a computationally efficient low-resolution (horizontal resolution of the atmospheric model is 5.625°) global climate model. The model was spun-up to equilibrium using an initial atmospheric CO2 concentration of 207.5 ppm, orbital forcing appropriate for 50 ka BP and an estimate of the 50 ka BP ice-sheet orography and albedo obtained from an off-line ice-sheet model simulation68. The transient run was initialised from this spin-up and run from 50 to 30 ka. Orbital, greenhouse gas, and ice-sheet forcings were updated continuously during the transient simulation; orbital parameters were derived following ref. 69, greenhouse-gas concentrations were from ice-core records, and the ice-sheet was from the off-line ice-sheet model simulation. In order to trigger D-O events, a time-series of freshwater inputs was derived by optimising freshwater fluxes such that simulated sea-surface temperature (SST) in the eastern subtropical North Atlantic were congruent with alkenone-based reconstructions of SST in that region.

We used simulated atmospheric temperature from the LOVECLIM experiments. Analyses in the original paper36 indicate that the simulations reproduce broadscale features of climate change during the D-O cycles well, and there is a good match with quantitative estimates for specific D-O events (e.g. D-O 8) from the Iberian margin and western Mediterranean regions, where highly resolved SST records are available36,37. The simulated air temperature changes over Greenland are somewhat smaller than those inferred from Greenland ice-core records36.

The LOVECLIM simulations were run with fixed cloud cover in these hindcast experiments. Studies examining the impact of using fixed clouds, albeit in a different model62, suggest that changes in cloud cover accentuate the temperature changes: it gets colder in the Northern Hemisphere, particularly in the North Atlantic region, but warmer in the Southern Hemisphere. However, the enhanced Northern Hemisphere cooling and Southern Hemisphere warming compensate each other so that the impact on global mean temperature is small. We assume that the same would be true in the LOVECLIM simulations.

To further evaluate the reliability of the LOVECLIM simulations, we compared the simulated temperature changes to reconstructions from the Voelker data set38, the only global data set that currently exists for MIS 3. Since this data set only contains a few records with quantitative estimates at high enough resolution to identify the temperature change for each D-O event, we compared the geographic patterns in warming or cooling trends globally (Supplementary Fig. 2.1–2.8; Supplementary Table 4). This analysis shows that: (a) the D-O events are registered as warmings over nearly all of the land areas in the world; (b) the geographic patterns of warming or cooling trends are consistent between observations and simulations, accepting that there may not be an exact geographic mapping because of the low resolution of the model; and (c) where there is quantitative information, the results are broadly consistent with the magnitude of the simulated changes.

Menviel et al.37 also showed that LOVECLIM surface air temperature and sea-surface temperature anomalies for D-O stadials and Heinrich stadials are consistent with observations, which provides further confidence that the LOVECLIM model captures global temperature change patterns linked to these events.

Identification of minimum and maximum

We used the start date of each D-O event provided by ref. 14. We then calculated binned averages of the CO2, CH4, N2O records and LOVECLIM simulated global mean temperature anomaly to 30 ka, centred on each D-O start date, using 25-year bins.

There is some uncertainty in the chronology of the start dates of each D-O event, and further uncertainty may be caused by the conversion from the GICC05 to the AICC2012 timescale (Supplementary Table 1). We therefore used a 200-year interval before and after the assumed D-O start date to identify the minima for CO2, CH4, N2O and LOVECLIM simulated global mean temperature anomaly to 30 ka. We assumed the maxima must occur within 500 years for CO2, CH4 and LOVECLIM simulated global mean temperature anomaly to 30 ka, and 600 years for N2O. The different lengths of time considered reflect the time needed to reach equilibrium and are also influenced by the resolution of the records. See Supplementary Fig. 1.1–1.8 for details for each D-O event.

Calculation of radiative forcing and propagation of uncertainties

We calculated binned values of each gas as follows:

$$\begin{array}{c}{c}_{{gas}}\,=\,\frac{\mathop{\sum }\limits_{k\,=\,1}^{m}{c}_{{gas},k}}{m}\end{array}$$
(5.1)
$$\begin{array}{c}{\sigma }_{{c}_{{gas}}}\,=\,\sqrt{\frac{\mathop{\sum }\limits_{k\,=\,1}^{m}{\sigma }_{{c}_{{gas},k}}^{2}}{{m}^{2}}}\end{array}$$
(5.2)

where cgas,k denotes a value in the bin with its standard error σcgas,k; m denotes the total number of values in this bin; cgas denotes the average value in this bin with its propagated standard error σcgas.

We calculated the radiative forcing34 associated with the change between minimum and maximum values for each event, as follows:

CO2:

$$\begin{array}{c}{\Delta R}_{C}\,=\,({a}_{1}{\left(C\,-\,{C}_{0}\right)}^{2}\,+\,{b}_{1}\left|C\,-\,{C}_{0}\right|+{c}_{1}\bar{N}\,+\,5.36)\,\times\, {{{{{\rm{ln}}}}}}\left(\frac{C}{{C}_{0}}\right)\end{array}$$
(5.3)
$$\begin{array}{c}{\sigma }_{{\Delta R}_{C}}\,=\,\sqrt{{\left(\frac{\partial {\Delta R}_{C}}{\partial C}\right)}^{2}{\sigma }_{C}^{2}\,+\,{\left(\frac{\partial {\Delta R}_{C}}{\partial {C}_{0}}\right)}^{2}{\sigma }_{{C}_{0}}^{2}}\end{array}$$
(5.4)

where a1 = −2.4 × 10−7 W m−2 ppm−1, b1 = 7.2 × 10−4 W m−2 ppm−1, c1 = −2.1 × 10−4 W m−2 ppb−1

CH4:

$$\begin{array}{c}{\Delta R}_{M}\,=\,({a}_{2}\bar{M}\,+\,{b}_{2}\bar{N}\,+\,0.043)\,\times\, \left(\sqrt{M}\,-\,\sqrt{{M}_{0}}\right)\end{array}$$
(5.5)
$$\begin{array}{c}{\sigma }_{{\triangle R}_{M}}\,=\,\sqrt{{\left(\frac{\partial {\triangle R}_{M}}{\partial M}\right)}^{2}{\sigma }_{M}^{2}\,+\,{\left(\frac{\partial \triangle {R}_{M}}{\partial {M}_{0}}\right)}^{2}{\sigma }_{{M}_{0}}^{2}}\end{array}$$
(5.6)

where a2 = −1.3 × 10−6 W m−2 ppb−1, b2 = −8.2 × 10−6 W m−2 ppb−1

N2O:

$$\begin{array}{c}{\triangle R}_{N}\,=\,({a}_{3}\bar{C}\,+\,{b}_{3}\bar{N}\,+\,{c}_{3}\bar{M}\,+\,0.117)\,\times\, \left(\sqrt{N}\,-\,\sqrt{{N}_{0}}\right)\end{array}$$
(5.7)
$$\begin{array}{c}{\sigma }_{\triangle {R}_{N}}\,=\,\sqrt{{\left(\frac{\partial \triangle {R}_{N}}{\partial N}\right)}^{2}{\sigma }_{N}^{2}\,+\,{\left(\frac{\partial {\triangle R}_{N}}{\partial {N}_{0}}\right)}^{2}{\sigma }_{{N}_{0}}^{2}}\end{array}$$
(5.8)

where a3 = −8.0 × 10−6 W m−2 ppm−1, b3 = 4.2 × 10−6 W m−2 ppb−1, c3 = −4.9 × 10−6 W m−2 ppb−1 C, M, N denote the maximum values identified for CO2, CH4 and N2O, respectively; C0, M0, N0 denote the minimum values identified for CO2, CH4 and N2O, respectively; \(\bar{C}\), \(\bar{M}\), \(\bar{N}\) denote the mean values identified for CO2, CH4 and N2O, respectively; ∆RC, ∆RM, ∆RN denote the radiative forcing brought about by CO2, CH4 and N2O, with their corresponding standard errors, σ∆RC, σ∆RM, σ∆RN, respectively.

Calculation of temperature increase and propagation of uncertainties

LOVECLIM provides yearly outputs, we used the standard deviation in each 25-year bin to approximate the standard error of the binned average.

The global mean temperature and its standard error was calculated as follows:

$$\begin{array}{c}{T}_{{mean}\;{global}}\,=\,\frac{\sum {T}_{{each}}{w}_{{each}}}{\sum {w}_{{each}}}\end{array}$$
(6.1)
$$\begin{array}{c}{\sigma }_{{T}_{{mean}\;{global}}}\,=\,\frac{\sqrt{\sum {\sigma }_{{T}_{{each}}}^{2}{w}_{{each}}^{2}}}{\sum {w}_{{each}}}\end{array}$$
(6.2)

where the weight of each grid (weach) is the cosine value of the latitude (in radian) of that grid.

We converted the data to anomaly to 30 ka as in the original paper36:

$$\begin{array}{c}{T}_{{mean}\;{global}\;{anomaly}}\,=\,{T}_{{mean}\;{global}}\,-\,{T}_{{mean}\;{global},30\,{ka}}\end{array}$$
(6.3)
$$\begin{array}{c}{\sigma }_{{T}_{{mean}\;{global}\;{anomaly}}}\,=\,\sqrt{{\sigma }_{{T}_{{mean}\;{global}}}^{2}\,+\,{\sigma }_{{T}_{{mean}\;{global},30\,{ka}}}^{2}}\end{array}$$
(6.4)

The global mean temperature change and its standard error were calculated using the minimum and maximum identified for each D-O event:

$$\begin{array}{c}\triangle {T}_{{mean}\;{global}}\,=\,{T}_{{mean}\;{global}\;{anomaly},{\max }}\,-\,{T}_{{mean}\;{global}\;{anomaly},{\min }}\end{array}$$
(6.5)
$$\begin{array}{c}{\sigma }_{\triangle {T}_{{mean}\;{global}}}\,=\,\sqrt{{\sigma }_{{T}_{{mean}\;{global}\;{anomaly},{\max }}}^{2}\,+\,{\sigma }_{{T}_{{mean}\;{global}\;{anomaly},{\min }}}^{2}}\end{array}$$
(6.6)

Calculation of gain and propagation of uncertainties

$$\begin{array}{c}g\,=\,c{\lambda }_{0}\end{array}$$
(7.1)
$$\begin{array}{c}{\sigma }_{g}^{2}\,=\,\sqrt{{c}^{2}{\sigma }_{{\lambda }_{0}}^{2}\,+\,{\lambda }_{0}^{2}{\sigma }_{c}^{2}}\end{array}$$
(7.2)

where g is the gain; σg is the standard error of the gain; c is the maximum likelihood estimated slope from the Deming package70, using radiative forcing (∆RC, ∆RM, ∆RN, ∆RC + ∆RM + ∆RN) and temperature increase (∆Tmean global) with the standard error of radiative forcing (σ∆RC, σ∆RM, σ∆RN, √(σ2∆RC + σ2∆RM + σ2∆RN)) and the standard error of temperature increase (σ∆Tmean global) as inputs; σc is the standard error of c obtained using the Deming package70; λ0 is the climate sensitivity parameter; σλ0 is the standard error of λ0. We derive two estimates of the climate sensitivity parameter λ0: using an ECS of 3.23 ± 0.66 K yields a value of λ0 of 0.87 K W−1 m2, and a value of σλ0 of 0.09 K W−1 m2; using an ECS of 3.0 ± 0.5 K yields a value of λ0 of 0.81 K W−1 m2, and a value of σλ0 of 0.07 K W−1 m2.

Calculation of gains from previous estimates

Some of the previous estimates give gains directly8,10; some give the amplifications51, 1/(1 − gain), which can be converted to gains easily; some provide values of c (radiative forcing per degree)7,9,11, which can be converted to gains using Eqs. 7.1, 7.2; some give the CO2 concentration gradient52, which can be approximated to gains using Eqs. 8.1, 8.2.

$$\begin{array}{c}g\,\approx\, d\alpha \end{array}$$
(8.1)
$$\begin{array}{c}{\sigma }_{g}^{2}\,\approx\, \sqrt{{d}^{2}{\sigma }_{\alpha }^{2}\,+\,{\alpha }^{2}{\sigma }_{d}^{2}}\end{array}$$
(8.2)

where g is the gain; σg is the standard error of the gain; d is the gradient (ppm CO2/K); σd is the standard error of the gradient; α is the climate sensitivity parameter (in the unit of K/ppm CO2), this paper uses 0.0115 K/ppm CO2; σα is the standard error of α, this paper uses 0.0012 K/ppm CO2.

Some previous estimates give the minimum and maximum concentration of gases and northern hemisphere temperature change during Little Ice Age (LIA)50, which can be converted to corresponding radiative forcing and global mean temperature change, assuming the northern hemisphere temperature to be 2/3 of the global mean temperature. There is only one estimate for the LIA, so maximum likelihood estimation of the slope is not available. Instead, we use Eqs. 8.3, 8.4 to derive c values, then use Eqs. 7.1, 7.2 to derive gains.

$$\begin{array}{c}c\,=\,\triangle R/\triangle T\end{array}$$
(8.3)
$$\begin{array}{c}{\sigma }_{c}\,=\,c\sqrt{{\left(\frac{{\sigma }_{\triangle R}}{\triangle R}\right)}^{2}\,+\,{\left(\frac{{\sigma }_{\triangle T}}{\triangle T}\right)}^{2}}\end{array}$$
(8.4)

where c is radiative forcing per kelvin; σc is the standard error of c; ∆R is the radiative forcing; σ∆R is the standard error of ∆R; ∆T is global mean temperature change; σ∆T is the standard error of ∆T.

Previous LIA estimates use either the Moberg et al.55 or Mann and Jones54 climate reconstructions; neither can now be assumed to be accurate. We therefore recalculated the feedbacks using the full 7000-member ensemble across all methods of the PAGES2k Consortium 2019 global mean temperature reconstructions56. Assuming the 95% range as an approximation of the 95% confidence interval, we derive a global mean temperature change (∆T) with a standard error (σ∆T). We identified the minimum and maximum concentration of the greenhouse gases with their standard errors from the same data source as in the LIA feedback papers, and converted these to radiative forcing using the same method as used for D-O events. Finally, we used Eqs. 8.3, 8.4 to obtain c values, then used Eqs. 7.1, 7.2 to obtain gains.