A 156 kyr smoothed history of the atmospheric greenhouse gases CO2, CH4, and N2O and their radiative forcing

Abstract

Continuous records of the atmospheric greenhouse gases (GHGs) CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O are necessary input data for transient climate simulations, and their associated radiative forcing represents important components in analyses of climate sensitivity and feedbacks. Since the available data from ice cores are discontinuous and partly ambiguous, a well-documented decision process during data compilation followed by some interpolating post-processing is necessary to obtain those desired time series. Here, we document our best possible data compilation of published ice core records and recent measurements on firn air and atmospheric samples spanning the interval from the penultimate glacial maximum ($\sim$ 156 kyr BP) to the beginning of the year 2016 CE. We use the most recent age scales for the ice core data and apply a smoothing spline method to translate the discrete and irregularly spaced data points into continuous time series. These splines are then used to compute the radiative forcing for each GHG using well-established, simple formulations. We compile only a Southern Hemisphere record of CH${}_{\mathrm{4}}$ and discuss how much larger a Northern Hemisphere or global CH${}_{\mathrm{4}}$ record might have been due to its interpolar difference. The uncertainties of the individual data points are considered in the spline procedure. Based on the given data resolution, time-dependent cutoff periods of the spline, defining the degree of smoothing, are prescribed, ranging from 5000 years for the less resolved older parts of the records to 4 years for the densely sampled recent years. The computed splines seamlessly describe the GHG evolution on orbital and millennial timescales for glacial and glacial–interglacial variations and on centennial and decadal timescales for anthropogenic times. Data connected with this paper, including raw data and final splines, are available at10.1594/PANGAEA.871273.

Introduction

Our knowledge of changes in the atmospheric mixing ratios of the important greenhouse gases (GHGs) CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O beyond the instrumental record is mainly based on discrete data points derived from gas extractions in polar ice cores. While there are recent developments towards continuous CH${}_{\mathrm{4}}$ records using gas extraction and measurement systems coupled to continuous-flow analysis systems , this approach has not yet been developed for the other two important GHGs, CO${}_{\mathrm{2}}$ and N${}_{\mathrm{2}}$O. To obtain the continuous GHG records, necessary for transient climate simulations, these discrete data have to be processed in order to extract those variabilities that have climatological significance and to account for measurement uncertainties.

All three GHG records have special features which need some attention during data compilation:

• For some of the CO${}_{\mathrm{2}}$ records obtained from different ice cores, there exist significant and as yet unexplained offsets . These offsets need to be addressed in our data compilation.

• Due to the dominance of CH${}_{\mathrm{4}}$ sources in the Northern Hemisphere, the CH${}_{\mathrm{4}}$ concentrations are higher in records from Greenland than from Antarctica referred to as interpolar difference; e.g..

• In situ production of N${}_{\mathrm{2}}$O connected to high mineral dust values leads to unreliable N${}_{\mathrm{2}}$O concentrations e.g., particularly during glacial peak times and in records from Greenland, for which special care has to be taken during data selection.

Rapid changes are most pronounced in CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O (and to some extent also in CO${}_{\mathrm{2}}$) during millennial-scale climate variability, or the so-called Dansgaard–Oeschger (D/O) events. Therefore, only well synchronised ice cores from Greenland and Antarctica can be used if records from the Northern and the Southern Hemisphere are to be merged into one global record. However, even with the recent efforts on ice core age scale development, there remain issues with this north–south synchronisation. For example, inconsistencies in the timing of abrupt changes in CH${}_{\mathrm{4}}$ concentration in the North Greenland Ice Core Project (NGRIP), EPICA Dronning Maud Land (EDML), and Talos Dome (TALDICE) ice cores have been identified for several D/O event transitions if based on AICC2012, the Antarctic Ice Core Chronology of four major Antarctic ice cores . Furthermore, when comparing data from the West Antarctic Ice Sheet Divide ice core (WDC) on its most recent age scale, WD2014, with data from Greenlandic ice cores, the chronology of the latter (GICC05) has been stretched by 0.63 % in order to find the best match to the absolute U/Th-dated paleo record of Hulu Cave .

In order for these issues to be overcome, careful data selection and processing are required. Here, we document our assumptions during data compilation and calculate continuous time series of CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O via spline-smoothing with a nominal temporal resolution $\mathrm{\Delta }t$ of 1 year from the penultimate glacial maximum until present, the time window of interest for PALMOD, the German Paleo Modelling Project (www.palmod.de). Note, however, that this $\mathrm{\Delta }t$ represents not the true resolution but only the typical spline average for each year and that the ice core information represents a low-pass filtered signal of atmospheric variability concentrations by the slow bubble enclosure process. Furthermore, the resulting spline is of restricted use for in-depth analysis with a focus on the rates of changes in the three GHGs, since the spline smoothing suppresses the most abrupt changes in the GHGs. Here, we extend the ice-core-based paleo records using instrumental data up until the beginning of the year 2016 CE, including several decades of overlap between the ice core and instrumental data. The resulting continuous GHG records might also be of interest and may be used in the Last Deglaciation experiment within PMIP4 (Paleoclimate Modelling Intercomparison Project phase 4) . Note that different GHG data sets have so far been chosen to force transient simulations for the last 21 kyr in , but well-motivated different set-ups (e.g. using the GHG splines compiled here) are possible within PMIP4.

Previous splines (similar to our approach here but not identical in detail) have also been proposed to be used in interglacial experiments of the Holocene within PMIP4 . Within the most recent model intercomparison project, the Coupled Model Intercomparison Project Phase 6 (CMIP6), a slightly different compilation of GHGs for historical times, or the Common Era, has been presented . While this alternative approach has its focus on the time since 1850 CE, its data compilation nevertheless extends back until the year 0 CE, based solely on the Law Dome ice core in non-instrumental times . We will finally compare our splines with these forcing data sets proposed by to be used within CMIP6.

As will be seen in detail in the next section, the mathematical formulation of the spline smoothing method needs information on the uncertainties or errors in the data points supporting the spline. These data uncertainties represent the precisions of individual measurements ($\mathrm{1}\mathit{\sigma }$ errors) and are of the order of a few parts per million for CO${}_{\mathrm{2}}$ or a few parts per billion for CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O. The uncertainty in the final spline, however, is larger, since the applied smoothing, which depends on the chosen cutoff periods, adds some additional uncertainty. Furthermore, the estimates of the radiative forcing based on these three GHGs given here are even more uncertain, since the calculations of the radiative forcing themselves are based on models with an embedded intrinsic uncertainty of about $\sim$ 10 % . Note that the calculations of the GHG radiative forcing provided here are just a first-order approximation, since we use the simplified expressions of , while full climate models calculate radiative forcing internally, when forced with variable GHG concentrations.

In the following, ages are either given in years CE (Common Era) or in years BP (before present), where present is defined as 1950 CE. We define the onset of anthropogenic activities at 1750 CE (or 200 BP), based on the timing of the increase in CO${}_{\mathrm{2}}$ and CH${}_{\mathrm{4}}$ in our final splines, although we acknowledge that the onset of the Anthropocene is still debated e.g..

Details on the spline smoothing method

The numerical code for spline smoothing is based on , but see also and for further details, discussions, and applications. It offers the possibility to select different cutoff periods for different time intervals or parts of the input data set, which is needed when data spacing is variable throughout the data set.

In a smoothing spline a cost function is minimised. This cost function includes two terms: (i) the error-weighted deviation between the spline value and the actual data value and (ii) the curvature of the spline, represented by its second derivative. A parameter $\mathit{\lambda }$ defines how much weight is given to the curvature. For a large $\mathit{\lambda }$, the optimisation results in low curvature, i.e. a very smooth spline and relatively large deviations from the original data. Similarly, increasing errors in the data results in a smoother spline for a given $\mathit{\lambda }$. In other words, the smoothing of the spline depends on both the assumed errors in the data and the parameter $\mathit{\lambda }$.

According to Fourier, each time series can be represented by a sum of sine functions. Since a smoothing spline acts as a low-pass filter, high frequencies are dampened in the spline. The period at which the amplitude is attenuated to 50 % is defined as the cutoff period $P_{\mathrm{c}}$ e.g.. The parameter $\mathit{\lambda }$ is linked to $P_{\mathrm{c}}$ as described in detail in Eq. () below.

Let us assume input data are $t_j$, $y_j$, and $v_j$ corresponding to time, value, and error ($\mathrm{1}\mathit{\sigma }$). For a given interval of the input data, an average error, $v$, and an average data spacing, $\mathrm{\Delta }t$, can be computed. The link between the cutoff period ($P_{\mathrm{c}}$), the data spacing ($\mathrm{\Delta }t$), and the $\mathrm{1}\mathit{\sigma }$ error in the input data ($v$) is $$P_{\mathrm{c}}=\mathrm{2}\mathit{\pi }\cdot {\left(\mathit{\lambda },\cdot ,\mathrm{\Delta },t,\cdot ,v^{\mathrm{2}}\right)}^{\mathrm{0.25}}.$$

In the following, we prescribe $P_{\mathrm{c}}$ and can calculate $\mathit{\lambda }$ following the given relationship in Eq. (). We choose a $P_{\mathrm{c}}$ value such that it is much larger than the temporal resolution of the data, $\mathrm{\Delta }t$, to avoid overfitting. However, since the choice of $P_{\mathrm{c}}$ is also partially subjective, we investigate its influence on the final spline by sensitivity studies, in which $P_{\mathrm{c}}$ is varied by $\pm$50 %. One aspect of Eq. () is that $P_{\mathrm{c}}$ depends only weakly on $\mathrm{\Delta }t$.

Locations of the different data sources, ordered north to south. Individual sites of the NOAA observational network are not explicitly mentioned here, when they only contribute to global mean calculations. SH CH${}_{\mathrm{4}}$: Southern Hemisphere CH${}_{\mathrm{4}}$.

