A 156 kyr smoothed history of the atmospheric greenhouse gases CO 2 , CH 4 , and N 2 O and their radiative forcing

. Continuous records of the atmospheric greenhouse gases (GHGs) CO 2 , CH 4 , and N 2 O are necessary input data for transient climate simulations, and their associated radiative forcing represents important compo-nents 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 interpo-lating 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 ﬁrn air and atmospheric samples spanning the interval from the penultimate glacial maximum ( ∼ 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 4 and discuss how much larger a Northern Hemisphere or global CH 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, deﬁning 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 ﬁnal splines, are available at https://doi.org/10.1594/PANGAEA.871273.


Introduction
Our knowledge of changes in the atmospheric mixing ratios of the important greenhouse gases (GHGs) CO 2 , CH 4 , and N 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 4 records using gas extraction and measurement systems coupled to continuous-flow analysis systems (Schüpbach et al., 2009;Chappellaz et al., 2013;Rhodes et al., 2013, 364 P. Köhler et al.: 156 kyr smoothed history of CO 2 , CH 4 , and N 2 O -For some of the CO 2 records obtained from different ice cores, there exist significant and as yet unexplained offsets (Ahn et al., 2012;Bereiter et al., 2012;Marcott et al., 2014;Bauska et al., 2015). These offsets need to be addressed in our data compilation.
-Due to the dominance of CH 4 sources in the Northern Hemisphere, the CH 4 concentrations are higher in records from Greenland than from Antarctica (referred to as interpolar difference; e.g. Baumgartner et al., 2012).
-In situ production of N 2 O connected to high mineral dust values leads to unreliable N 2 O concentrations (e.g. Schilt et al., 2010a), 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 4 and N 2 O (and to some extent also in CO 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 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 (Baumgartner et al., 2014) 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 (WAIS Divide Project Members, 2015).
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 2 , CH 4 , and N 2 O via splinesmoothing (Enting, 1987;Bruno and Joos, 1997) with a nominal temporal resolution 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 t represents not the true resolution but only the typical spline average for each year and that the ice core information represents a lowpass 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) (Ivanovic et al., 2016). Note that different GHG data sets have so far been chosen to force transient simulations for the last 21 kyr in Ivanovic et al. (2016), but well-motivated different setups (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 (Otto-Bliesner et al., 2016). 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 (Meinshausen et al., 2017). 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 (MacFarling-Meure et al., 2006;Rubino et al., 2013). We will finally compare our splines with these forcing data sets proposed by Meinshausen et al. (2017) 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 (1σ errors) and are of the order of a few parts per million for CO 2 or a few parts per billion for CH 4 and N 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 (Myhre et al., 1998) with an embedded intrinsic uncertainty of about ∼ 10 % (Forster et al., 2007). Note that the calculations of the GHG radiative forcing provided here are just a first-order approximation, since we use the simplified expressions of Myhre et al. (1998), 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 2 and CH 4 in our final splines, although we acknowledge that the onset of the Anthropocene is still debated (e.g. Lewis and Maslin, 2015;Steffen et al., 2016;Williams et al., 2016). P. Köhler et al.: 156 kyr smoothed history of CO 2 , CH 4 , and N 2 O 365 2 Details on the spline smoothing method The numerical code for spline smoothing is based on Enting (1987), but see also Bruno and Joos (1997) and Enting et al. (2006) 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 λ defines how much weight is given to the curvature. For a large λ, 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 λ. In other words, the smoothing of the spline depends on both the assumed errors in the data and the parameter λ.
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 c (e.g. Enting et al., 2006). The parameter λ is linked to P c as described in detail in Eq. (1) below.
Let us assume input data are t j , y j , and v j corresponding to time, value, and error (1σ ). For a given interval of the input data, an average error, v, and an average data spacing, t, can be computed. The link between the cutoff period (P c ), the data spacing ( t), and the 1σ error in the input data (v) is In the following, we prescribe P c and can calculate λ following the given relationship in Eq. (1). We choose a P c value such that it is much larger than the temporal resolution of the data, t, to avoid overfitting. However, since the choice of P c is also partially subjective, we investigate its influence on the final spline by sensitivity studies, in which P c is varied by ±50 %. One aspect of Eq. (1) is that P c depends only weakly on t.
Let us now assume we have a data set with variable data spacing, for which we would like to apply different smoothing depending on t. We proceed by modifying λ to follow the predefined individual P 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. λ is computed using Eq. (1) for the given cutoff period, average data spacing, and average error for this first interval.
-Other intervals: A modified λ = λ · s 2 , with λ taken from the reference interval, is used for other intervals, implying that for the reference interval s = 1 and λ = λ. The scaling factor s is chosen to gain the desired P c after where P c , t , and v 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 is calculated, in which the revised uncertainty v j is defined by Eq.
(2) using the cutoff-related scaling factor s. From this intermediate product, the final spline with time-dependent P 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 c1 to P c2 a transition window around the time of change, t change , exists (defined as t change ± P ct , with P ct being the smaller of P c1 and P c2 ), in which the variability of the spline transits from one cutoff period to the other and does not follow the prescribed P c exactly. The summaries of the spline calculation contained in Tables 3, 6, and 8 show the effect of this transition in a column of averaged realised cutoff periods, which are always slightly lower than the prescribed cutoff periods.
The uncertainties of the final splines are calculated from the square root of the sum of squares of three individual errors (σ = σ 2 1 + σ 2 2 + σ 2 3 ).
-P c error (σ 1 ): Mean difference from the standard spline by smoothing with cutoff periods P c which are varied by ±50 %.
-Data resolution error (σ 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 (σ 3 ): Repeated (n = 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 1. Please note that CH 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 MacFarling-Meure et al. (2006) and Rubino et al. (2013) 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.

Atmospheric CO 2
There are small offsets of a few parts per million in measured CO 2 concentration between records obtained from different ice cores (e.g. Ahn et al., 2012;Bereiter et al., 2012;Marcott et al., 2014;Bauska et al., 2015). These offsets may be related to inter-laboratory differences in the calibration or po- Error bars around the ice core data points are ±2σ . WDC data have been adjusted to reduce offsets; see text for details. In (a) the right axis contains the resulting radiative forcing R [CO 2 ] = 5.35 · ln(CO 2 /(278 ppm)) W m −2 calculated after Myhre et al. (1998). (b) Total uncertainty of the spline based on three individual error sources; see text for details. (c) Temporal resolution ( t) of the CO 2 data points underlying the spline on a log scale. Additionally, the prescribed time-dependent cutoff period P c is plotted, including its variation by ±50 %, which has been used to determine σ 1 . Table 3. Statistics of the CO 2 spline. Interval: i CO 2 ; s: scaling factor to fulfil the constraints given by the prescribed cutoff period P c ; P c : average realised cutoff period; t: mean data spacing; v: mean 1σ error -exact time framing is given by the age of the first (t start ) and last (t 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 2 for details and citations. tentially due to in situ artefact production of CO 2 in the ice archive. For a detailed discussion, see Ahn et al. (2012) and the supplement to Bereiter et al. (2012). 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 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 2. 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 Veres et al., 2013). The CO 2 record from WDC is used here on its more recent age scale WD2014 to have the timing of CO 2 and the other two GHGs consistently on the same chronology. Using WD2014 instead of the original chronology WDC06A-7 shifts the WDC CO 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 2 data compilation extends to ∼ 156 kyr BP, at which point in time well-resolved CO 2 records stop. The full CO 2 spline covering the whole time window from 2016 CE to 156 kyr BP is plotted in Fig. 1, including its uncertainty estimate (b) and the temporal resolution, t, of the compilation of data points (c). The 11-point running mean of 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, t decreases to an average of 200 years.
The CO 2 data contributing to this spline are described below (further details in Table 2):  (Fig. 2a). There is a small interpolar difference in CO 2 concentrations, with higher concentrations in the north than in the south; e.g. the 10year averages from 2006 CE to 2015 CE were 3.5 ppm   (Rubino et al., 2013) -overlap consistently with direct atmospheric measurements. We therefore take these data as our reference time series for the Common Era ( Fig. 2b) 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.
3. Data from the WDC ice core exist for the times of 11-1210 BP, or 1939-740 CE (Ahn et al., 2012;Bauska et al., 2015) and for Termination I (see point no. 5 below). These WDC data overlap with the Law Dome data (MacFarling-Meure et al., 2006;Rubino et al., 2013); however, the available high-resolution CO 2 records from different ice cores (Law Dome, WDC, EDML) show some apparent offsets (Ahn et al., 2012). Whilst the CO 2 data in all three ice cores converge on similar concentrations during the anthropogenic rise in CO 2 after 1750 CE, the WDC CO 2 concentrations are slightly higher than CO 2 in the other two ice cores prior to 1750 CE. In pre-anthropogenic times, EDC data not contained in the comparison of Ahn et al. (2012) 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 −3.13 ppm to bring them into agreement with the Law Dome CO 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 2 records within this time interval are 8 and 13 years for WDC and Law Dome, respectively. The amplitude of the CO 2 minima around 300-400 BP is controversial (Bauska et al., 2015). In our final spline, little of the large negative anomaly in CO 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. 2b). 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. 2b). 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) (Monnin et al., 2001(Monnin et al., , 2004 and further back in time (see point no. 7 below). They overlap with the Law Dome data between 350 and 1950 BP (Fig. 2b) 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. 2c). cott et al., 2014). WDC data are available for 8.8-22.9 kyr BP and are adjusted by −6.06 ppm (Fig. 2c). This difference corresponds to the duration-weighted mean offset between the WDC and EDC records during three intervals of relatively constant CO 2 (22.3-18.5 kyr BP: WDC (n = 29) 194.75 ± 2.44 ppm; EDC (n = 21) 188.22 ± 2.32 ppm; 14.5-13.0 kyr BP: WDC (n = 45) 243.02 ± 2.44 ppm; EDC (n = 9) 237.57 ± 1.42 ppm; 11.5-9.0 kyr BP: WDC (n = 36) 269.97 ± 3.67 ppm; EDC (n = 27) 264.24 ± 1.88 ppm). The intervals have been selected to minimise the influence of potential age scale differences between the two records.

Termination I is best covered by data from WDC (Mar
Only those EDC studies focusing on CO 2 measurements (Monnin et al., 2001(Monnin et al., , 2004 have been considered here, rather than those with a main focus on δ 13 C (Lourantou et al., 2010b;Schmitt et al., 2012), which have a lower precision in CO 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 2 concentration uncertainty of about 5 ppm (accuracy). The corresponding uncertainty in the radiative forcing R [CO 2 ] following a simplified expression of Myhre et al. (1998), is 0.15 or 0.09 W m −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.
6. 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 Bereiter et al., 2012). 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 2 during glacial times ( Fig. 2df) 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 (Bereiter et al., 2012), and the EDML record between 43.2 and 113.4 kyr BP (Lüthi et al., 2010;Bereiter et al., 2012). Talos Dome CO 2 data include some outliers in the interval 34-38 kyr BP that disagree with CO 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.
7. 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. 2f) -the EDC CO 2 record is used in the compilation (Schneider et al., 2013;Lourantou et al., 2010a).
For every supporting data point j a 1σ uncertainty or error v j has to be assigned in order to be able to calculate the smoothing spline (see Sect. 2 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 MacFarling-Meure et al. (2006), 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 = 2152 data points including 20 ages with duplicate entries. These duplicates are averaged (reducing n to 2132) and the assigned uncertainties based on this averaging.
To account for the variable temporal resolution of the data points (Fig. 1c) whilst preserving as much of the abrupt changes in CO 2 during Termination I as possible (e.g. Marcott et al., 2014), 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 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 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 3.
The total 1σ uncertainty of the spline is < 2 ppm on average (Fig. 1b). 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 2 record of the last 2 kyr to be used within CMIP6 (Meinshausen et al., 2017) is nearly indistinguishable from our spline across the instrumental period (Fig. 2a); however, CO 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. 2b). 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 2 values chosen as boundary conditions for several time slice experiments within PMIP4 (Ivanovic et al., 2016;Otto-Bliesner et al., 2016) can be compared with snapshots from our splines (Table 4). 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 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 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 (Monnin et al., 2001;Marcott et al., 2014), an alternative continuous CO 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 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 (Marcott et al., 2014). For the abrupt rise in CO 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 2 potentially caused by permafrost thawing during the northern hemispheric warming into the Bølling-Allerød interstadial (Köhler et al., 2014(Köhler et al., , 2015. Transient simulations investigating these abrupt jumps in CO 2 concentration should use a CO 2 times series that contains greater details than our lowfrequency spline.

Atmospheric CH 4
Our data compilation of CH 4 data and the consistently calculated CH 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 (Baumgartner et al., 2014), there remains some mismatch in the timing of the NH and the SH CH 4 records in the most recent chronology AICC2012. NH CH 4 , and consequently global CH 4 concentrations, should, according to the estimates of the interpolar difference, be larger than our SH CH 4 values. Therefore, our SH CH 4 spline represents the lower bound of global CH 4 concentration. Baumgartner et al. (2012) found that NH CH 4 was up to +4 % (+14 ppb) and up to +10 % (+60 ppb) larger than the SH CH 4 during glacial times and the Holocene, respectively. However, new and as yet unpublished results Our data compilation starts with the beginning of the year 2016 CE (−66.0 BP) and stops around 156 kyr BP to cover the same time interval as for CO 2 (Fig. 3a). The 11point running mean temporal resolution between neighbouring data points, t, is less than 100 years for most of the last 67 kyr, increasing to ∼ 700 years between 67 and 156 kyr BP (Fig. 3c). 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.  (Hansen et al., 2005;Köhler et al., 2010). 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 ( t) of the CH 4 data points underlying the spline on a log scale. Additionally, the prescribed time-dependent cutoff period P c is plotted, including its variation by ±50 %, which has been used to determine σ 1 .   Rubino et al., 2013) with an overlap of more than 2 decades with the instrumental measurements (Fig. 4a, b). The CH 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 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. 4b). However, since WDC is the more Table 6. Statistics of the CH 4 spline. Interval: i CH 4 ; s: scaling factor to fulfil the constraints given by the prescribed cutoff period P c ; P c : average realised cutoff period; t: mean data spacing; v: mean 1σ error -exact scaling factor to fulfil the constraints given by the prescribed cutoff period P c ; P c : average realised cutoff period. t: mean data spacing; v: mean 1σ error; exact time framing is given by the age of the first (t start ) and last (t 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 5 for details and citations.  Marcott et al., 2014;Buizert et al., 2015;Mitchell et al., 2013Mitchell et al., , 2011Sigl et al., 2016). Starting with the year 9821 BP, continuous CH 4 data from WDC with higher temporal resolution are now available and are used to support our spline . These continuous CH 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 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 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 4 WDC data from OSU and PSU has been identified. To account for this offset the PSU CH 4 data have been adjusted by +9.9 ppb (Supplementary Information to Rhodes et al., 2015). To date WDC CH 4 are the temporally most highly resolved data of the last glacial, and therefore they are our reference record (Fig. 4b-e). The data not only contain the wellknown abrupt CH 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. Mitchell et al., 2013;Rhodes et al., 2015).
4. We extend our SH CH 4 data compilation beyond WDC with data from EDC, spanning the period from 67 to 156 kyr BP (Loulergue et al., 2008) (Fig. 4e-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 4 record older than 156 kyr is not considered here. CH 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 (Ruth et al., 2007). However, a well-documented EDML CH 4 record is not available to date, and therefore none of the published EDML CH 4 data for this interval Schilt et al., 2010b) are considered here.
Compiled data contain 30 214 data points, among which duplicate entries exist for 39 ages. These duplicates are averaged giving n = 30 175.
The whole data set is divided into seven intervals with different assigned cutoff periods. P 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 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. 3c). More details are shown in Table 6.
The total 1σ uncertainty of our final CH 4 spline is around 3-10 ppb in the Holocene, ∼ 2 ppb in the time window supported by the continuous WDC CH 4 data (9.8-67 kyr BP), and around 10 ppb in earlier parts. During some short time windows, σ CH 4 reaches a maximum of 20 ppb (Fig. 3b). 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 σ CH 4 reached its local maxima. An abrupt jump in σ CH 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 σ 3 ) at the same time as t, and therefore P c , increased by 2 orders of magnitude (changing σ 1 ). The SH CH 4 record to be used within CMIP6 (Meinshausen et al., 2017) largely agrees with our SH CH 4 spline (Fig. 4a, b). However, during instrumental times the CMIP6 SH CH 4 record is consistently larger than our SH CH 4 spline by about 10-15 ppb, probably caused by the inclusion of different stations in the calculation of the SH CH 4 record within CMIP6, while we rely on South Pole data. Prior to the instrumental CH 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 4 record (used here but ignored within CMIP6) leads to some CH 4 variations within our SH CH 4 spline on the order 10-25 ppb that are not captured by the CMIP6 SH CH 4 record.
A comparison of our final spline with the GHG values chosen for the PMIP4 time slice experiments (Ivanovic et al., 2016;Otto-Bliesner et al., 2016) is not straightforward, since we only compile SH CH 4 data, while the PMIP4 experiments use global values. Taking the two records at face value, one finds that our SH CH 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 4 , which leads to a local minimum in SH CH 4 around 6 ka. Similarly, our SH CH 4 spline is 7 ppb higher than the global CH 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 4 record, which shows a local maximum at 21 kyr BP. A hundred years later, our SH CH 4 spline has a local minimum which is 11 ppb smaller than the global CH 4 values taken for PMIP4 (Table 4).