Site	Latitude	Longitude	Data used here
NGRIP	75.10${}^{\circ }$ N	42.32${}^{\circ }$ W	N${}_{\mathrm{2}}$O
GRIP	72.583${}^{\circ }$ N	37.633${}^{\circ }$ W	comparing to SH CH${}_{\mathrm{4}}$
Barrow	71.3230${}^{\circ }$ N	156.6114${}^{\circ }$ W	comparing to SH CH${}_{\mathrm{4}}$
Mauna Loa	19.5362${}^{\circ }$ N	155.5763${}^{\circ }$ W	CO${}_{\mathrm{2}}$
Law Dome${}^{\mathrm{1}}$	$\sim$ 66.73${}^{\circ }$ S	$\sim$ 112.83${}^{\circ }$ E	CO${}_{\mathrm{2}}$, SH CH${}_{\mathrm{4}}$, N${}_{\mathrm{2}}$O
Talos Dome (TALDICE)	72.817${}^{\circ }$ S	159.183${}^{\circ }$ E	CO${}_{\mathrm{2}}$, N${}_{\mathrm{2}}$O
EPICA Dronning Maud Land (EDML)	75.0${}^{\circ }$ S	0.067${}^{\circ }$ E	CO${}_{\mathrm{2}}$
EPICA Dome C (EDC)	75.1${}^{\circ }$ S	123.35${}^{\circ }$ E	CO${}_{\mathrm{2}}$, SH CH${}_{\mathrm{4}}$, N${}_{\mathrm{2}}$O
Taylor Glacier${}^{\mathrm{2}}$	$\sim$ 77.77${}^{\circ }$ S	$\sim$ 161.7${}^{\circ }$ E	N${}_{\mathrm{2}}$O
WAIS Divide Ice Core (WDC)	79.468${}^{\circ }$ S	112.086${}^{\circ }$ W	CO${}_{\mathrm{2}}$, SH CH${}_{\mathrm{4}}$
Siple Dome	81.66${}^{\circ }$ S	148.82${}^{\circ }$ W	CO${}_{\mathrm{2}}$
South Pole${}^{\mathrm{1}}$	90${}^{\circ }$ S	59${}^{\circ }$ E	CO${}_{\mathrm{2}}$, SH CH${}_{\mathrm{4}}$

Notes: ${}^{\mathrm{1}}$ The data compilation of and on CO${}_{\mathrm{2}}$, SH CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O uses data from the Law Dome deep ice core and from various shallow ice and firn cores in its vicinity but also from atmospheric data from Cape Grim and firn core data from the South Pole. While we here state all the relevant positions, the original source of the individual data points is not marked in Tables , , and or in the data files uploaded to PANGAEA, where data are only labelled with “Law Dome” as their source. Please see the original references for further details. ${}^{\mathrm{2}}$ Data taken from Taylor Glacier are based on a “horizontal ice core”, which does not have a point location like all other sites do.

Let us now assume we have a data set with variable data spacing, for which we would like to apply different smoothing depending on $\mathrm{\Delta }t$. We proceed by modifying $\mathit{\lambda }$ to follow the predefined individual $P_{\mathrm{c}}$ for each interval of the input data set as follows.

• Reference interval: We take the most recent time window, consisting of instrumental measurements, as reference interval. $\mathit{\lambda }$ is computed using Eq. () for the given cutoff period, average data spacing, and average error for this first interval.

• Other intervals: A modified ${\mathit{\lambda }}^{\prime }=\mathit{\lambda }\cdot s^{\mathrm{2}}$, with $\mathit{\lambda }$ taken from the reference interval, is used for other intervals, implying that for the reference interval $s=\mathrm{1}$ and ${\mathit{\lambda }}^{\prime }=\mathit{\lambda }$. The scaling factor $s$ is chosen to gain the desired $P_{\mathrm{c}}^{\prime }$ after

  $$s={\displaystyle \frac{(\frac{P_{\mathrm{c}}^{\prime }}{\mathrm{2}\mathit{\pi }}{)}^{\mathrm{2}}}{\sqrt{(\mathit{\lambda }\cdot \mathrm{\Delta }{t}^{\prime })}\cdot v^{\prime }}},$$

  where $P_{\mathrm{c}}^{\prime }$, $\mathrm{\Delta }{t}^{\prime }$, and $v^{\prime }$ are the cutoff period, the mean data spacing, and the mean error for the interval under consideration.

An intermediate product with $t_j,y_j$, and $v_j^{\prime }$ is calculated, in which the revised uncertainty $v_j^{\prime }$ is defined by Eq. () using the cutoff-related scaling factor $s$. From this intermediate product, the final spline with time-dependent $P_{\mathrm{c}}$ is calculated. In doing so, the approach abstains from any further merging of partial time series to a final spline. The resulting spline follows the prescribed cutoff periods throughout the whole time series. However, for every change in cutoff period from $P_{\mathrm{c}\mathrm{1}}$ to $P_{\mathrm{c}\mathrm{2}}$ a transition window around the time of change, $t_{\mathrm{change}}$, exists (defined as $t_{\mathrm{change}}\pm P_{\mathrm{ct}}$, with $P_{\mathrm{ct}}$ being the smaller of $P_{\mathrm{c}\mathrm{1}}$ and $P_{\mathrm{c}\mathrm{2}}$), in which the variability of the spline transits from one cutoff period to the other and does not follow the prescribed $P_{\mathrm{c}}$ exactly. The summaries of the spline calculation contained in Tables , , and show the effect of this transition in a column of averaged realised cutoff periods, which are always slightly lower than the prescribed cutoff periods.

Data used to construct the CO${}_{\mathrm{2}}$ spline.

Time (in BP)	Time (in CE)	Source	Age scale	Citation
$-$66 to $-$8	2016 to 1958	Mauna Loa${}^{\mathrm{1}}$ (monthly)	–
$-$10 to 1949	1960 to 1	Law Dome${}^{\mathrm{2}}$	as in references.
200 to 1210	1750 to 740	WDC${}^{\mathrm{3}}$	WD2014
1902 to 10 954	48 to before CE	EDC${}^{\mathrm{4}}$	AICC2012
8807 to 22 909	–	WDC${}^{\mathrm{6}}$	WD2014
21 926 to 48 720	–	Siple Dome	GICC05
38 127 to 69 672	–	Talos Dome${}^{\mathrm{5}}$	AICC2012
43 205 to 113 429	–	EDML	AICC2012
104 331 to 156 306	–	EDC	AICC2012
124 859 to 153 135	–	EDC	AICC2012

Notes: ${}^{\mathrm{1}}$ Data taken from ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt. ${}^{\mathrm{2}}$ Law Dome data are taken from various sources; see references for details. They are available from 2001 CE to 1 CE; only data before 1960 CE are taken for the spline. ${}^{\mathrm{3}}$ WDC data are available from 10 BP to 1217 BP, but not all were used here. Published WDC data have been shifted by $-$3.13 ppm (see Fig. A1). ${}^{\mathrm{4}}$ EDC data are available from 350 BP to 22 236 BP, but not all were used here. ${}^{\mathrm{5}}$ Talos Dome data exist from 34 360 BP but contain some outliers before 38 kyr BP. ${}^{\mathrm{6}}$ Published WDC data have been shifted by $-$6.06 ppm (see Fig. A2).

The uncertainties of the final splines are calculated from the square root of the sum of squares of three individual errors ($\mathit{\sigma }=\sqrt{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{3}}^{\mathrm{2}}}$).

• $P_{\mathrm{c}}$ error (${\mathit{\sigma }}_{\mathrm{1}}$): Mean difference from the standard spline by smoothing with cutoff periods $P_{\mathrm{c}}$ which are varied by $\pm \mathrm{50}$ %.

• Data resolution error (${\mathit{\sigma }}_{\mathrm{2}}$): The importance of uncertainty of the individual data points $v_i$ for the spline smoothing by setting all $v_i$ to 0.01.

• Monte Carlo error (${\mathit{\sigma }}_{\mathrm{3}}$): Repeated ($n=\mathrm{500}$) realisation of the data sets $y_i$ by randomly choosing data points out of the normally distributed data using the given uncertainty ranges $v_i$.

Greenhouse gas data compilations and final splines

Our GHG data compilations are based on various data sets from 13 global distributed locations. An overview of the locations, including latitude and longitude, is provided in Table . Please note that CH${}_{\mathrm{4}}$ data are only included from Southern Hemisphere records. These pointwise data sets are supplemented for the instrumental period by some global mean data from the National Oceanic and Atmospheric Administration (NOAA) observational network, including Radiatively Important Trace Species (RITS) nitrous oxide data from the Earth System Research Laboratory (NOAA/ESRL) halocarbons program and nitrous oxide data from the NOAA/ESRL halocarbons in situ program, which consists of globally distributed measurements. Individual data uploaded to the database PANGAEA, based on and are all labelled as “Law Dome” data for simplicity, although these two studies contain data from the Law Dome deep ice core, data from various shallow cores, and atmospheric data from Cape Grim and the South Pole. Please refer to the original publications for a precise characterisation of the sample origins.

CO${}_{\mathrm{2}}$ spline covering all data: 2016 CE–156 307 BP. Error bars around the ice core data points are $\pm \mathrm{2}\mathit{\sigma }$. WDC data have been adjusted to reduce offsets; see text for details. In (a) the right axis contains the resulting radiative forcing $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{5.35}\cdot ln({\mathrm{CO}}_{\mathrm{2}}/(\mathrm{278}\mathrm{ppm}\left)\right)$ W m${}^{-\mathrm{2}}$ calculated after . (b) Total uncertainty of the spline based on three individual error sources; see text for details. (c) Temporal resolution ($\mathrm{\Delta }t$) of the CO${}_{\mathrm{2}}$ data points underlying the spline on a log scale. Additionally, the prescribed time-dependent cutoff period $P_{\mathrm{c}}$ is plotted, including its variation by $\pm$50 %, which has been used to determine ${\mathit{\sigma }}_{\mathrm{1}}$.

[Figure omitted. See PDF]

Statistics of the CO${}_{\mathrm{2}}$ spline. Interval: $i_{{\mathrm{CO}}_{\mathrm{2}}}$; $s$: scaling factor to fulfil the constraints given by the prescribed cutoff period $P_{\mathrm{c}}$; $\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$: average realised cutoff period; $\mathrm{\Delta }t$: mean data spacing; $v$: mean $\mathrm{1}\mathit{\sigma }$ error – exact time framing is given by the age of the first ($t_{\mathrm{start}}$) and last ($t_{\mathrm{stop}}$) data point of the interval (in years BP); $N$: number of data points within interval. In the last column the underlying data source is briefly mentioned; see Table  for details and citations.

$i_{{\mathrm{CO}}_{\mathrm{2}}}$	$s$	$P_{\mathrm{c}}$	$\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$	$\mathrm{\Delta }t$	$v$	$t_{\mathrm{start}}$	$t_{\mathrm{stop}}$	$N$	Data source
	–	yr	yr	yr	ppm	yr BP	yr BP	–
1	1.00	4.0	4.0	0.1	0.3	$-$66.0	$-$8.1	699	Mauna Loa, Law Dome
2	2.65	20.0	18.5	0.4	1.3	$-$8.0	19.2	69	Law Dome
3	7.49	40.0	37.5	1.0	1.2	20.6	117.1	96	Law Dome
4	79.91	160.0	151.0	4.3	0.9	123.1	997.8	206	Law Dome, WDC
5	388.87	500.0	468.8	13.0	1.0	1006.0	1796.5	62	Law Dome, WDC
6	5377.63	3000.0	2883.1	93.3	1.0	1814.0	8994.9	78	Law Dome, WDC, EDC
7	1532.28	1600.0	1519.4	48.8	1.3	9060.2	10 962.5	40	WDC, EDC
8	357.66	600.0	567.4	28.1	1.0	11 060.3	18 463.6	264	WDC
9	1563.85	2000.0	1806.4	176.2	1.1	18 559.8	109 840.0	519	WDC, Siple D, Talos Dome, EDML, EDC
10	1690.90	3000.0	2593.3	383.9	1.5	110 555.4	127 829.0	46	EDML, EDC
11	225.60	1000.0	921.5	257.3	1.5	128 024.5	134 970.7	28	EDC
12	530.27	2000.0	1853.3	871.7	1.4	135 387.0	156 306.8	25	EDC

Atmospheric CO${}_{\mathrm{2}}$

There are small offsets of a few parts per million in measured CO${}_{\mathrm{2}}$ concentration between records obtained from different ice cores e.g.. These offsets may be related to inter-laboratory differences in the calibration or potentially due to in situ artefact production of CO${}_{\mathrm{2}}$ in the ice archive. For a detailed discussion, see and the supplement to . In addition to these offsets, the amplitudes of GHG variations can differ from one core to the next due to the site-dependent bubble enclosure characteristics, which act as low-pass filtering. Offsets require the adjustment of individual records to avoid spurious CO${}_{\mathrm{2}}$ changes when linking different records from different laboratory and ice cores. Ice core data are considered here on the best (most recent) age model available, whose details are contained in Table . AICC2012 refers to the most recent Antarctic Ice Cores Chronology, providing age models for EPICA Dome C (EDC), EDML, Talos Dome, Vostok, and the NGRIP record . The CO${}_{\mathrm{2}}$ record from WDC is used here on its more recent age scale WD2014 to have the timing of CO${}_{\mathrm{2}}$ and the other two GHGs consistently on the same chronology. Using WD2014 instead of the original chronology WDC06A-7 shifts the WDC CO${}_{\mathrm{2}}$ time series towards younger ages: by about 100 years during Termination I and by about 10 years during the last 1.2 kyr.

Our CO${}_{\mathrm{2}}$ data compilation extends to $\sim$ 156 kyr BP, at which point in time well-resolved CO${}_{\mathrm{2}}$ records stop. The full CO${}_{\mathrm{2}}$ spline covering the whole time window from 2016 CE to 156 kyr BP is plotted in Fig. , including its uncertainty estimate (b) and the temporal resolution, $\mathrm{\Delta }t$, of the compilation of data points (c). The 11-point running mean of $\mathrm{\Delta }t$ is around 100 years in the Holocene, between 20 and 50 years during Termination I, varies between 40 and 200 years between 20 and 70 kyr BP, and rises to 1000 years prior to 70 kyr BP. Across Termination II, $\mathrm{\Delta }t$ decreases to an average of 200 years.

The CO${}_{\mathrm{2}}$ data contributing to this spline are described below (further details in Table ):

Comparison of our final spline data with values used for PMIP4 experiments for 21 kyr and 1850 CE, 6 kyr, and 127 kyr . Please note that the PMIP4 values should be millennial-scale mean numbers to serve as forcing values for time slice experiments, while the values given from our study are snapshots of the given points in time. Furthermore, we calculate SH CH${}_{\mathrm{4}}$ values, while in PMIP4 the global CH${}_{\mathrm{4}}$ is given.

GHG	Unit	1850 CE	6 kyr	21 kyr	127 kyr
Our study:
CO${}_{\mathrm{2}}$	ppm	286.1	264.4	187	274
SH CH${}_{\mathrm{4}}$	ppb	795	553	382	660
N${}_{\mathrm{2}}$O	ppb	271	261	206	257
PMIP4:
CO${}_{\mathrm{2}}$	ppm	284.6	264.4	190	275
global CH${}_{\mathrm{4}}$	ppb	808	597	375	685
N${}_{\mathrm{2}}$O	ppb	273	262	200	255

Details of the CO${}_{\mathrm{2}}$ spline. Light and dark grey bands around the spline represent $\pm \mathrm{1}\mathit{\sigma }$ and $\pm \mathrm{2}\mathit{\sigma }$, respectively; error bars around the ice core data points are $\pm \mathrm{2}\mathit{\sigma }$. (a) Instrumental times (1950–2016 CE); (b) 0–2000 BP; (c) Termination I; (d) 0–40 kyr BP without the Law Dome data showing the anthropogenic rise; (e) 40–90 kyr BP; (f) 90–160 kyr BP. WDC data have been adjusted to reduce offsets; see text for details. Dashed line labelled CMIP6 in panels (a) and (b) is the compiled CO${}_{\mathrm{2}}$ record to be used in CMIP6 experiments for the last 2 kyr .

[Figure omitted. See PDF]

• Our CO${}_{\mathrm{2}}$ data compilation uses instrumental monthly CO${}_{\mathrm{2}}$ data taken from the NOAA network up to the beginning of the year 2016 CE, or $-\mathrm{66.0}$ BP . We here choose to take only the data of the original so-called “Keeling curve” started by C.D. Keeling at NOAA's Mauna Loa Observatory in 1958 CE, and since 1974 CE independently measured by both the Scripps Institution of Oceanography (scrippsco2.ucsd.edu) and NOAA (www.esrl.noaa.gov/gmd/ccgg/trends/) (Fig. a). There is a small interpolar difference in CO${}_{\mathrm{2}}$ concentrations, with higher concentrations in the north than in the south; e.g. the 10-year averages from 2006 CE to 2015 CE were 3.5 ppm lower at the South Pole than at Mauna Loa and 1.4 ppm higher at Barrow (Alaska) than at Mauna Loa . We therefore assume that CO${}_{\mathrm{2}}$ data from Mauna Loa are a good approximation of the global average concentration. However, in practice this interpolar difference cannot be determined prior to the instrumental records since CO${}_{\mathrm{2}}$ is only measured on ice cores from Antarctica, as the higher impurity content gives rise to artefacts in any CO${}_{\mathrm{2}}$ measurement based on Greenlandic ice cores e.g..

• The firn and ice data compilation of Law Dome, which also contains some contributions from Cape Grim and the South Pole – available for the time from 1996 CE to 1 CE ($-\mathrm{46}$ BP to 1949 BP) and 2001 CE to 154 CE ($-\mathrm{51}$ BP to 1796 BP) – overlap consistently with direct atmospheric measurements. We therefore take these data as our reference time series for the Common Era (Fig. b) but include only the data from year 1960 CE and older in our spline compilation. In doing so, we use the more precise and temporally more highly resolved instrumental data for later times.

• Data from the WDC ice core exist for the times of 11–1210 BP, or 1939–740 CE and for Termination I (see point no. 5 below). These WDC data overlap with the Law Dome data ; however, the available high-resolution CO${}_{\mathrm{2}}$ records from different ice cores (Law Dome, WDC, EDML) show some apparent offsets . Whilst the CO${}_{\mathrm{2}}$ data in all three ice cores converge on similar concentrations during the anthropogenic rise in CO${}_{\mathrm{2}}$ after 1750 CE, the WDC CO${}_{\mathrm{2}}$ concentrations are slightly higher than CO${}_{\mathrm{2}}$ in the other two ice cores prior to 1750 CE. In pre-anthropogenic times, EDC data not contained in the comparison of also agree more with the Law Dome data than with those of WDC. We therefore choose to take WDC data only prior to the anthropogenic rise (200–1210 BP or 1750–740 CE). Furthermore, WDC data are adjusted by $-\mathrm{3.13}$ ppm to bring them into agreement with the Law Dome CO${}_{\mathrm{2}}$ record. The data from Law Dome and the adjusted data from WDC contribute to our data compilation between 200–1210 BP. The mean temporal resolutions of both ice core CO${}_{\mathrm{2}}$ records within this time interval are 8 and 13 years for WDC and Law Dome, respectively. The amplitude of the CO${}_{\mathrm{2}}$ minima around 300–400 BP is controversial . In our final spline, little of the large negative anomaly in CO${}_{\mathrm{2}}$ contained in the Law Dome data is preserved, since we smoothed the ice data in this time window with a cutoff period of 160 years (Fig. b). The time between the start of the anthropogenic rise (1750 CE) and the start of the instrumental record (1958 CE) is only supported by the Law Dome data in our compilation (Fig. b). Further details on this adjustment of the WDC data are covered in Fig. A1 in the Appendix.