Atmospheric N 2 O
For the data compilation of the third GHG, N 2 O, one has to be aware that during times of high dust input, in situ production of N 2 O occurs, leading to artefacts in the paleo record (Schilt et al., 2010a). Furthermore, the precise synchronisation of Northern and Southern Hemisphere records, as already explained for CH 4 , is crucial to accurately obtain the changes in N 2 O during millennial-scale D/O events.
The compiled record starts at the beginning of the year 2016 CE (−66.0 BP) but extends back in time only until ∼134.5 kyr BP (Fig. 5a) because the ice cores on which the N 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 2 O data containing artefacts across the penultimate glacial maximum (Schilt et al., 2010a). 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 2 O of EDML and EDC, the data from EDML have a lower temporal resolution than those of EDC (Schilt et al., 2010a).
The data sets contributing to the N 2 O stack are listed below.  Myhre et al. (1998), neglecting interacting effects of CH 4 and N 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 ( t) of the N 2 O data points underlying the spline on a log scale. Additionally, the prescribed time-dependent cutoff period P c is plotted, including its variation by ±50 %, which has been used to determine σ 1 . Note that due to the long atmospheric lifetime of N 2 O, any interpolar difference can be safely neglected.  -Meure et al., 2006) and correlate well with the instrumental data in overlapping intervals (Fig. 6a,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.

Law
3. In the Holocene, N 2 O was measured at EDC  from 334 BP until 11.5 kyr BP. For the last two millennia, the EDC N 2 O data points are sparser than the Law Dome data; therefore, the EDC N 2 O data are only considered for times older than what is cov- ; s: scaling factor to fulfil the constraints given by the prescribed cutoff period P c ; P c : average realised cutoff period; t: mean data spacing; v: mean 1σ error -exact time framing is given by the age of the first (t start ) and last (t 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 7 for details and citations.  (Fig. 6b, d).
4. The most highly resolved N 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 4 (Baggenstos et al., 2017). The Taylor Glacier N 2 O record used in our spline covers the interval 9.6 to 15.8 kyr BP (Fig. 6c).
5. The last glacial interval is well resolved by N 2 O data from the NGRIP record (Flückiger et al., 2004;Schilt et al., 2010bSchilt et al., , 2013. While the NGRIP N 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 2 O data during Termination I (Fig. 6c-f). Five data points near the bedrock in the bottom part of the NGRIP records apparently have higher N 2 O values than found in ice cores from the Southern Hemisphere. These data points are rejected here, leading to the oldest NGRIP N 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. 3.2 for details), the usage of N 2 O data from NGRIP might introduce erroneous phasing between our global N 2 O record and the purely SHbased CH 4 spline, particularly during abrupt change connected to D/O events. However, N 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.
6. Additional N 2 O data going back to 134.4 kyr BP are obtained from the Talos Dome ice core (Schilt et al., 2010b) and from further data of the EDC ice core (compilation found in Schilt et al., 2010a; data source between 29.0-134.5 kyr BP: Stauffer et al., 2002;Spahni et al., 2005). Since -besides EDML, which correlates well with EDC (Schilt et al., 2010a) -these are the only N 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 2 O data seems to be systematically lower than NGRIP N 2 O data, with Taylor Glacier data in between both (Fig. 6c). We therefore believe that a mixture of all three records (N 2 O from NGRIP, EDC, and Talos Dome) most likely represents a reasonable mean global N 2 O value ( Fig. 6c-f). The relatively large difference in N 2 O from different ice cores during the last glacial times indicates that the uncertainty (accuracy) in N 2 O is probably higher than the reported measurement errors (precision) of up to 7 ppb.
The generally assigned 1σ uncertainty of each data point is 7 ppb for the Law Dome ice core (MacFarling-Meure et al., 2006). The uncertainty of individual data points in other ice cores was in general less than 7 ppb Schilt et al., 2014). 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 2 O data to n = 2344. Using an estimate of the radiative forcing of N 2 O, which neglects the interacting effects of CH 4 and N 2 O, (Myhre et al., 1998), we estimate that the 1σ error in N 2 O is related to an uncertainty in the radiative forcing of about 0.04 W m −2 , slightly larger than the uncertainty in R related to the CH 4 data. Comparing the different values of N 2 O in Talos Dome and NGRIP for same intervals reveals differences on the order of about 10 ppb (e.g. Fig. 6c-f), suggesting that the ice-core-specific values of N 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 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 ∼ 500 years (Fig. 5c). Based on this distribution of 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 c between 500 and 2000 years is prescribed. More details on the spline are found in Table 8.
The total 1σ uncertainty of the spline varies between 1 and 6 ppb (Fig. 5b). This uncertainty is mainly based on σ 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 σ 1 , the error related to the cutoff period.
If compared with the N 2 O compilation used within CMIP6 (Meinshausen et al., 2017) both approaches largely agree for instrumental times (Fig. 6a). Further back in time during the last 2 kyr, both approaches rely on the same data: the published Law Dome/Cape Grim N 2 O data (MacFarling-Meure et al., 2006). Interestingly, both time series differ by up to 6 ppb between 0.7 and 2.0 kyr BP (Fig. 6b). 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 2 O concentration than all data from Law Dome or other SH ice cores (Fig. 6b), for some unknown reason.
The N 2 O data used as starting values in the PMIP4 experiments 1850 CE, 6 kyr, 21 kyr, and 127 kyr (Ivanovic et al., 2016;Otto-Bliesner et al., 2016) 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 4).