• EDC data exist between 350 BP and the Last Glacial Maximum (LGM) and further back in time (see point no. 7 below). They overlap with the Law Dome data between 350 and 1950 BP (Fig. b) without any apparent offset, and therefore no adjustment is applied to the EDC data. However, EDC data are only included in our compilation for the interval 1.9–11 kyr BP because Law Dome and WDC data provide a better resolution for the interval younger than 1.9 kyr BP, whilst the WDC data are the more highly resolved record for the interval older than 11 kyr BP (Fig. c).

• Termination I is best covered by data from WDC . WDC data are available for 8.8–22.9 kyr BP and are adjusted by $-$6.06 ppm (Fig. c). This difference corresponds to the duration-weighted mean offset between the WDC and EDC records during three intervals of relatively constant CO${}_{\mathrm{2}}$ (22.3–18.5 kyr BP: WDC ($n=\mathrm{29}$) $\mathrm{194.75}\pm \mathrm{2.44}$ ppm; EDC ($n=\mathrm{21}$) $\mathrm{188.22}\pm \mathrm{2.32}$ ppm; 14.5–13.0 kyr BP: WDC ($n=\mathrm{45}$) $\mathrm{243.02}\pm \mathrm{2.44}$ ppm; EDC ($n=\mathrm{9}$) $\mathrm{237.57}\pm \mathrm{1.42}$ ppm; 11.5–9.0 kyr BP: WDC ($n=\mathrm{36}$) $\mathrm{269.97}\pm \mathrm{3.67}$ ppm; EDC ($n=\mathrm{27}$) $\mathrm{264.24}\pm \mathrm{1.88}$ ppm). The intervals have been selected to minimise the influence of potential age scale differences between the two records. Only those EDC studies focusing on CO${}_{\mathrm{2}}$ measurements have been considered here, rather than those with a main focus on ${\mathit{\delta }}^{\mathrm{13}}$C , which have a lower precision in CO${}_{\mathrm{2}}$ concentrations. More details on this adjustment of the WDC data during Termination I are found in Fig. A2. Our offset corrections imply an absolute CO${}_{\mathrm{2}}$ concentration uncertainty of about 5 ppm (accuracy). The corresponding uncertainty in the radiative forcing $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ following a simplified expression of ,

  $$\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{5.35}\cdot \ln ({\mathrm{CO}}_{\mathrm{2}}/(\mathrm{278}\mathrm{ppm}\left)\right)\mathrm{W}{\mathrm{m}}^{-\mathrm{2}},$$

  is 0.15 or 0.09 W m${}^{-\mathrm{2}}$ for a reference concentration of 180 or 280 ppm, respectively. This uncertainty is larger than the measurement uncertainty (precision) of the order of 1 ppm attached to individual data points which is used to determine the smoothing spline through the data.

• Further back in time all ice core records used have some data overlap with their successive records. There are some small offsets between the different records for details, see. We treat them all alike, so the spline averages over all cores, and we select a large cutoff period of 2000 years for the interval 18.5–110 kyr BP to account for those uncertainties. Rapid variations in CO${}_{\mathrm{2}}$ during glacial times (Fig. d–f) are best recorded in the Siple Dome record between 21.9 and 48.7 kyr BP , the Talos Dome record between 34.4 and 69.7 kyr BP , and the EDML record between 43.2 and 113.4 kyr BP . Talos Dome CO${}_{\mathrm{2}}$ data include some outliers in the interval 34–38 kyr BP that disagree with CO${}_{\mathrm{2}}$ records from other ice cores by more than 10 ppm. Therefore, Talos Dome data are only considered for the times older than 38.0 kyr BP.

• From 104.3 to 156.3 kyr BP – the interval spanning the last glacial inception, the last interglacial, Termination II, and the penultimate glacial maximum (Fig. f) – the EDC CO${}_{\mathrm{2}}$ record is used in the compilation .

For every supporting data point $j$ a $\mathrm{1}\mathit{\sigma }$ uncertainty or error $v_j$ has to be assigned in order to be able to calculate the smoothing spline (see Sect.  for details). A nominal uncertainty of 0.3 ppm is assigned to the Mauna Loa data, which is for conservative reasons slightly higher than the generally stated measurement uncertainty of 0.2 ppm (https://www.esrl.noaa.gov/gmd/ccgg/about/co2_measurements.html). Uncertainties for the ice and firn data are taken either as reported or set to 0.5 ppm if the reported standard deviation is missing or less than 0.5 ppm. For Law Dome data published in , we take the reported uniform uncertainty of 1.1 ppm. Note that the adjusted ice core offset between WDC and the other ice cores is not considered in our uncertainty calculation as this represents a systematic error.

The data selection as described above then leads to $n=\mathrm{2152}$ data points including 20 ages with duplicate entries. These duplicates are averaged (reducing $n$ to $\mathrm{2132}$) and the assigned uncertainties based on this averaging.

To account for the variable temporal resolution of the data points (Fig. c) whilst preserving as much of the abrupt changes in CO${}_{\mathrm{2}}$ during Termination I as possible e.g., the spline is divided into 12 intervals with different nominal cutoff periods that vary between 4 years (for instrumental times) and 3000 years (for the Holocene). A low $P_{\mathrm{c}}$ of 600 years was assigned to the high-resolution interval of Termination I (11–18.5 kyr BP). For the glacial interval between 18.5 and 110 kyr BP, $P_{\mathrm{c}}$ of 2000 years was chosen. For the warm interglacial between 110 and 128 kyr BP, we assign a cutoff period of 3000 years similar to the Holocene. Across Termination II (128–135 kyr BP), we use a 1000-year cutoff period, whilst for the penultimate glacial maximum a cutoff period of 2000 years was used. A summary of all details on the calculated spline is found in Table .

The total $\mathrm{1}\mathit{\sigma }$ uncertainty of the spline is $<$ 2 ppm on average (Fig. b). Across some short time windows, it rises up to 6 ppm, and around 108 kyr BP, it reaches a maximum of 11 ppm. The three different error sources contribute equally to the total uncertainty; however, time windows with large uncertainties are often dominated by one error source.

The CO${}_{\mathrm{2}}$ record of the last 2 kyr to be used within CMIP6 is nearly indistinguishable from our spline across the instrumental period (Fig. a); however, CO${}_{\mathrm{2}}$ concentrations during the pre-anthropogenic interval of the last 2 kyr are partly larger by a few parts per million than our spline (Fig. b). This difference is readily explained by the underlying data and the different filtering. We use a combination of Law Dome and WDC data between 200 and 1210 BP, whilst only Law Dome data are considered for CMIP6.

The CO${}_{\mathrm{2}}$ values chosen as boundary conditions for several time slice experiments within PMIP4 can be compared with snapshots from our splines (Table ). However, one needs to be aware that some short-term fluctuations in our spline might offset the values from long-term averages and lead to differences between our final splines and the PMIP4 forcing data. For the mid-Holocene (6 kyr experiment), both our spline and data used in PMIP4 are based on the same EDC data and processed with the identical spline routines and cutoff frequencies, leading to identical values. Values differ by a few parts per million for the experiments 1850 CE, 21 kyr, and 127 kyr.

Since spline smoothing is a low-pass filter, abrupt changes in CO${}_{\mathrm{2}}$ are always smaller in the spline than in the original data sets. Therefore, if one wants to investigate the impact of abrupt increases in CO${}_{\mathrm{2}}$ concentration on the climate system that have been identified during three intervals (around 11.6, 14.7 or 16.2 kyr BP) across Termination I , an alternative continuous CO${}_{\mathrm{2}}$ record needs to be constructed. One approach might be to reduce the cutoff period so that the spline would include these pronounced jumps. For example, one might want to capture the rise in CO${}_{\mathrm{2}}$ of 12 and 13 ppm across a single century at 16.2 and 11.6 kyr BP, respectively, as identified in the WDC record . For the abrupt rise in CO${}_{\mathrm{2}}$ around 14.7 kyr BP, even an increase of 15 ppm in 200 years, slightly larger than the 12 ppm of the WDC record, has been suggested to represent atmospheric changes in CO${}_{\mathrm{2}}$ potentially caused by permafrost thawing during the northern hemispheric warming into the Bølling–Allerød interstadial . Transient simulations investigating these abrupt jumps in CO${}_{\mathrm{2}}$ concentration should use a CO${}_{\mathrm{2}}$ times series that contains greater details than our low-frequency spline.

Data used to construct (or compare to) the Southern Hemisphere CH${}_{\mathrm{4}}$ spline.

Time (in BP)	Time (in CE)	Source	Age scale	Spline${}^{\mathrm{5}}$	Citation
$-$66 to $-$34	2016 to 1984	NOAA network${}^{\mathrm{1}}$ (annual)	–	no
$-$66 to $-$33	2016 to 1983	South Pole${}^{\mathrm{2}}$ (monthly)	–	yes
$-$66 to $-$33	2016 to 1983	Barrow${}^{\mathrm{3}}$ (monthly)	–	no
$-$32 to 168	1982 to 1782	Law Dome${}^{\mathrm{4}}$	as in references	yes
169 to 4669	1781 to before CE	WDC discrete, OSU${}^{\mathrm{6}}$	WD2014	yes
4689 to 9798	–	WDC discrete, PSU${}^{\mathrm{7}}$	WD2014	yes
9821 to 67 233	–	WDC continuous	WD2014	yes
192 to 100 469	–	GRIP	GICC05ext	no
67 401 to 15 6211	–	EDC	AICC2012	yes

Notes: ${}^{\mathrm{1}}$ Global annual mean of the NOAA network. Data taken from ftp://aftp.cmdl.noaa.gov/products/trends/ch4/ch4_annmean_gl.txt. ${}^{\mathrm{2}}$ Data taken from ftp://aftp.cmdl.noaa.gov/data/trace_gases/ch4/flask/surface/ch4_spo_surface-flask_1_ccgg_month.txt. ${}^{\mathrm{3}}$ Data taken from ftp://aftp.cmdl.noaa.gov/data/trace_gases/ch4/flask/surface/ch4_brw_surface-flask_1_ccgg_month.txt. ${}^{\mathrm{4}}$ Law Dome data are taken from various sources; see references for details. They exist for 2005 CE to 14 CE (or $-$55 BP to 1936 BP), but only the part bridging the gap between instrumental data and WDC is taken for the calculation of the spline (1982 CE to 1782 CE or $-$32 BP to 168 BP). ${}^{\mathrm{5}}$ Indicates if the data are used to construct the spline. ${}^{\mathrm{6}}$ Measured in laboratory at Oregon State University (OSU). ${}^{\mathrm{7}}$ Measured in laboratory at Penn State University (PSU); measured data shifted by $+$9.9 ppb to account for unexplained OSU and PSU laboratory offset .

Atmospheric CH${}_{\mathrm{4}}$

Our data compilation of CH${}_{\mathrm{4}}$ data and the consistently calculated CH${}_{\mathrm{4}}$ spline is restricted to the Southern Hemisphere (SH). Northern hemispheric (NH) data are shown for comparison but are not included in the spline, since for such efforts chronologies of ice cores from both hemispheres have to match perfectly during abrupt climate changes of the D/O events. However, as has been shown , there remains some mismatch in the timing of the NH and the SH CH${}_{\mathrm{4}}$ records in the most recent chronology AICC2012. NH CH${}_{\mathrm{4}}$, and consequently global CH${}_{\mathrm{4}}$ concentrations, should, according to the estimates of the interpolar difference, be larger than our SH CH${}_{\mathrm{4}}$ values. Therefore, our SH CH${}_{\mathrm{4}}$ spline represents the lower bound of global CH${}_{\mathrm{4}}$ concentration. found that NH CH${}_{\mathrm{4}}$ was up to $+$4 % ($+$14 ppb) and up to $+$10 % ($+$60 ppb) larger than the SH CH${}_{\mathrm{4}}$ during glacial times and the Holocene, respectively. However, new and as yet unpublished results point to in situ CH${}_{\mathrm{4}}$ production in Greenland ice cores during times of high dust fluxes, calling for a revision of the interpolar difference in CH${}_{\mathrm{4}}$ during glacial times. For this reason, we refrain from calculating an NH or global CH${}_{\mathrm{4}}$ spline. As CH${}_{\mathrm{4}}$ is only of secondary importance for the total greenhouse gas radiative forcing, this systematic error is of little relevance for climate simulation studies. Studies are under way to improve our knowledge of the NH CH${}_{\mathrm{4}}$ value for glacial times, too. Using an approximation of the radiative forcing $$\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}\sim \mathrm{1.4}\cdot \mathrm{0.036}\cdot (\sqrt{{\mathrm{CH}}_{\mathrm{4}}/\mathrm{ppb}}-\sqrt{\mathrm{742}})\mathrm{W}{\mathrm{m}}^{-\mathrm{2}},$$ which neglects the interacting effects of CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O but which considers the approximate increase in $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$ by 40 % through indirect effects of CH${}_{\mathrm{4}}$ on stratospheric H${}_{\mathrm{2}}$O and tropospheric O${}_{\mathrm{3}}$ , we estimate that our restriction of CH${}_{\mathrm{4}}$ to the SH only would underestimate the radiative forcing of CH${}_{\mathrm{4}}$ by less than 0.05 W m${}^{-\mathrm{2}}$.

Our data compilation starts with the beginning of the year 2016 CE ($-\mathrm{66.0}$ BP) and stops around 156 kyr BP to cover the same time interval as for CO${}_{\mathrm{2}}$ (Fig. a). The 11-point running mean temporal resolution between neighbouring data points, $\mathrm{\Delta }t$, is less than 100 years for most of the last 67 kyr, increasing to $\sim$ 700 years between 67 and 156 kyr BP (Fig. c). Our strategy here is to select one data set for each point in time and use overlapping intervals only for confirmation of data consistency. The following data sets are considered here.

CH${}_{\mathrm{4}}$ spline covering all data: 2016 CE–156 211 BP. Details on plotted data are explained in the text. The maximum ice core data uncertainty ($\pm \mathrm{2}\mathit{\sigma }$) is given in the lower left corner. In (a) the right axis contains the resulting radiative forcing approximated with $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}\sim \mathrm{1.4}\cdot \mathrm{0.036}\cdot (\sqrt{{\mathrm{CH}}_{\mathrm{4}}/\mathrm{ppb}}-\sqrt{\mathrm{742}})$ W m${}^{-\mathrm{2}}$ based on , but neglecting interacting effects of CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O and considering indirect effects of CH${}_{\mathrm{4}}$ on stratospheric H${}_{\mathrm{2}}$O and tropospheric O${}_{\mathrm{3}}$ . The latitudinal origin of data is indicated by NH and SH, indicating Northern and Southern Hemisphere, respectively. (b) Total uncertainty of the spline based on three individual error sources; see text for details. (c) Temporal resolution ($\mathrm{\Delta }t$) of the CH${}_{\mathrm{4}}$ data points underlying the spline on a log scale. Additionally, the prescribed time-dependent cutoff period $P_{\mathrm{c}}$ is plotted, including its variation by $\pm$50 %, which has been used to determine ${\mathit{\sigma }}_{\mathrm{1}}$.

[Figure omitted. See PDF]

Details of the southern hemispheric CH${}_{\mathrm{4}}$ spline. Light and dark grey bands around the spline represents $\pm \mathrm{1}\mathit{\sigma }$ and $\pm \mathrm{2}\mathit{\sigma }$, respectively. (a) Instrumental times (1950–2016 CE); (b) 0–2000 BP; (c) Termination I; (d) 0–40 kyr BP without the Law Dome data showing the anthropogenic rise; (e) 40–90 kyr BP; (f) 90–160 kyr BP. Hemispheric origin of data is indicated by NH (north) and SH (south). WDC PSU data are adjusted by $+\mathrm{9.9}$ ppb. GRIP+: Greenland composite; GICC05+: GICC05 model extended. See text for details. Dashed line labelled CMIP6 in panels (a) and (b) is the compiled CH${}_{\mathrm{4}}$ record to be used in CMIP6 experiments for the last 2 kyr .

[Figure omitted. See PDF]

Statistics of the CH${}_{\mathrm{4}}$ spline. Interval: $i_{{\mathrm{CH}}_{\mathrm{4}}}$; $s$: scaling factor to fulfil the constraints given by the prescribed cutoff period $P_{\mathrm{c}}$; $\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$: average realised cutoff period; $\mathrm{\Delta }t$: mean data spacing; $v$: mean $\mathrm{1}\mathit{\sigma }$ error – exact scaling factor to fulfil the constraints given by the prescribed cutoff period $P_{\mathrm{c}}$; $\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$: average realised cutoff period. $\mathrm{\Delta }t$: mean data spacing; $v$: mean $\mathrm{1}\mathit{\sigma }$ error; exact time framing is given by the age of the first ($t_{\mathrm{start}}$) and last ($t_{\mathrm{stop}}$) data point of the interval (in years BP); $N$: number of data points within interval. In the last column the underlying data source is briefly mentioned; see Table  for details and citations.

$i_{{\mathrm{CH}}_{\mathrm{4}}}$	$s$	$P_{\mathrm{c}}$	$\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$	$\mathrm{\Delta }t$	$v$	$t_{\mathrm{start}}$	$t_{\mathrm{stop}}$	$N$	Data source
	–	yr	yr	yr	ppb	yr BP	yr BP	–
1	1.00	4.0	4.0	0.1	2.0	$-$66.0	$-$33.1	396	South Pole
2	0.68	10.0	9.5	1.7	4.0	$-$31.8	165.1	114	Law Dome
3	13.12	50.0	49.7	8.2	2.4	169.0	2591.0	296	WDC discrete
4	33.63	100.0	98.8	20.0	2.4	2602.0	9798.0	361	WDC discrete
5	3.09	20.0	20.0	2.0	3.3	9821.0	67 233.0	28 707	WDC continuous
6	56.26	500.0	497.0	257.1	10.0	67 401.0	127 831.0	236	EDC
7	171.96	1000.0	982.0	440.4	10.0	128 026.0	156 211.0	65	EDC

• From the NOAA network, the annual global mean concentration of CH${}_{\mathrm{4}}$ from 2016 CE to 1984 CE is available (www.esrl.noaa.gov/gmd/ccgg/trends/). These global mean concentrations lie in between the seasonally resolved CH${}_{\mathrm{4}}$ concentration measured at Barrow, Alaska (NH), and at the South Pole (SH), both reaching back in time until 1983 CE (Fig. a). The interpolar difference between the NH (Barrow) and the SH (South Pole) was $+$161 ppb or $+$9 % at the beginning of 2016 CE. In absolute CH${}_{\mathrm{4}}$ concentration, this most recent interpolar difference is $\sim$ 100 ppb larger than the interpolar difference in the Holocene, while the relative interpolar difference during both time intervals is comparable . An estimation of the radiative forcing of this interpolar difference reveals that, for the time covered by the NOAA network, the $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$ for the NH was $<\mathrm{0.1}$ W m${}^{-\mathrm{2}}$ larger than for the SH. This estimate of the radiative forcing of the interpolar difference is obtained from Eq. (), based on Barrow and the South Pole CH${}_{\mathrm{4}}$ data. For our SH compilation we used the South Pole data.