Data availability
Data connected with this paper are available in the scientific database PANGAEA (https://doi.org/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 1σ 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 t = 1 year, including the 1σ total uncertainty.
-Corresponding radiative forcing based on the simplified Eqs.
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 2 O) and to 156 kyr BP (for CO 2 and CH 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 (Ivanovic et al., 2016) or in other future model intercomparison projects.
The resulting radiative forcing of the three GHGs calculated here with the simplified non-interacting Eqs. (3)-(5) has relative uncertainties of ∼ 10 % (Forster et al., 2007). These equations have also been used in Köhler et al. (2010) and Joos and Spahni (2008), while Schilt et al. (2010b) used a different equation for R [CO 2 ] . The latter two studies furthermore include the interacting effects of CH 4 and N 2 O (which we ignore here) but neglect the 40 % increase in CH 4 radiative forcing due to indirect effects of CH 4 (Hansen et al., 2005). 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 (Etminan et al., 2016), with differences between old and new expressions in R [CO 2 ] , R [CH 4 ] , and R [N 2 O] of less than 0.01, 0.04, and 0.02 W m −2 , respectively. While the differences in R [CO 2 ] and R [N 2 O] lie within their uncertainty bands, R [CH 4 ] is slightly higher, leading to a revised relative uncertainty of 14 % (Etminan et al., 2016). We have refrained from applying the new equations throughout the study, since the amplification in R [CH 4 ] by 40 % through indirect effects of CH 4 (Hansen et al., 2005) 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 2 , which is responsible for about two-thirds of the total radiative forcing R [GHG] during both the anthropogenic rise (Fig. 7a) and the reduction during the LGM (Fig. 7b). High-resolution variability in CH 4 (captured due to smaller cutoff periods during spline calculations than for CO 2 ) also imposes some fine-scale structure on the overall GHG radiative forcing (Fig. 7c, d); however, the dominant features are still driven by CO 2 . Figure A2. Details of differences between CO 2 of EDC (Monnin et al., 2001(Monnin et al., , 2004 and WDC (Marcott et al., 2014) during Termination I, showing how the adjustment of the WDC has been calculated. Grey areas mark the three time windows with relatively stable CO 2 from which the difference in CO 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).