• Ice core and firn air data from Law Dome and Cape Grim (SH) exist from 2005 CE back to 14 CE ($=$ 1936 BP) with an overlap of more than 2 decades with the instrumental measurements (Fig. a, b). The CH${}_{\mathrm{4}}$ data from Law Dome and Cape Grim used in our compilation span the period from 1982 CE to 1782 CE ($=$ 168 BP), bridging the gap between instrumental data and CH${}_{\mathrm{4}}$ from WDC. Where the Law Dome data overlap with the data from either the South Pole or WDC, no apparent systematic offsets between the different data sets have been identified. WDC and Law Dome data differ slightly across short intervals between 1000 and 2000 BP (Fig. b). However, since WDC is the more highly resolved record, it alone is included in the spline; no data adjustment is necessary here.

• The discrete CH${}_{\mathrm{4}}$ data from WDC (SH) span the interval 169 BP to 67 kyr BP . Starting with the year 9821 BP, continuous CH${}_{\mathrm{4}}$ data from WDC with higher temporal resolution are now available and are used to support our spline . These continuous CH${}_{\mathrm{4}}$ data have already been post-processed, including the support from some discrete WDC data points to improve the data set, whenever larger gaps in the continuous record appeared. The data product of the continuous CH${}_{\mathrm{4}}$ WDC data obtained at NOAA is splined to a constant temporal resolution of 2 years. For the missing part of the Holocene not contained in the continuous WDC data, discrete WDC CH${}_{\mathrm{4}}$ data are used. They have been measured in two different laboratories: at Oregon State University (OSU; 169–4669 BP) and at Pennsylvania State University (PSU; 4689–9798 BP). An unexplained inter-laboratory offset between the discrete CH${}_{\mathrm{4}}$ WDC data from OSU and PSU has been identified. To account for this offset the PSU CH${}_{\mathrm{4}}$ data have been adjusted by $+\mathrm{9.9}$ ppb Supplementary Information to. To date WDC CH${}_{\mathrm{4}}$ are the temporally most highly resolved data of the last glacial, and therefore they are our reference record (Fig. b–e). The data not only contain the well-known abrupt CH${}_{\mathrm{4}}$ changes at the onset and end of the millennial-scale D/O events in high resolution and accuracy but also centennial-scale features that are understood to be of climatic origin e.g..

• We extend our SH CH${}_{\mathrm{4}}$ data compilation beyond WDC with data from EDC, spanning the period from 67 to 156 kyr BP (Fig. e–f). These EDC data actually extend back to 800 kyr BP, but since our focus here is on the time since the penultimate glacial maximum (i.e. the last 156 kyr), the CH${}_{\mathrm{4}}$ record older than 156 kyr is not considered here. CH${}_{\mathrm{4}}$ data from EDML might be in part more highly resolved than in EDC because of a higher annual layer thickness between 67 kyr BP (the end of WDC) and 80 kyr BP . However, a well-documented EDML CH${}_{\mathrm{4}}$ record is not available to date, and therefore none of the published EDML CH${}_{\mathrm{4}}$ data for this interval are considered here.

• The NH Greenland composite of CH${}_{\mathrm{4}}$ is only plotted for comparison to the SH data (Fig. b–f).

Data used to construct the N${}_{\mathrm{2}}$O spline.

Time (in BP)	Time (in CE)	Source	Age scale	Citation
$-$66 to $-$49	2016 to 1999	NOAA network (monthly)	–	nitrous oxide data from the NOAA/ESRL halocarbons in situ program${}^{\mathrm{1}}$
$-$50 to $-$38	2000 to 1988	NOAA network (monthly)	–	RITS nitrous oxide data from the NOAA/ESRL halocarbons program${}^{\mathrm{2}}$
$-$33 to 1937	1983 to 13	Law Dome${}^{\mathrm{3}}$	as in references
1975 to 11 502	–	EDC${}^{\mathrm{4}}$	AICC2012
29 065 to 134 519	–	EDC	AICC2012
9858 to 15 843	–	Taylor Glacier	WD2014${}^{\mathrm{5}}$
15 000 to 118 602	–	NGRIP${}^{\mathrm{6}}$	AICC2012
15 000 to 134 418	–	Talos Dome${}^{\mathrm{7}}$	AICC2012

Notes: ${}^{\mathrm{1}}$ ftp://ftp.cmdl.noaa.gov/hats/n2o/insituGCs/CATS/global/insitu_global_N2O.txt. ${}^{\mathrm{2}}$ ftp://ftp.cmdl.noaa.gov/hats/n2o/insituGCs/RITS/global/RITS_global_N2O.txt. ${}^{\mathrm{3}}$ Law Dome data are taken from various sources; see references for details. They exist from 2004 CE to 13 CE (or $-$54 BP to 1937 BP), but only those older than the instrumental record (1986 CE and older) are taken here. ${}^{\mathrm{4}}$ Data exist from 334 BP (or 1616 CE) until 11 502 BP, but only data not yet covered by the Law Dome records (13 CE or 1937 BP and older) are considered here. ${}^{\mathrm{5}}$ WD2014 age model for Taylor Glacier, published in . ${}^{\mathrm{6}}$ Data exist for 11 068 BP–119 555 BP, but only those older than 15 kyr BP are considered here. Five data points in the oldest part considerably disagree with SH records and are therefore rejected. ${}^{\mathrm{7}}$ Data exist for 217 BP–134 418 BP, but only those older than 15 kyr BP are considered here.

The assigned data uncertainty ($\mathrm{1}\mathit{\sigma }$ error) is 2.0, 4.0, 2.4, 3.3, and 10 ppb for instrumental data, Law Dome, discrete WDC, continuous WDC, and EDC, respectively . Using the approximation of $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$, given the above $\mathrm{1}\mathit{\sigma }$ errors the uncertainty in the radiative forcing is $<$0.01 W m${}^{-\mathrm{2}}$.

Compiled data contain 30 214 data points, among which duplicate entries exist for 39 ages. These duplicates are averaged giving $n=\mathrm{30}\mathrm{175}$.

The whole data set is divided into seven intervals with different assigned cutoff periods. $P_{\mathrm{c}}$ ranges from 4 years (for the interval covered by instrumental data) to 20 years (for the interval covered by the continuous WDC record). Due to lower data coverage during the Holocene and further back in time, $P_{\mathrm{c}}$ is increased to 100 years (0.2–9.8 kyr BP), 500 years (60–128 kyr BP), and 1000 years (128–156 kyr BP) (Fig. c). More details are shown in Table .

The total $\mathrm{1}\mathit{\sigma }$ uncertainty of our final CH${}_{\mathrm{4}}$ spline is around 3–10 ppb in the Holocene, $\sim$ 2 ppb in the time window supported by the continuous WDC CH${}_{\mathrm{4}}$ data (9.8–67 kyr BP), and around 10 ppb in earlier parts. During some short time windows, ${\mathit{\sigma }}_{{\mathrm{CH}}_{\mathrm{4}}}$ reaches a maximum of 20 ppb (Fig. b). The uncertainty is dominated by the Monte Carlo error before 9.8 kyr BP and by the error in the cutoff period in the Holocene and during those short events in which ${\mathit{\sigma }}_{{\mathrm{CH}}_{\mathrm{4}}}$ reached its local maxima. An abrupt jump in ${\mathit{\sigma }}_{{\mathrm{CH}}_{\mathrm{4}}}$ appears at 67 kyr BP (transition from continuous WDC to discrete EDC data), when individual data point uncertainty rose from 3.3 to 10 ppb (changing ${\mathit{\sigma }}_{\mathrm{3}}$) at the same time as $\mathrm{\Delta }t$, and therefore $P_{\mathrm{c}}$, increased by 2 orders of magnitude (changing ${\mathit{\sigma }}_{\mathrm{1}}$).

N${}_{\mathrm{2}}$O spline covering all data: 2016 CE–134 519 BP. Details on plotted data are explained in the text. The maximum ice core data uncertainty ($\pm \mathrm{2}\mathit{\sigma }$) is sketched in the lower left corner. In (a) the right axis contains the resulting radiative forcing approximated with $\mathrm{\Delta }{R}_{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}\sim \mathrm{0.12}\cdot (\sqrt{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/\mathrm{ppb}}-\sqrt{\mathrm{272}})$ W m${}^{-\mathrm{2}}$ after , neglecting interacting effects of CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O. Filled symbols: data taken for spline; open symbols: data not taken for spline. (b) Total uncertainty of the spline based on three individual error sources; see text for details. (c) Temporal resolution ($\mathrm{\Delta }t$) of the N${}_{\mathrm{2}}$O data points underlying the spline on a log scale. Additionally, the prescribed time-dependent cutoff period $P_{\mathrm{c}}$ is plotted, including its variation by $\pm$50 %, which has been used to determine ${\mathit{\sigma }}_{\mathrm{1}}$.

[Figure omitted. See PDF]

Details of the N${}_{\mathrm{2}}$O spline. Light and dark grey bands around the spline represent $\pm \mathrm{1}\mathit{\sigma }$ and $\pm \mathrm{2}\mathit{\sigma }$, respectively. (a) Instrumental times (1950–2016 CE); (b) 0–2000 BP; (c) Termination I; (d) 0–40 kyr BP without the Law Dome data showing the anthropogenic rise; (e) 40–90 kyr BP; (f) 90–140 kyr BP. Filled symbols: data taken for spline; open symbols: data not taken for spline. See text for further details. Dashed line labelled CMIP6 in panels (a) and (b) is the compiled N${}_{\mathrm{2}}$O record to be used in CMIP6 experiments for the last 2 kyr .

[Figure omitted. See PDF]

Statistics of N${}_{\mathrm{2}}$O spline. Interval: $i_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$; $s$: scaling factor to fulfil the constraints given by the prescribed cutoff period $P_{\mathrm{c}}$; $\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$: average realised cutoff period; $\mathrm{\Delta }t$: mean data spacing; $v$: mean $\mathrm{1}\mathit{\sigma }$ error – exact time framing is given by the age of the first ($t_{\mathrm{start}}$) and last ($t_{\mathrm{stop}}$) data point of the interval (in years BP); $N$: number of data points within interval. In the last column the underlying data source is briefly mentioned; see Table  for details and citations.

$i_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$	$s$	$P_{\mathrm{c}}$	$\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$	$\mathrm{\Delta }t$	$v$	$t_{\mathrm{start}}$	$t_{\mathrm{stop}}$	$N$	Data source
	–	yr	yr	yr	ppb	yr BP	yr BP	–
1	1.00	4.0	3.4	0.1	0.8	$-$66.0	$-$38.3	334	NOAA network
2	3.01	50.0	48.7	2.5	7.2	$-$33.7	95.0	53	Law Dome
3	19.08	200.0	190.2	16.8	7.0	104.0	389.6	18	Law Dome
4	321.85	1000.0	952.4	85.1	4.6	400.3	9425.6	107	Law Dome, EDC
5	77.33	500.0	469.0	58.2	5.8	9517.2	15 974.7	112	EDC, Taylor Glacier, NGRIP, Talos Dome
6	1443.10	2000.0	1915.3	60.4	4.9	16 003.0	116 900.0	1672	EDC, NGRIP, Talos Dome
7	4236.11	5000.0	4792.5	370.0	4.2	117 130.0	134 519.0	48	EDC, NGRIP, Talos Dome

The SH CH${}_{\mathrm{4}}$ record to be used within CMIP6 largely agrees with our SH CH${}_{\mathrm{4}}$ spline (Fig. a, b). However, during instrumental times the CMIP6 SH CH${}_{\mathrm{4}}$ record is consistently larger than our SH CH${}_{\mathrm{4}}$ spline by about 10–15 ppb, probably caused by the inclusion of different stations in the calculation of the SH CH${}_{\mathrm{4}}$ record within CMIP6, while we rely on South Pole data. Prior to the instrumental CH${}_{\mathrm{4}}$ data around 1980 CE, the maximum difference between both approaches is 30 ppb. This difference might be caused by the statistical routines within CMIP6 to account for missing stations. Further back in time (around 1150 BP, 1300 BP, and 1900 BP), higher-frequency variation contained in the WDC CH${}_{\mathrm{4}}$ record (used here but ignored within CMIP6) leads to some CH${}_{\mathrm{4}}$ variations within our SH CH${}_{\mathrm{4}}$ spline on the order 10–25 ppb that are not captured by the CMIP6 SH CH${}_{\mathrm{4}}$ record.

A comparison of our final spline with the GHG values chosen for the PMIP4 time slice experiments is not straightforward, since we only compile SH CH${}_{\mathrm{4}}$ data, while the PMIP4 experiments use global values. Taking the two records at face value, one finds that our SH CH${}_{\mathrm{4}}$ is 13, 44, and 25 ppb smaller than the global mean value used in PMIP4 for 1850 CE, 6 kyr, and 127 kyr, respectively. In particular, the large SH-global difference of 44 ppb around 6 kyr seems to be rather large but is readily explained by the centennial variability contained in the WDC CH${}_{\mathrm{4}}$, which leads to a local minimum in SH CH${}_{\mathrm{4}}$ around 6 ka. Similarly, our SH CH${}_{\mathrm{4}}$ spline is 7 ppb higher than the global CH${}_{\mathrm{4}}$ value chosen within PMIP4 for the 21 kyr experiment. This difference can again be explained by the centennial-scale variability contained in the WDC CH${}_{\mathrm{4}}$ record, which shows a local maximum at 21 kyr BP. A hundred years later, our SH CH${}_{\mathrm{4}}$ spline has a local minimum which is 11 ppb smaller than the global CH${}_{\mathrm{4}}$ values taken for PMIP4 (Table ).

Atmospheric N${}_{\mathrm{2}}$O

For the data compilation of the third GHG, N${}_{\mathrm{2}}$O, one has to be aware that during times of high dust input, in situ production of N${}_{\mathrm{2}}$O occurs, leading to artefacts in the paleo record . Furthermore, the precise synchronisation of Northern and Southern Hemisphere records, as already explained for CH${}_{\mathrm{4}}$, is crucial to accurately obtain the changes in N${}_{\mathrm{2}}$O during millennial-scale D/O events.

The compiled record starts at the beginning of the year 2016 CE ($-\mathrm{66.0}$ BP) but extends back in time only until $\sim$134.5 kyr BP (Fig. a) because the ice cores on which the N${}_{\mathrm{2}}$O compilation is based in the older parts, Talos Dome, EDC, and NGRIP, have either no data points between 134.5 and 156 kyr BP or unreliable N${}_{\mathrm{2}}$O data containing artefacts across the penultimate glacial maximum . The latter is also the case for EDML, whose data have not been taken to support the spline because despite the agreement of the N${}_{\mathrm{2}}$O of EDML and EDC, the data from EDML have a lower temporal resolution than those of EDC .

The data sets contributing to the N${}_{\mathrm{2}}$O stack are listed below.

• There are two contributions of N${}_{\mathrm{2}}$O data based on instrumental measurements to the NOAA network or ESRL halocarbon program: (a) in situ N${}_{\mathrm{2}}$O data are available from 2016 CE back until 1999 CE, and (b) the RITS N${}_{\mathrm{2}}$O data from 2000 CE go back until 1988 CE. Both represent global mean monthly values (Fig. a). Note that due to the long atmospheric lifetime of N${}_{\mathrm{2}}$O, any interpolar difference can be safely neglected.

• Law Dome and Cape Grim N${}_{\mathrm{2}}$O data exist from 2004 CE back until 13 CE (1937 BP) and correlate well with the instrumental data in overlapping intervals (Fig. a,b). Here, the Law Dome data contribute to the spline only for those years not covered by the instrumental record, i.e. 1983 CE and earlier.

• In the Holocene, N${}_{\mathrm{2}}$O was measured at EDC from 334 BP until 11.5 kyr BP. For the last two millennia, the EDC N${}_{\mathrm{2}}$O data points are sparser than the Law Dome data; therefore, the EDC N${}_{\mathrm{2}}$O data are only considered for times older than what is covered in the Law Dome N${}_{\mathrm{2}}$O record, i.e. before 1975 BP (Fig. b, d).

• The most highly resolved N${}_{\mathrm{2}}$O record across large parts of Termination I is provided by the horizontal ice core from Taylor Glacier which has been linked to the chronology of the WDC ice core (WD2014) via CH${}_{\mathrm{4}}$ . The Taylor Glacier N${}_{\mathrm{2}}$O record used in our spline covers the interval 9.6 to 15.8 kyr BP (Fig. c).

• The last glacial interval is well resolved by N${}_{\mathrm{2}}$O data from the NGRIP record . While the NGRIP N${}_{\mathrm{2}}$O data cover the times between 11 kyr BP and 119.6 kyr BP, we only take those data older than 15 kyr BP due to the more highly resolved Taylor Glacier N${}_{\mathrm{2}}$O data during Termination I (Fig. c–f). Five data points near the bedrock in the bottom part of the NGRIP records apparently have higher N${}_{\mathrm{2}}$O values than found in ice cores from the Southern Hemisphere. These data points are rejected here, leading to the oldest NGRIP N${}_{\mathrm{2}}$O data point at 118.6 kyr BP. We are aware that due to the imperfect north–south synchronisation of gas records in AICC2012 (see Sect.  for details), the usage of N${}_{\mathrm{2}}$O data from NGRIP might introduce erroneous phasing between our global N${}_{\mathrm{2}}$O record and the purely SH-based CH${}_{\mathrm{4}}$ spline, particularly during abrupt change connected to D/O events. However, N${}_{\mathrm{2}}$O data coverage in the SH is very sparse and a spline only based on SH data would be even less reliable. This potential synchronisation problem is also addressed by large cutoff periods of the spline of 2000 to 5000 years beyond 16 kyr BP.

• Additional N${}_{\mathrm{2}}$O data going back to 134.4 kyr BP are obtained from the Talos Dome ice core and from further data of the EDC ice core (compilation found in ; data source between 29.0–134.5 kyr BP: ). Since – besides EDML, which correlates well with EDC – these are the only N${}_{\mathrm{2}}$O records with reliable data going back to the penultimate glacial maximum, we consider all data points from the Talos Dome and EDC ice cores here before 15 kyr BP. The data points of Talos Dome and EDC in general agree with the NGRIP data over the last glacial cycle, but NGRIP diverges from the SH records towards higher (probably biased) values in the warm previous interglacial around 115 kyr BP. As already explained above, these five NGRIP data points are rejected. However, across Termination I, Talos Dome N${}_{\mathrm{2}}$O data seems to be systematically lower than NGRIP N${}_{\mathrm{2}}$O data, with Taylor Glacier data in between both (Fig. c). We therefore believe that a mixture of all three records (N${}_{\mathrm{2}}$O from NGRIP, EDC, and Talos Dome) most likely represents a reasonable mean global N${}_{\mathrm{2}}$O value (Fig. c–f). The relatively large difference in N${}_{\mathrm{2}}$O from different ice cores during the last glacial times indicates that the uncertainty (accuracy) in N${}_{\mathrm{2}}$O is probably higher than the reported measurement errors (precision) of up to 7 ppb.

Calculated radiative forcing of CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, N${}_{\mathrm{2}}$O, and their sum ($\mathrm{\Delta }{R}_{\left[\mathrm{GHG}\right]}$), including $\mathrm{1}\mathit{\sigma }$-uncertainty bands. In addition to the uncertainty of the spline, further uncertainties need to be considered: 10 % relative uncertainty contained in the simplified Eqs. ()–(); $\mathrm{0.2}$ W m${}^{-\mathrm{2}}$ uncertainty in $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$ and $\mathrm{\Delta }{R}_{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}$ due to the omission of interaction effects; 5 % uncertainty in efficiency of CH${}_{\mathrm{4}}$ for more details, see. The total uncertainty of the $\mathrm{\Delta }{R}_{\left[\mathrm{GHG}\right]}$ is calculated as the square root of the sum of squares of the uncertainties of the individual GHGs. The calculations are based on Eqs. ()–() given in the text following. Sub-figures focus on specific time windows: (a) anthropogenic rise since 1750 CE; (b) Termination I; (c) 20–90 kyr BP, including the abrupt changes during D/O event; (d) full record from 2016 CE to 156 kyr BP – here N${}_{\mathrm{2}}$O was kept constant beyond 134 kyr BP.

[Figure omitted. See PDF]

The generally assigned $\mathrm{1}\mathit{\sigma }$ uncertainty of each data point is 7 ppb for the Law Dome ice core . The uncertainty of individual data points in other ice cores was in general less than 7 ppb . For the instrumental measurements, we take the reported uncertainties of around 1 ppb. For 58 times, more than one data point for the same age exists. These duplicates are averaged reducing the number of N${}_{\mathrm{2}}$O data to $n=\mathrm{2344}$. Using an estimate of the radiative forcing of N${}_{\mathrm{2}}$O, which neglects the interacting effects of CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O, $$\mathrm{\Delta }{R}_{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}\sim \mathrm{0.12}\cdot (\sqrt{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/\mathrm{ppb}}-\sqrt{\mathrm{272}})\mathrm{W}{\mathrm{m}}^{-\mathrm{2}},$$ , we estimate that the $\mathrm{1}\mathit{\sigma }$ error in N${}_{\mathrm{2}}$O is related to an uncertainty in the radiative forcing of about 0.04 W m${}^{-\mathrm{2}}$, slightly larger than the uncertainty in $\mathrm{\Delta }R$ related to the CH${}_{\mathrm{4}}$ data. Comparing the different values of N${}_{\mathrm{2}}$O in Talos Dome and NGRIP for same intervals reveals differences on the order of about 10 ppb (e.g. Fig. c–f), suggesting that the ice-core-specific values of N${}_{\mathrm{2}}$O contain an intrinsic uncertainty which is comparable to the measurement error.

The mean temporal resolution (11-point running mean) of the underlying N${}_{\mathrm{2}}$O data is around 50 years across large parts of the last glacial cycle (15–60 kyr BP), with slightly lower resolution of 100 years in the Holocene and between 60 and 115 kyr BP. In MIS5.5, Termination II, and the penultimate glacial maximum, the mean temporal resolution rises to $\sim$ 500 years (Fig. c). Based on this distribution of $\mathrm{\Delta }t$ the prescribed cutoff periods for the spline vary for seven different intervals between 4 (for the instrumental period) and 5000 years (for data older than 117 kyr BP). For the majority of the data (400 yr BP to 117 kyr BP), a $P_{\mathrm{c}}$ between 500 and 2000 years is prescribed. More details on the spline are found in Table .

The total $\mathrm{1}\mathit{\sigma }$ uncertainty of the spline varies between 1 and 6 ppb (Fig. b). This uncertainty is mainly based on ${\mathit{\sigma }}_{\mathrm{3}}$ (the error related to the Monte Carlo statistics) for periods younger than the LGM. For intervals older than the LGM, the main uncertainty is ${\mathit{\sigma }}_{\mathrm{1}}$, the error related to the cutoff period.

If compared with the N${}_{\mathrm{2}}$O compilation used within CMIP6 both approaches largely agree for instrumental times (Fig. a). Further back in time during the last 2 kyr, both approaches rely on the same data: the published Law Dome/Cape Grim N${}_{\mathrm{2}}$O data . Interestingly, both time series differ by up to 6 ppb between 0.7 and 2.0 kyr BP (Fig. b). This difference is in the range of the ice core data uncertainty, and therefore still small, but we have no ready explanation. The records used in CMIP6 have a higher N${}_{\mathrm{2}}$O concentration than all data from Law Dome or other SH ice cores (Fig. b), for some unknown reason.

The N${}_{\mathrm{2}}$O data used as starting values in the PMIP4 experiments 1850 CE, 6 kyr, 21 kyr, and 127 kyr agree within 1 or 2 ppb with values based on our calculated spline; only for 21 kyr is the offset with 6 ppb greater (Table ).

Data connected with this paper are available in the scientific database PANGAEA (10.1594/PANGAEA.871273).

In detail, for each of three GHGs the following data are available.

• Final, compiled raw data ($t_j$, $y_j$, and $v_j$, corresponding to time, value, and assumed $\mathrm{1}\mathit{\sigma }$ error), including the data source, as described in the article.

• Preprocessed raw data (averaging of duplicate entries for similar times).

• Calculated splines with time steps of $\mathrm{\Delta }t=\mathrm{1}$ year, including the $\mathrm{1}\mathit{\sigma }$ total uncertainty.

• Corresponding radiative forcing based on the simplified Eqs. ()–(), including the $\mathrm{1}\mathit{\sigma }$ total uncertainty.

When using these data, please consider citing the original publications from which the data underlying this compilation have been taken.

Conclusions

We have compiled available greenhouse gas records and, by calculating a smoothing spline, we were able to provide continuous records over the last glacial cycle, starting from the beginning of the year 2016 CE and going back to 134 kyr BP (for N${}_{\mathrm{2}}$O) and to 156 kyr BP (for CO${}_{\mathrm{2}}$ and CH${}_{\mathrm{4}}$). These records should serve as boundary conditions to calculate the greenhouse gas radiative forcing in transient climate simulations as planned, for example, in the German project PALMOD, or they might be used in the Last Deglaciation experiment within PMIP4 or in other future model intercomparison projects. The resulting radiative forcing of the three GHGs calculated here with the simplified non-interacting Eqs. ()–() has relative uncertainties of $\sim$ 10 % . These equations have also been used in and , while used a different equation for $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$. The latter two studies furthermore include the interacting effects of CH${}_{\mathrm{4}}$ and N${}_{\mathrm{2}}$O (which we ignore here) but neglect the 40 % increase in CH${}_{\mathrm{4}}$ radiative forcing due to indirect effects of CH${}_{\mathrm{4}}$ . Our results are very similar to recent calculations based on a complete and revised set of simplified equations, which also consider interacting effects between the three GHGs , with differences between old and new expressions in $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]},\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$, and $\mathrm{\Delta }{R}_{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}$ of less than 0.01, 0.04, and 0.02 W m${}^{-\mathrm{2}}$, respectively. While the differences in $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ and $\mathrm{\Delta }{R}_{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}$ lie within their uncertainty bands, $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$ is slightly higher, leading to a revised relative uncertainty of $\mathrm{14}$ % . We have refrained from applying the new equations throughout the study, since the amplification in $\mathrm{\Delta }{R}_{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}$ by 40 % through indirect effects of CH${}_{\mathrm{4}}$ is not considered. We prefer to estimate the radiative forcing attributable to one of the three GHGs individually, without interacting effects. Our forcing calculations clearly illustrate the dominant contribution of CO${}_{\mathrm{2}}$, which is responsible for about two-thirds of the total radiative forcing $\mathrm{\Delta }{R}_{\left[\mathrm{GHG}\right]}$ during both the anthropogenic rise (Fig. a) and the reduction during the LGM (Fig. b). High-resolution variability in CH${}_{\mathrm{4}}$ (captured due to smaller cutoff periods during spline calculations than for CO${}_{\mathrm{2}}$) also imposes some fine-scale structure on the overall GHG radiative forcing (Fig. c, d); however, the dominant features are still driven by CO${}_{\mathrm{2}}$.

Details of differences between CO${}_{\mathrm{2}}$ of Law Dome , EDC , and WDC during the last 1500 years, showing how the adjustment of the WDC has been calculated. The grey area marks the pre-anthropogenic time window (before 1750 CE) covered in WDC (200–1210 BP), from which the difference in CO${}_{\mathrm{2}}$ to WDC and Law Dome records has been determined. Horizontal lines mark the mean values for the different ice cores (cyan: Law Dome (all data); magenta: WDC). The mean offset between the WDC and Law Dome of 3.13 ppm is subtracted from the WDC data in (b).

[Figure omitted. See PDF]

Details of differences between CO${}_{\mathrm{2}}$ of EDC and WDC during Termination I, showing how the adjustment of the WDC has been calculated. Grey areas mark the three time windows with relatively stable CO${}_{\mathrm{2}}$ from which the difference in CO${}_{\mathrm{2}}$ to both records has been determined. Horizontal lines mark the mean values for the different ice cores (green: EDC; magenta: WDC). The duration-weighted mean offset between the WDC and EDC of 6.06 ppm is subtracted from the WDC data in (b).

[Figure omitted. See PDF]

PK initiated the work, compiled the data, calculated the spline, and led the writing of the paper. CNA, JS, TFS, and HF contributed specific insights on the data selection and advised on the spline smoothing. All co-authors commented on and improved the initial draft.

The authors declare that they have no conflict of interest.

Acknowledgements

We thank NOAA for the availability of the instrumental GHG data, specifically Ed Dlugokencky and Pieter Tans, NOAA/ESRL (www.esrl.noaa.gov/gmd/ccgg/trends/), and Ralph Keeling, Scripps Institution of Oceanography (scrippsco2.ucsd.edu/). Furthermore, we thank Fortunat Joos for some earlier thoughts on similar GHG data compilations and for providing the Fortran routines to calculate the spline fits; Jenny Roberts for improving the English; Thomas Bauska and two anonymous reviewers for their comments; and Rainer Sieger (PANGAEA) for comments on the uploaded data, which led to the preparation of Table . This is a contribution to the German BMBF project PALMOD. Long-term support of ice core research at the University of Bern by the Swiss National Science Foundation (SNF) is gratefully acknowledged. Edited by: Alexander Kokhanovsky Reviewed by: Thomas K. Bauska and two anonymous referees