Bipolar carbon and hydrogen isotope constraints on the Holocene methane budget

Abstract

Atmospheric methane concentration shows a well-known decrease over the first half of the Holocene following the Northern Hemisphere summer insolation before it started to increase again to preindustrial values. There is a debate about what caused this change in the methane concentration evolution, in particular, whether an early anthropogenic influence or natural emissions led to the reversal of the atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ concentration evolution. Here, we present new methane concentration and stable hydrogen and carbon isotope data measured on ice core samples from both Greenland and Antarctica over the Holocene. With the help of a two-box model and the full suite of ${\mathrm{CH}}_{\mathrm{4}}$ parameters, the new data allow us to quantify the total methane emissions in the Northern Hemisphere and Southern Hemisphere separately as well as their stable isotopic signatures, while interpretation of isotopic records of only one hemisphere may lead to erroneous conclusions. For the first half of the Holocene our results indicate an asynchronous decrease in Northern Hemisphere and Southern Hemisphere${\mathrm{CH}}_{\mathrm{4}}$ emissions by more than 30 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$ in total, accompanied by a drop in the northern carbon isotopic source signature of about $-\mathrm{3}$ ‰. This cannot be explained by a change in the source mix alone but requires shifts in the isotopic signature of the sources themselves caused by changes in the precursor material for the methane production. In the second half of the Holocene, global ${\mathrm{CH}}_{\mathrm{4}}$ emissions increased by about 30 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$, while preindustrial isotopic emission signatures remained more or less constant. However, our results show that this early increase in methane emissions took place in the Southern Hemisphere, while Northern Hemisphere emissions started to increase only about 2000 years ago. Accordingly, natural emissions in the southern tropics appear to be the main cause of the ${\mathrm{CH}}_{\mathrm{4}}$ increase starting 5000 years before present, not supporting an early anthropogenic influence on the global methane budget by East Asian land use changes.

Introduction

Atmospheric methane (${\mathrm{CH}}_{\mathrm{4}}$) is a potent greenhouse gas and its concentrations are strongly coupled to the Earth's climate system. Due to the human influence on the Earth system, the ${\mathrm{CH}}_{\mathrm{4}}$ concentration ([${\mathrm{CH}}_{\mathrm{4}}$]) in the atmosphere has increased by a factor of 2.5 (relative to the preindustrial level) over the last centuries and contributes significantly to the human-induced radiative forcing (Etheridge et al., 1998; Dlugokencky et al., 2005). Today both anthropogenic ${\mathrm{CH}}_{\mathrm{4}}$ sources (rice agriculture, livestock, fossil fuel production, anthropogenic biomass burning, and landfills) and natural ${\mathrm{CH}}_{\mathrm{4}}$ sources (natural wetlands, wildfires, geologic emissions, wild animals (including termites), and marine ${\mathrm{CH}}_{\mathrm{4}}$ hydrates) contribute to the global ${\mathrm{CH}}_{\mathrm{4}}$ emissions (Kirschke et al., 2013). The main mechanism removing ${\mathrm{CH}}_{\mathrm{4}}$ from the atmosphere is the chemical reaction of ${\mathrm{CH}}_{\mathrm{4}}$ with OH radicals in the troposphere. Together with the other sink processes (such as stratospheric loss, reaction with Cl radicals in the marine boundary layer, and soil uptake) the OH sink determines the recent atmospheric lifetime of ${\mathrm{CH}}_{\mathrm{4}}$, which is $\mathrm{9.1}\pm \mathrm{0.9}$ years (Prather et al., 2012).

${\mathrm{CH}}_{\mathrm{4}}$ concentration data compilation. Discrete [${\mathrm{CH}}_{\mathrm{4}}$] data over the Holocene used for the ${\mathrm{CH}}_{\mathrm{4}}$ data compilation on synchronised gas age scales and corrected for offsets among data sets (references are listed in Table 1). The $\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}$-${\mathrm{H}}_{\mathrm{2}}$O data from the layer-counted NGRIP ice core have been used as an absolute time marker for fast climatic changes (e.g. the Younger Dryas–Holocene transition and the 8.2 ka event) assuming an essentially synchronous change in rapid temperature variations and atmospheric ${\mathrm{CH}}_{\mathrm{4}}$. The markers at the bottom of the plot indicate the tie points of the synchronised age scales (see Table 3).

[Figure omitted. See PDF]

Despite the fact that the total ${\mathrm{CH}}_{\mathrm{4}}$ emissions today are well known, there is still considerable debate about the individual contributions of (in particular small) ${\mathrm{CH}}_{\mathrm{4}}$ sources and their variability, e.g. evidenced by the mismatch of bottom-up and top-down estimates of the total ${\mathrm{CH}}_{\mathrm{4}}$ emissions (Kirschke et al., 2013; Crill and Thornton, 2017; Saunois et al., 2016). Moreover, a long-standing debate exists whether long-term Holocene ${\mathrm{CH}}_{\mathrm{4}}$ trends are naturally caused by wetland emissions (Singarayer et al., 2011) or whether an anthropogenic influence started already 4000–5000 years ago with early human land use changes (Ruddiman, 2003) (see below). To better assess the human influence on the global methane cycle, robust estimates of the natural (pre-anthropogenic) ${\mathrm{CH}}_{\mathrm{4}}$ emissions and their climate-coupled variations are required. Only polar ice cores can provide this information as they represent a direct, albeit (through the bubble enclosure process in the firn) low-pass-filtered, archive of trapped air and therefore offer the possibility to investigate the atmospheric composition of the past. With such ice core measurements the millennial to centennial evolution of the atmospheric [${\mathrm{CH}}_{\mathrm{4}}$] has been determined back to 800 000 years before present (ka BP), where present refers to the year 1950, using the Antarctic EPICA Dome C (EDC) ice core showing glacial–interglacial variations in atmospheric [${\mathrm{CH}}_{\mathrm{4}}$] by a factor of about 2 (Loulergue et al., 2008). More important for the assessment of the human influence on ${\mathrm{CH}}_{\mathrm{4}}$ levels are their variations during the Holocene, i.e. our current interglacial that started about 10 000 years ago. In this time period [${\mathrm{CH}}_{\mathrm{4}}$] decreased by about 100 parts per billion (ppb) until about 5000 years ago (Fig. 1) after an early maximum during the Preboreal period (about 11.5–10 ka BP) (Blunier et al., 1995; Brook et al., 1996; Flückiger et al., 2002; Schilt et al., 2010). ${\mathrm{CH}}_{\mathrm{4}}$ values started to slowly increase again in the mid-Holocene until the fast anthropogenic increase resulting from industrialisation started around 200 years ago (MacFarling Meure et al., 2006). High-resolution ${\mathrm{CH}}_{\mathrm{4}}$ measurements on the West Antarctic Ice Sheet (WAIS) and Greenland GISP2 (Greenland Ice Sheet Project 2) ice cores (Mitchell et al., 2013) over the last 2800 years also show centennial ${\mathrm{CH}}_{\mathrm{4}}$ variations of 20–30 ppb superimposed on the long-term increase, reflecting natural variability in ${\mathrm{CH}}_{\mathrm{4}}$ emissions.

With the trend reversal in the mid-Holocene leading to a “bowl shape” (Blunier et al., 1995) of [${\mathrm{CH}}_{\mathrm{4}}$], the Holocene differs from many previous interglacials, where the [${\mathrm{CH}}_{\mathrm{4}}$] steadily dropped to glacial levels concurrently with the Northern Hemisphere summer insolation. One hypothesis postulates that ${\mathrm{CH}}_{\mathrm{4}}$ emissions caused by early human land use were responsible for the ${\mathrm{CH}}_{\mathrm{4}}$ turnaround at about 5 ka BP (Ruddiman and Thomson, 2001; Ruddiman, 2003). According to these authors, early farming activities, mainly rice agriculture in the eastern part of China, led to significant releases of ${\mathrm{CH}}_{\mathrm{4}}$ long before the industrial era (Ruddiman et al., 2008) and thus the authors propose that the bowl was shaped by human activities. Another explanation supported by a model-based study explains the late Holocene ${\mathrm{CH}}_{\mathrm{4}}$ rise by calling upon increased emissions from southern tropical wetlands due to an insolation-driven strengthening of the monsoon in western Amazonia (Singarayer et al., 2011). Accordingly, this model is able to produce a global [${\mathrm{CH}}_{\mathrm{4}}$] rise in the Holocene driven by changes in the orbital forcing only. Note that due to the different orbital configuration during the last interglacial (MIS5.5), this model is also able to reproduce the steadily declining [${\mathrm{CH}}_{\mathrm{4}}$] during MIS5.5, when strongly declining boreal ${\mathrm{CH}}_{\mathrm{4}}$ emissions overcompensate for increasing southern tropical wetland emissions.

An important point is that the two scenarios explain the Holocene ${\mathrm{CH}}_{\mathrm{4}}$ rise with additional emissions taking place in different geographic regions. By measuring the [${\mathrm{CH}}_{\mathrm{4}}$] on ice cores from both polar regions (Greenland and Antarctica), the interpolar difference (IPD) can be calculated. This quantity can be used to draw conclusions about the hemispheric distribution of the ${\mathrm{CH}}_{\mathrm{4}}$ emissions (Brook et al., 2000; Chappellaz et al., 1997; Mitchell et al., 2013; Yang et al., 2017). Another tool to investigate the processes involved in the ${\mathrm{CH}}_{\mathrm{4}}$ cycle are measurements of the stable isotope signature in ice core ${\mathrm{CH}}_{\mathrm{4}}$. The ${\mathrm{CH}}_{\mathrm{4}}$ released by the various sources is associated with typical hydrogen and carbon isotopic signatures $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ (D representing deuterium or ${}^{\mathrm{2}}$H) (Quay et al., 1999; Walter et al., 2008; Whiticar, 1993; Whiticar and Schaefer, 2007) and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ (Whiticar, 1993; Quay et al., 1999; Whiticar and Schaefer, 2007; Etiope et al., 2008; Walter et al., 2008). In addition, the sink processes lead to a strong fractionation (Snover and Quay, 2000; Feilberg et al., 2005; Whiticar and Schaefer, 2007; Levine et al., 2011b), which leaves atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ strongly enriched in heavy isotopes relative to the isotopic emission signatures. Changes of the relative contributions in sources lead to an alteration of the atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ isotopic composition. ${\mathrm{CH}}_{\mathrm{4}}$ isotope measurements in ice cores therefore allow us to deduce the role the different processes related to the ${\mathrm{CH}}_{\mathrm{4}}$ cycle played in the past. By measuring both $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathit{\delta }}^{\mathrm{13}}$C-${\mathrm{CH}}_{\mathrm{4}}$ from the Greenland GISP2 ice core over the Holocene, Sowers (2010) was able to provide first evidence of possible changes in the ${\mathrm{CH}}_{\mathrm{4}}$ precursor material due to a shift in the C3 to C4 plant ratio but could not unambiguously answer the question about the Holocene ${\mathrm{CH}}_{\mathrm{4}}$ anomaly.

In this study we combined for the first time comprehensive information on the ${\mathrm{CH}}_{\mathrm{4}}$ IPD and measurements of the bipolar carbon and hydrogen stable isotopic signature of ${\mathrm{CH}}_{\mathrm{4}}$ by measuring all three parameters – [${\mathrm{CH}}_{\mathrm{4}}$], $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$, and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ – on ice samples from both polar regions with high analytical precision and improved resolution. Furthermore, we use a two-box model approach, which allows us to deconvolve the atmospheric signal and to calculate the total emissions and their isotopic signatures in each hemisphere.

As mentioned above, changes in the sinks leave an imprint on the isotopic composition of ${\mathrm{CH}}_{\mathrm{4}}$ in the atmosphere. However, while minor changes in the relative importance of the sink processes cannot be ruled out, in general, atmospheric chemistry models point to small changes in total lifetime (Kaplan et al., 2006; Levine et al., 2011a) and we also assume the relative contribution of individual sinks to be rather constant over the Holocene when relatively stable climate conditions prevailed. Therefore, the discussion of this study centres around ${\mathrm{CH}}_{\mathrm{4}}$ source processes and hence emissions of different ${\mathrm{CH}}_{\mathrm{4}}$ source categories.

Method

Measurements

The stable isotope data ($\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$) presented in this study were measured at the University of Bern on samples from the NGRIP (North Greenland Ice Core Project) ice core from Greenland. Antarctic samples from the EDML (EPICA Dronning Maud Land) and the TALDICE (TALos Dome Ice CorE) ice cores were analysed for $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$, respectively (Bock et al., 2017). Note that for both parameters the samples from Greenland and Antarctica were measured during the same time interval with the same measurement system by the same operators and in randomised order. This procedure is crucial to avoid any systematic error (e.g. due to a long-term drift of the system) and ensures unbiased results in our source inversion.

The two isotopic parameters were measured at the University of Bern on two independent measurement systems designed for this specific purpose as described in detail in previous publications (Bock et al., 2010, 2014; Schmitt et al., 2014). They follow the same principal procedure: the air is (i) extracted from the ice sample by purging the meltwater with He; (ii) water vapour, ${\mathrm{CO}}_{\mathrm{2}}$, and the major air compounds (${\mathrm{N}}_{\mathrm{2}}$, ${\mathrm{O}}_{\mathrm{2}}$, and Ar) are removed from the gas sample; (iii) ${\mathrm{CH}}_{\mathrm{4}}$ is separated from other trace gases using gas chromatography (GC); (iv) ${\mathrm{CH}}_{\mathrm{4}}$ is pyrolysed to ${\mathrm{H}}_{\mathrm{2}}$ or combusted to ${\mathrm{CO}}_{\mathrm{2}}$ in the $\mathit{\delta }$D and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ systems, respectively; and (v) the isotopic composition of ${\mathrm{H}}_{\mathrm{2}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ is measured using an isotope-ratio mass spectrometer (IRMS). For both systems effects of interfering masses, specifically krypton, in the mass spectrometer are avoided by technical adaptations (Schmitt et al., 2013). In the case of $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$, the focussed ${\mathrm{H}}_{\mathrm{2}}$ peak is separated from other components using a second GC separation (post-pyrolysis trapping) (Bock et al., 2014). In the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ system, the ${\mathrm{CH}}_{\mathrm{4}}$-derived ${\mathrm{CO}}_{\mathrm{2}}$ is captured and delayed using a fused silica trap (at $-\mathrm{196}$ ${}^{\circ }$C) after the combustion oven and thereby separated from krypton (Schmitt et al., 2014). Note that many $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ values based on such a gas chromatography–mass spectrometry method reported in the literature are not corrected for this effect, leading to significant offsets, especially for low ${\mathrm{CH}}_{\mathrm{4}}$ concentrations.

The $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ values are determined relative to our primary air standard Air Controlé, a recent clean air ([${\mathrm{CH}}_{\mathrm{4}}$] $=\mathrm{1971}\pm \mathrm{7}$ ppb), which has been cross-referenced to $-\mathrm{93.6}\pm \mathrm{2.2}$ ‰ with respect to Vienna Standard Mean Ocean Water (VSMOW) using Alert 2002/11, an air collected at Alert Station, Canada (Bock et al., 2014). The latter was referenced by Poss (2003) to be $-\mathrm{82.15}\pm \mathrm{0.3}$ ‰ with respect to VSMOW, which anchors our $\mathit{\delta }$D values to the isotope scale of the University of Heidelberg (Bock et al., 2010, 2014). Note that so far no internationally uniform $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ standardisation scale has been established, leading to systematic inter-laboratory offsets as described in Umezawa et al. (2018). The $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ measurement precision ($\mathrm{1}\mathit{\sigma }$) of 2.3 ‰ was derived by calculating the pooled standard deviation of multiple NGRIP replicate (vertically neighbouring samples) measurements. The reference gas for the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ measurements is Boulder, an ambient air diluted with ${\mathrm{CH}}_{\mathrm{4}}$-free air ([${\mathrm{CH}}_{\mathrm{4}}$] $=\mathrm{1508.18}\pm \mathrm{0.17}$ ppb) calibrated by the National Oceanic and Atmospheric Administration (NOAA) to a ${\mathit{\delta }}^{\mathrm{13}}$C-${\mathrm{CH}}_{\mathrm{4}}$ value of $-\mathrm{47.34}\pm \mathrm{0.02}$ ‰ with respect to Vienna Pee Dee Belemnite (VPDB) (Schmitt et al., 2014; Umezawa et al., 2018). The $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ measurement precision ($\mathrm{1}\mathit{\sigma }$) derived from replicate measurements is better than 0.15 ‰ for ice core samples.

Monte Carlo spline approximation of millennial ${\mathrm{CH}}_{\mathrm{4}}$ changes

To compare values from different ice core data sets (with different and variable temporal resolutions and measurement precisions) interpolation of the data to a common age scale is needed. Additionally, the records have to be low-pass filtered to the same (millennial) timescales that can be resolved in the isotope records. Therefore, in this study smoothing splines following Enting (1987) are used to calculate a continuous low-pass-filtered record from the discretely measured values. The stiffness of the spline curve is characterised by its cut-off period. A Python (http://www.python.org/, last access: 11 October 2018) routine was developed to calculate splines in a Monte Carlo manner, in which the data points are varied within a normal distribution within their standard deviation defined by the individual uncertainty of the measured values. This procedure allows for the calculation of the average of the splines representing the best-guess mean evolution of the data on the timescale resolved by the spline. Furthermore, using the standard deviation of all splines, the spread of the individual splines can be translated to an uncertainty of the mean spline. This spline calculation is used for various purposes in this study. The choice of the cut-off period has a large influence on the outcome of the spline and is therefore always documented wherever it is selected manually. For some applications (e.g. in the routines described in Sect. 2.3.2 and 2.3.3), however, we apply an objective criterion to define an appropriate cut-off period for a data set: wherever the cut-off period is not specified explicitly, it is set to 8 times the median of the data resolution.

Data sets used for ${\mathrm{CH}}_{\mathrm{4}}$ data compilation. The ${\mathrm{CH}}_{\mathrm{4}}$ offset correction is added to the data to correct for inter-laboratory and time-period-related offsets (see Sect. 2.3.2). Additional data sets used for the time matching (e.g. covering the rapid ${\mathrm{CH}}_{\mathrm{4}}$ increase after 1850 CE or the Younger Dryas–Holocene transition) are listed.

Data used for ${\mathrm{CH}}_{\mathrm{4}}$ compilation
Data set name	Data type	Reference	${\mathrm{CH}}_{\mathrm{4}}$ offset (ppb)
NGRIP_Be18	NGRIP ice core	This study	0
GRIP_Fl01	GRIP ice core	GRIP data Flückiger (2001, unpublished)${}^{\ast }$	17.69
GISP2_Br96	GISP2 ice core	Brook et al. (1996)	34.13
GISP2_Br09	GISP2 ice core	Brook (2009)	$-\mathrm{1.81}$
GISP2_Mi13	GISP2 ice core	Mitchell et al. (2013)	$-\mathrm{6.15}$
TALDICE_Bo17	TALDICE ice core	Bock et al. (2017)	0
TALDICE_Sc10	TALDICE ice core	Schilt et al. (2010)	$-\mathrm{2.09}$
EDML_Sc10	EDML ice core	Schilt et al. (2010)	$-\mathrm{4.63}$
WAIS_Mi13a	WAIS (WDC05A) ice core	Mitchell et al. (2013)	$-\mathrm{6.15}$
WAIS_Mi13b	WAIS (WDC06A) ice core	Mitchell et al. (2013)	$-\mathrm{6.15}$
EDC_Fl02	EDC ice core	Flückiger et al. (2002)	4.23
Additional data used for time matching
Data set name	Data type	Reference
NGRIP_d18O	NGRIP ice core $\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}$(${\mathrm{H}}_{\mathrm{2}}$O)	Rasmussen et al. (2006); Vinther et al. (2006)
NGRIP_Ba12	NGRIP ice core	Baumgartner et al. (2014)
NEEM_Bu12	NEEM firn model	Buizert et al. (2012)
NOAA_Alert	Alert flasks	Dlugokencky et al. (2017)
LawDome_Ma06	Law Dome firn and ice core	MacFarling Meure et al. (2006)
CapeGrim_csiro	Cape Grim flasks	CSIRO (27 November 2017)

${}^{\ast }$ GRIP data compilation performed by Jacqueline Flückiger in 2001 containing data from various publications, e.g. Blunier et al. (1995) and Chappellaz et al. (1993, 1997).

If the sampling rate in some parts of a data set is much higher than in others, this part has a larger impact on the spline than other more coarsely resolved sections. For short periods of large signal change (e.g. the 8.2 ka event for ${\mathrm{CH}}_{\mathrm{4}}$) such an oversampling relative to other parts of the data set may lead to significant artefacts in the calculated spline. Wherever this might be an issue, an additional step was introduced into the spline calculation routine, in which the data were down-sampled in a randomised way within the Monte Carlo loop (i.e. for each individual spline) to end up with a subset of a prescribed number of data points, for which the chance of data points to get picked is inversely proportional to the local data density. The down-sampled data are still unequally distributed, but on average the data density is similar over the whole time period studied.

Data correction and compilation

For further analysis of our new data, we complemented our [${\mathrm{CH}}_{\mathrm{4}}$] records with data from other studies from the GRIP (Greenland Ice Core Project), GISP2, TALDICE, EDML, WAIS (West Antarctic Ice Sheet), and EDC ice cores shown in Table 1. However, combining different data sets introduces various challenges. In the following paragraphs we describe how these difficulties have been addressed.

Timescale synchronisation

All ice cores have individually determined gas age scales, which allow us to assign an age to each sample. Since we compare values from different ice cores (and from both hemispheres), small deviations in the age calculations have significant consequences for the interpretation (e.g. artefacts in the box model inversion; see Sect. 2.4). In an attempt to synchronise the different data sets, we linearly stretched and compressed the individual age scales between manually chosen tie points in three steps.

• Wherever possible, we tied the data sets of each of the cores to the same reference age scale. For ${\mathrm{CH}}_{\mathrm{4}}$ data sets that reach far enough to the present, the fast ${\mathrm{CH}}_{\mathrm{4}}$ increase of the last centuries was used to synchronise the timescales. As a reference we used [${\mathrm{CH}}_{\mathrm{4}}$] data from the Law Dome ice core (MacFarling Meure et al., 2006) and an atmospheric [${\mathrm{CH}}_{\mathrm{4}}$] history derived from NEEM firn air (Buizert et al., 2012) representing the ${\mathrm{CH}}_{\mathrm{4}}$ evolution in Antarctica and in Greenland, respectively. Both records cover the entire anthropogenic ${\mathrm{CH}}_{\mathrm{4}}$ increase and overlap with recent atmospheric measurements (CSIRO, 2017; Dlugokencky et al., 2017). For the early Holocene the layer-counted NGRIP $\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}$-${\mathrm{H}}_{\mathrm{2}}$O record was used as reference timescale. Where the ${\mathrm{CH}}_{\mathrm{4}}$ data show strong signals such as at the 8.2 ka event and the Younger Dryas–Holocene transition, these were synchronised to the corresponding $\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}$ signal, assuming that there is no significant lag between the temperature changes recorded in $\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}$ and the response of ${\mathrm{CH}}_{\mathrm{4}}$ source to this change (Baumgartner et al., 2014).

• With the use of splines of different cut-off periods and their time derivatives we were looking for typical features and signals in the [${\mathrm{CH}}_{\mathrm{4}}$] evolution of the different data sets within the same hemisphere to define additional tie points and synchronise the data sets for each hemisphere. However, within each hemisphere the temporal alignments of all the distinguishable ${\mathrm{CH}}_{\mathrm{4}}$ features were already within the limitation of the data resolution. Therefore, no additional tie points could be defined at this stage, within each hemisphere.

• In the third synchronisation step we aligned the data of both hemispheres after applying the offset corrections (see Sect. 2.3.2). Additional time corrections (for the Southern Hemisphere data) could be defined by looking at the fine structure of the [${\mathrm{CH}}_{\mathrm{4}}$] evolution (e.g. the mid-Holocene [${\mathrm{CH}}_{\mathrm{4}}$] minimum (around 5400 a BP) or the local [${\mathrm{CH}}_{\mathrm{4}}$] maximum before the early Holocene decrease (around 9850 a BP)). Over the late Holocene period Mitchell et al. (2013) provided already synchronised data for the GISP2 and the WAIS ice core based on decadal to centennial variability; we retained their timescale over this period.

Table 2 lists the original age scales of the different ice cores. The tie points (original age and age correction) used to achieve an optimal temporal synchronisation among the ${\mathrm{CH}}_{\mathrm{4}}$ data sets are shown in Table 3. Temporal shifts to data with ages younger than the first tie point and older than the last tie point were chosen to be identical to those of the first and last tie points, respectively. Please note that for the ice cores, where the primary age scales provide an age uncertainty, all but three tie point corrections (GRIP: 273 years; TALDICE: 127 years and 6042 years) lie within the $\mathrm{1}\mathit{\sigma }$ range of the uncertainty given for the individual age scales.

Ice cores and age scales. Short names of the ice cores and initial gas age scales used for the ${\mathrm{CH}}_{\mathrm{4}}$ compilation.

Ice core	Initial age scale	References
NGRIP	AICC2012 ice core chronology	Bazin et al. (2013); Veres et al. (2013)
GRIP	GICC05	Rasmussen et al. (2014); Seierstad et al. (2014)
GISP2	Mitchell, Meese/Sowers	Meese et al. (1994); Mitchell et al. (2013)
TALDICE	AICC2012 ice core chronology	Bazin et al. (2013); Veres et al. (2013)
EDML	AICC2012 ice core chronology	Bazin et al. (2013); Veres et al. (2013)
EDC	AICC2012 ice core chronology	Bazin et al. (2013); Veres et al. (2013)
WAIS	Mitchell, Buizert WD2004	Buizert et al. (2015); Mitchell et al. (2013)

Time synchronisation. Tie points of the reference timescale for the individual ice cores (left column, in yr BP) and time corrections (right column, in years) applied with respect to the core-dependent timescale listed in Table 2. The correction values are linearly interpolated among the tie points and held constant before the earliest and after the latest tie points.

NGRIP		GRIP		GISP2		TALDICE		EDML		EDC		WAIS
294	0	273	$-\mathrm{40}$	2752	$-\mathrm{5}$	127	$-\mathrm{134}$	2747	0	349	$-\mathrm{225}$	2752	$-\mathrm{5}$
11 800	$-\mathrm{150}$	8157	30	8172	$-\mathrm{15}$	2835	$-\mathrm{88}$	3754	39	2959	$-\mathrm{212}$	3753	40
		11 606	34	11 771	$-\mathrm{165}$	3825	$-\mathrm{32}$	4538	63	3961	$-\mathrm{168}$	4532	69
						4596	5	5279	130	4741	$-\mathrm{140}$	5269	140
						5325	84	6008	215	5478	$-\mathrm{69}$	5993	230
						6042	181	6361	129	6204	19	6344	146
						6389	101	8107	36	6555	$-\mathrm{65}$	8079	64
						8106	37	8164	22	8292	$-\mathrm{149}$	8228	28
						8252	4	8254	2	8353	$-\mathrm{164}$	9869	$-\mathrm{13}$
						9876	$-\mathrm{20}$	9875	$-\mathrm{19}$	8438	$-\mathrm{182}$	10 391	110
						10 393	108	10 391	110	9397	$-\mathrm{179}$	11 591	49
						11 580	60	11 576	64	9989	$-\mathrm{133}$
										10 462	39
										11 547	93

Offset correction

Another issue that can complicate compiling multiple data sets is concentration offsets. The ${\mathrm{CH}}_{\mathrm{4}}$ data sets we use to complement our new data were measured with different methods, in different labs, and at different times over the last decades. If we use the information about the relative evolution of the methane concentration documented by other measurement series, we have to be aware that the absolute level might show an offset from our data due to calibration and extraction issues. Therefore, we correct all data sets for such potential offsets. To do so, we calculated residuals defined as the difference between the spline values of the ${\mathrm{CH}}_{\mathrm{4}}$ data from other labs evaluated at the age of our data and the ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratios of our data. A least-squares optimisation routine was used to find the offsets for the spline values of the different data sets. The offset corrections applied to the data are shown in Table 1. The GISP2 and WAIS data sets by Mitchell et al. (2013) are treated as a single data set, as they were measured in the same lab and during the same measurement campaign, and their difference (Greenland minus Antarctica) represents the true IPD. In Fig. 1 the synchronised and offset-corrected ${\mathrm{CH}}_{\mathrm{4}}$ time series as used in this study are shown.

Outlier detection

The GISP2 data sets from Brook et al. (1996) and Brook (2009) show an increased incidence of anomalously high [${\mathrm{CH}}_{\mathrm{4}}$] values (relative to other data sets) in the mid-Holocene. As this ice is from the brittle ice zone (the zone where air bubbles and clathrates coexist in the ice, and which is prone to bad core quality due to the pressure and temperature relaxation process after the core was brought to the surface) this is probably due to infiltration of modern air into the ice samples in (micro)cracks in the ice after core retrieval (Gow et al., 1997; Neff, 2014). To avoid artificial elevation of the splined signal, which would result in a significant increase in the IPD, we identified and removed such outliers within the depth interval of 650–1400 m of the GISP2 ice core (equivalent to the time window 2.7–8.1 ka BP). To do so we used an iterative elimination approach. In each step, first the two GISP2 data sets were corrected for offsets relative to each other (Brook et al., 1996, was shifted to fit with Brook, 2009, by 36 ppb) and combined to one data set. Thereafter, for each GISP2 data point the difference between its value and the spline of all other GISP2 data in the brittle zone was calculated and divided by the $\mathrm{1}\mathit{\sigma }$ uncertainty of the Monte Carlo spline at its age. If the largest (positive) relative residual exceeded the threshold value of 4, the corresponding data point was removed. This procedure was repeated until no outliers could be found any more. Table 4 shows the data points that were identified and eliminated through this method (also shown in Fig. 1). In total 15 out of 95 data points (within this depth range) were removed, reflecting the bad core quality in this depth interval. Note that in this outlier detection routine we make the assumption that the [${\mathrm{CH}}_{\mathrm{4}}$] evolution in the atmosphere is relatively smooth. As we focus only on millennial changes in the data here, this assumption does not affect our conclusions. In the following discussion we assume that the corrected GISP2 data and the ${\mathrm{CH}}_{\mathrm{4}}$ data from the GRIP and NGRIP ice cores over this time interval are not affected anymore by modern air intrusion in the brittle ice. If this assumption does not hold, the true Northern Hemisphere ${\mathrm{CH}}_{\mathrm{4}}$ concentration would be lower and so would Northern Hemisphere emissions at that time.

GISP2 outliers. Data points that were removed from the two GISP2 data sets by Brook et al. (2009) and Brook (1996) due to elevated ${\mathrm{CH}}_{\mathrm{4}}$ values probably due to contamination related to the ice sample quality in the brittle zone. Please note that the age and the concentration values are corrected for the age shift and the ${\mathrm{CH}}_{\mathrm{4}}$ correction, respectively.

Removed data
Depth (m)	Age (a BP)	${\mathrm{CH}}_{\mathrm{4}}$ (ppb)	Data set
655.1	2731.9	665.0	GISP2_Br96
736.1	3199.1	647.0	GISP2_Br96
819.7	3694.8	680.7	GISP2_Br09
932.7	4411.4	663.6	GISP2_Br09
992.3	4804.6	679.2	GISP2_Br09
1019.0	4989.5	643.8	GISP2_Br09
1059.0	5272.2	598.0	GISP2_Br96
1240.0	6694.7	649.6	GISP2_Br09
1253.9	6810.3	654.0	GISP2_Br96
1274.8	6978.6	644.0	GISP2_Br96
1392.0	8011.9	684.0	GISP2_Br09

Isotope data corrections

All stable isotope data have been corrected for gravitational enrichment. This process takes place in the firn column and leads to an enrichment of the heavier isotopes at the depth at which the air bubbles are finally closed off and the gas is trapped in the ice (Schwander, 1996). Using ${\mathit{\delta }}^{\mathrm{136}}$Xe data, which come as a side product during the ${\mathit{\delta }}^{\mathrm{13}}$C-${\mathrm{CH}}_{\mathrm{4}}$ measurement of the Bern $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ system, and additional published $\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}$-${\mathrm{N}}_{\mathrm{2}}$ data (Landais et al., 2006; Capron et al., 2013; Eggleston et al., 2016), we experimentally derived this gravitational enrichment and directly correct the ${\mathrm{CH}}_{\mathrm{4}}$ stable isotope data, as this gravitational enrichment is only dependent on the mass difference between the two isotopes considered. The correction for the TALDICE $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ values decreases from the early Holocene to the present from 0.39 ‰ to 0.28 ‰. The NGRIP ${\mathit{\delta }}^{\mathrm{136}}$Xe data do not show a systematic change. Thus, a constant value of 0.32 ‰ is subtracted from the NGRIP isotope values. Only very few $\mathit{\delta }{}^{\mathrm{15}}\mathrm{N}$-${\mathrm{N}}_{\mathrm{2}}$ data points in the early Holocene were available to correct the EDML $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ record. Since in the case of $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ the measurement error and the observed signal are much larger than the gravitational enrichment, it is adequate to assume that this correction is also constant over time without any consequences for the interpretation of the data. For the EDML $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ record a correction of 0.46 ‰ was applied.

Due to the relatively constant climatic conditions and slow [${\mathrm{CH}}_{\mathrm{4}}$] changes occurring over the Holocene, thermal diffusion (Severinghaus et al., 1998) and diffusive fractionation (Buizert et al., 2013) – other processes potentially leading to an artefact in the stable isotope signature – are not considered in this study.

Box model inversion

Our box model inversion is a tool to disentangle the different processes involved in the atmospheric methane cycle. More precisely, it allows us to obtain quantitative information about the emission processes that led to the tropospheric ${\mathrm{CH}}_{\mathrm{4}}$ concentration changes recorded in the ice cores in the past. In other similar studies different approaches of varying complexity and with different boundary conditions have been presented (Fischer et al., 2008; Baumgartner et al., 2012; Mitchell et al., 2013). For this study we use a two-box model in which the boxes represent the northern and southern tropospheric hemispheres.

The following six inversion equations allow us to use the six measured parameters (the tropospheric [${\mathrm{CH}}_{\mathrm{4}}$] ($c_x$) and the two stable isotope ratios (${}^i{R}_x\left)\right)$ for each hemisphere to calculate the emissions ($E_x$) and their isotope ratios (${}^i{R}_{E_x}$) for both boxes analytically (with $x$ and $y$ representing north and south and $i$ the mass numbers $\mathrm{2}$ or $\mathrm{13}$ for deuterium and ${}^{\mathrm{13}}$C): $$\begin{array}{}& {\displaystyle }& {\displaystyle }{E}_x=m^{\ast }\cdot \left({\displaystyle \frac{r_x\cdot c_x}{{\mathit{\tau }}_x}}+{\displaystyle \frac{c_x}{\mathrm{2}\cdot \mathit{\theta }}}-{\displaystyle \frac{c_y}{\mathrm{2}\cdot \mathit{\theta }}}\right),& {\displaystyle }& {}^i{R}_{E_x}={\displaystyle \frac{m^{\ast }}{E_x}}\cdot \left({\displaystyle \frac{r_x\cdot c_x}{{\mathit{\tau }}_x}}{\cdot }^i{R}_x{\cdot }^i{\mathit{\alpha }}_x+{\displaystyle \frac{c_x}{\mathrm{2}\cdot \mathit{\theta }}}{\cdot }^i{R}_x-{\displaystyle \frac{c_y}{\mathrm{2}\cdot \mathit{\theta }}}{\cdot }^i{R}_y\right).\end{array}$$

The first term in Eqs. (1) and (2) represents the loss due to the ${\mathrm{CH}}_{\mathrm{4}}$ sinks, and the other two terms represent the interpolar mixing (out of and into the box, with $\mathit{\theta }$ being the mixing time between the poles as explained in more detail below). The parameter $m^{\ast }=\frac{m_{\mathrm{0}}}{c_{\mathrm{0}}}$ translates atmospheric concentrations (ppb) to total atmospheric inventories (Tg) using the mean atmospheric concentration value $c_{\mathrm{0}}=\mathrm{1650}$ ppb of the year 1987. Here we use the corresponding global ${\mathrm{CH}}_{\mathrm{4}}$ burden $m_{\mathrm{0}}=\mathrm{4800}$ Tg by Steele et al. (1992) that was also used in previous ice core studies (Fischer et al., 2008; Baumgartner et al., 2012). Note that this value is 7 %–10 % higher than more recent estimates (Dalsøren et al., 2016) and accordingly the derived absolute emissions may be 7 %–10 % overestimated. However, we focus our interpretation on relative emission changes, which are not affected by this scaling factor. The transport between the boxes is quantified by the interpolar exchange time $\mathit{\theta }=\mathrm{1.56}$ yr, which is determined through a model calibration using recent sulfur hexafluoride (SF${}_{\mathrm{6}}$) data (described in Sect. 2.4.1).

Values used in our box model inversion. Relative strengths, isotopic fractionations (${}^{\mathrm{2}}\mathit{\epsilon }$ and ${}^{\mathrm{13}}\mathit{\epsilon })$ (Cantrell et al., 1990; Tyler et al., 1994; Brenninkmeijer et al., 1995; Irion et al., 1996; Gierczak et al.,1997; Quay et al., 1999; Feilberg et al., 2005), the north–south ratio of the individual sink processes (ratio $n/s$), and the parameters for each hemisphere used in our ${\mathrm{CH}}_{\mathrm{4}}$ box model inversion.

Sink processes
	Rel. strength	${}^{\mathrm{2}}\mathit{\epsilon }\left(\mathrm{‰}\right)$	${}^{\mathrm{13}}\mathit{\epsilon }\left(\mathrm{‰}\right)$	Ratio $n/s$
$s_{\mathrm{OH}}$	0.820	$-\mathrm{231}$	$-\mathrm{5.4}$	246.0 $/$ 246.0
$s_{\mathrm{strat}}$	0.092	$-\mathrm{170}$	$-\mathrm{13.1}$	27.5 $/$ 27.5
$s_{\mathrm{soil}}$	0.047	$-\mathrm{80}$	$-\mathrm{22.0}$	19.6 $/$ 8.4
$s_{\mathrm{Cl}}$	0.042	$-\mathrm{459}$	$-\mathrm{58.0}$	9.6 $/$ 15.4
Numbers used for inversion equations
	${\mathit{\tau }}_x$(yr)	${}^{\mathrm{2}}{\mathit{\epsilon }}_x\left(\mathrm{‰}\right)$	${}^{\mathrm{13}}{\mathit{\epsilon }}_x\left(\mathrm{‰}\right)$	$r_x$
North box	7.60	$-\mathrm{222.9}$	$-\mathrm{7.57}$	0.456
South box	9.22	$-\mathrm{232.9}$	$-\mathrm{7.80}$	0.544

Note that the sizes of the two well-mixed boxes are not equal; the sizes are defined by the global mean annual latitude of the inter-tropical convergence zone (ITCZ) ${\mathit{\phi }}_{\mathrm{ITCZ}}=\mathrm{5}$${}^{\circ }$ N (Marshall et al., 2014), which acts as a barrier to the hemispheric exchange of air masses. The inequality of the box sizes is reflected in Eqs. (1) and (2) by $${\displaystyle }{r}_x={\displaystyle \frac{V_x}{V_{\mathrm{atm}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mp \sin \left({\mathit{\phi }}_{\mathrm{ITCZ}}\right)\right),$$ the ratio of the box volume relative to the total atmospheric volume ($r_{\mathrm{n}}=\mathrm{0.456}$, $r_{\mathrm{s}}=\mathrm{0.544}$ for ${\mathit{\phi }}_{\mathrm{ITCZ}}=\mathrm{5}$${}^{\circ }$ N). The global mean lifetime of ${\mathrm{CH}}_{\mathrm{4}}$ quantifies the total methane sink and is set to $\mathit{\tau }=\mathrm{8.4}$ yr (as in the four-box model used in Fischer et al., 2008, and in Bock et al., 2017). Again, using other estimates of the atmospheric lifetime (Prather et al., 2012; Naik et al., 2013) changes the absolute emission estimates but does not affect our conclusions on relative emission changes. According to the relative importance given by Kirschke et al. (2013) (average of bottom-up estimates of all three time intervals shown in their study) the sink is proportionally divided into the four different sink processes ($s_{\mathrm{OH}}$, $s_{\mathrm{strat}}$, $s_{\mathrm{soil}}$, $s_{\mathrm{Cl}}$) shown in Table 5. While $s_{\mathrm{OH}}$ and $s_{\mathrm{strat}}$ are thought to be equally distributed north and south of the ITCZ and therefore are the same in the two boxes, $s_{\mathrm{soil}}$ and $s_{\mathrm{Cl}}$ are partitioned according to the area ratios of land and open water, respectively, derived from a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ resolution land-cover database (Channan et al., 2014). The lifetime of ${\mathrm{CH}}_{\mathrm{4}}$ in each box, $${\displaystyle }{\mathit{\tau }}_x=\mathit{\tau }\cdot {\displaystyle \frac{\sum s_{\mathrm{tot}}}{\sum s_x}}\cdot r_x,$$ is determined by the north–south distribution of the sinks and the relative box sizes. The chosen values for the atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ lifetime and the strengths and the distribution of the sinks represent a best-guess based on the cited literature. Hence, other studies end up with slightly different values. Note also that we assume the model parameters to remain temporally unchanged over the Holocene.

The different sink processes lead to strong fractionations in both isotopes (Cantrell et al., 1990; Quay et al., 1999; Houweling et al., 2000; Snover and Quay, 2000, and references therein; Saueressig et al., 2001; Feilberg et al., 2005) and yield the hemispheric isotope fractionation factors ${}^i{\mathit{\alpha }}_x=\mathrm{1}{+}^i{\mathit{\epsilon }}_x$, with the fractionations ${}^i{\mathit{\epsilon }}_x$ shown in Table 5. As in our model with only two tropospheric boxes $s_{\mathrm{strat}}$ is applied to the total tropospheric ${\mathrm{CH}}_{\mathrm{4}}$ inventory instead of the isotopically more depleted stratospheric ${\mathrm{CH}}_{\mathrm{4}}$, the literature values for the stratospheric sink fractionation have to be adjusted to ${}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{strat}}=-\mathrm{170}$ ‰ and ${}^{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{strat}}=-\mathrm{13.1}$ ‰. These values have been empirically derived to match the tropospheric isotopic signature of the forward four-box model described by Fischer et al. (2008), which contains two stratospheric boxes. Using the identical ${\mathrm{CH}}_{\mathrm{4}}$ emission set-up in both models, the north–south air mass exchange of the two-box model has been increased to equal the tropospheric [${\mathrm{CH}}_{\mathrm{4}}$] values (IPD) of those of the four-box model. In a second step, the fractionations of the stratospheric sink ${}^i{\mathit{\epsilon }}_{\mathrm{strat}}$ in the two-box model have been adjusted to match the tropospheric $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ values of the two models.

The sink fractionation values used for the four-box model analysis by Bock et al. (2017) are very similar to the values we used in this study. Note that the numbers for the sink fractionation caused by stratospheric loss and by soils were accidentally swapped in the printed version of Table S3 by Bock et al. (2017) for both $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ and $\mathit{\delta }$D. However, their model runs and therefore the results are based on the correct values.

As the individual data points from ice core [${\mathrm{CH}}_{\mathrm{4}}$] and isotope measurements are subject to uncertainty in both time and magnitude even after our homogenisation procedure, we only use the millennial variations in ${\mathrm{CH}}_{\mathrm{4}}$, $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$, and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ in our box model inversion. To this end we calculated low-pass-filtered versions of the measured data using the Monte Carlo spline approximations described in Sect. 2.2, which allow us to quantify the mean and the uncertainty of the millennial variations. The uncertainty of the splines of the tropospheric values are introduced into the emission values using Gaussian error propagation. Due to the mentioned limitations we refrain from analysing shorter-scale variability.

SF${}_{\mathrm{6}}$ calibration of the two box model

The same approach used to set up the ${\mathrm{CH}}_{\mathrm{4}}$ inversion described in the previous section (Eqs. 1 and 2) allows us to formulate the general equations for the change of the concentration of any trace gas species in the two hemisphere boxes: $${\displaystyle \frac{\mathrm{d}{c}_x}{\mathrm{d}t}}\cdot r_x\cdot m^{\ast }=E_x-S_x-{\displaystyle \frac{m^{\ast }}{\mathrm{2}\cdot \mathit{\theta }}}\left(c_x-c_y\right),$$ with $\frac{\mathrm{d}{c}_x}{\mathrm{d}t}$ representing the rate of change of the hemispheric concentration, $m^{\ast }$ the translation factor linking a concentration to an atmospheric inventory, $E_x$ and $S_x$ the emissions and the sinks in box $x$, and $r_x$ the ratio of the box volumes depending on the mean annual latitude of the ITCZ ${\mathit{\phi }}_{\mathrm{ITCZ}}$ as described above, where $x$ and $y$ stand for north and south (or vice versa).

SF${}_{\mathrm{6}}$ is a strong greenhouse gas with a lifetime of more than 1000 years, is solely of anthropogenic origin, and is mainly emitted in the Northern Hemisphere ($r_{E{}_{\mathrm{n}}}=\mathrm{95}$ %, $r_{E{}_{\mathrm{s}}}=\mathrm{5}$ %) (Geller et al., 1997, p. 6; Kovács et al., 2017; Levin et al., 2010). NOAA/ESRL (2017) measurements taken at Alert Station (Canada), Summit (Greenland), Palmer Station (Antarctica), and South Pole (Antarctica) record the SF${}_{\mathrm{6}}$ evolution in both polar regions over the last decades. Since the target in this study is the long-time trend and not seasonal variations, the NOAA/ESRL data are smoothed using a smoothing spline with a cut-off period of 2 years. For every year in the recorded period (1996–2008 CE) the annual SF${}_{\mathrm{6}}$ increase ($\frac{\mathrm{d}{c}_x}{\mathrm{d}t}$) and the mean concentration ($c_x$ and $c_y$) are calculated. The records of the global emissions ($E_{\mathrm{tot}}$), the global atmospheric mean SF${}_{\mathrm{6}}$ concentration ($c_{\mathrm{tot}}$), and the total atmospheric inventory ($m_{\mathrm{tot}}$) published by Levin et al. (2010) are used to calculate the hemispheric emissions $E_x=r_{E_x}\cdot E_{\mathrm{tot}}$ and the translation factor $m^{\ast }=\mathrm{mean}\left(\frac{m_{\mathrm{tot}}}{c_{\mathrm{tot}}}\right)=\mathrm{25.7}\frac{Gg{\mathrm{SF}}_{\mathrm{6}}}{\mathrm{ppt}}$. With the approximation of a negligible sink term (since lifetime exceeds the timescale of interest by 2–3 orders of magnitude), Eq. (1) can be solved for the interpolar exchange time: $${\displaystyle }\mathit{\theta }={\displaystyle \frac{m^{\ast }\left(c_x-c_y\right)}{\mathrm{2}\left(\frac{\mathrm{d}{c}_x}{\mathrm{d}t}\cdot r_x\cdot m^{\ast }-r_{E_x}\cdot E_{\mathrm{tot}}\right)}},$$ where $c_x$ and $c_y$ are the SF${}_{\mathrm{6}}$ concentrations at high latitudes. Note that with this calibration we effectively define the exchange time between the northern and southern polar regions, which is somewhat longer than the mean interhemispheric exchange time. As we use $\mathit{\theta }$ later for our inversion of polar ${\mathrm{CH}}_{\mathrm{4}}$ concentrations measured in the ice cores, we effectively correct for the difference in polar and mean hemispheric ${\mathrm{CH}}_{\mathrm{4}}$ concentration using this SF${}_{\mathrm{6}}$ calibration. We use the term interpolar exchange time throughout the paper. Using Eq. (2) the interpolar exchange time can be calculated for each year. With ${\mathit{\phi }}_{\mathrm{ITCZ}}=\mathrm{5.0}$${}^{\circ }$ N we get the average value of $\mathit{\theta }=\mathrm{1.56}\pm \mathrm{0.17}$ yr.

Another valid strategy to cope with the ice core data representing polar tropospheric values would be to correct the data to represent mean hemispheric (${\mathit{\phi }}_{\mathrm{ITCZ}}=\mathrm{5.0}$${}^{\circ }$ N) values (similarly to how it has been carried out by Etheridge et al. (1998) and by Schaefer et al. (2016) to derive a global average signal from ice core ${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathit{\delta }}^{\mathrm{13}}$C-${\mathrm{CH}}_{\mathrm{4}}$ data) based on the knowledge about the modern ${\mathrm{CH}}_{\mathrm{4}}$ concentration distribution (e.g. as provided by Dlugokencky et al., 2017). However, using the corrected data for the box model inversion also requires the corresponding calculation of $\mathit{\theta }$ (to represent the hemispheric mixing time) using mean hemispheric SF${}_{\mathrm{6}}$ concentrations. The lower spatial resolution of the SF${}_{\mathrm{6}}$ concentration data (especially in the low latitudes) biases the mean hemispheric SF${}_{\mathrm{6}}$ values and thus does not allow a satisfying calculation of a hemispheric exchange time. Therefore, we decided not to follow this approach and use polar values in combination with the polar exchange time (as explained above) for our inversion study.

Artificial time series and inversion results. Artificial tropospheric [${\mathrm{CH}}_{\mathrm{4}}$] (a) and stable isotope time series (b, c) used as input for the box model inversion. The calculated mean hemispheric ${\mathrm{CH}}_{\mathrm{4}}$ emission strengths (d) and their isotopic signatures (e, f) show the dependence on the different input parameters. Three experiments are performed: all input values constant (solid lines), all constant but the northern tropospheric ${\mathrm{CH}}_{\mathrm{4}}$ isotope values varied (dashed lines), and all constant but [${\mathrm{CH}}_{\mathrm{4}}$] varied (dotted lines). To be consistent with the wording in the text, time runs from right to left. Note the inverted $y$ axes in panels (b), (c), (e), and (f).

[Figure omitted. See PDF]

Inversion of artificial time series

The interpretation of the output of our inversion model is challenging and not always intuitive. To get a better understanding of the model output, some simplified experiments using artificial time series have been performed and are shown in this study for didactic reasons. In a first approach all six parameters are held constant over time. In the panels on the left in Fig. 2 the tropospheric [${\mathrm{CH}}_{\mathrm{4}}$] (a) and the isotopic signatures $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ (b) and $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ (c) are shown for the northern box in red and for the southern box in blue. The constant values (solid lines) represent the levels of the splines of our measurements at 2 ka BP and serve as input to the inversion model. The panels on the right show the inversion results, i.e. the integrated emissions (d) and their mean isotopic signatures (e, f) that are emitted into the two hemisphere boxes according to our model. As the input parameters are constant, the emission signal must be constant in all parameters as well. However, although the Northern Hemisphere concentrations are only 6 % larger than the southern ones, a large difference in the emission strengths ($E_{\mathrm{n}}\approx \mathrm{2.2}\cdot E_{\mathrm{s}}$) is needed to maintain the IPD in the [${\mathrm{CH}}_{\mathrm{4}}$], compensating for the air mass exchange between the atmospheres. Also note that due to the sink fractionation the isotopic signature of atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ is strongly enriched in the heavy isotopes compared to the emissions (see Table 5).

For the second experiment (dashed lines in Fig. 2), the IPDs of the tropospheric isotope signatures (IPD${}_{\mathit{\delta }D}$ and IPD${}_{\mathit{\delta }\mathrm{13}\mathrm{C}}$) were changed, while the concentrations remained unchanged. By altering the northern stable isotope signal the isotopic IPDs are changed within the range observed in the splined ice core records (dashed lines in Fig. 2b and c). This leads to changes in the calculated isotope signatures of the emissions, whereas the emission strengths are unaffected. As shown (dashed lines) in Fig. 2e and f, the isotopic signature for both the Northern Hemisphere and the Southern Hemisphere emissions needs to be altered to yield the given atmospheric values of constant [${\mathrm{CH}}_{\mathrm{4}}$] and isotopic signatures changing only in one hemisphere. To achieve the declining isotope signal in the northern troposphere, the isotope signatures of emissions from the Northern Hemisphere have to become isotopically more depleted in heavy isotopes (later on conveniently referred to as isotopically “lighter”, equivalent to lower $\mathit{\delta }$ values) over time. At the same time the isotopically lighter ${\mathrm{CH}}_{\mathrm{4}}$ mixed from the northern into the southern box requires the southern box emissions to become more enriched in the heavy isotope (isotopically “heavier”, equivalent to higher $\mathit{\delta }$ values) over time. Since the southern source is much weaker than the northern source (more northern-sourced ${\mathrm{CH}}_{\mathrm{4}}$ is mixed from north to south than southern-sourced ${\mathrm{CH}}_{\mathrm{4}}$ is mixed from south to north), the change in the isotopic signatures of the emissions turns out to be even larger in the Southern Hemisphere.

The dotted line in panel (a) in Fig. 2 shows a third scenario, in which all tropospheric parameters are kept constant except for the [${\mathrm{CH}}_{\mathrm{4}}$] in the Northern Hemisphere. Similar to the former experiment, the [${\mathrm{CH}}_{\mathrm{4}}$] in the northern box is varied between values representing the minimal and the maximal measured IPD${}_{{\mathrm{CH}}_{\mathrm{4}}}$. The dotted curves in the right panels (d–f) in Fig. 2 indicate that this change has an influence on all six emission parameters of our model. The emission strength in the northern box increases to produce the rising [${\mathrm{CH}}_{\mathrm{4}}$]. Declining emission strength in the Southern Hemisphere compensates for the increasing ${\mathrm{CH}}_{\mathrm{4}}$ north-to-south flux due to the rising IPD${}_{{\mathrm{CH}}_{\mathrm{4}}}$. The ${\mathrm{CH}}_{\mathrm{4}}$ that enters a box by atmospheric mixing is already fractionated by the sinks and isotopically heavier than the newly emitted ${\mathrm{CH}}_{\mathrm{4}}$. If the relative amount of ${\mathrm{CH}}_{\mathrm{4}}$ input by mixing is decreased for the northern (increased for the southern) box, the emission signatures become less (more) depleted in heavy isotopes to maintain the tropospheric ${\mathrm{CH}}_{\mathrm{4}}$ isotopic signatures.

Data and inversion. The panels on the left (a–c) show the concentration and the two stable isotope signatures of ${\mathrm{CH}}_{\mathrm{4}}$ measured in ice cores representing the tropospheric values over the Holocene. The concentration data are complemented with data from other studies (see Table 1). All three parameters were measured on ice cores from Greenland (red) and from Antarctica (blue and cyan) representing the northern and the southern polar troposphere, respectively. A spline function with a cut-off period of 3000 years has been used to calculate the smoothed evolution of the tropospheric signals represented by the lines and the $\mathrm{1}\mathit{\sigma }$ and $\mathrm{2}\mathit{\sigma }$ error bands. The panels (g), (h), and (i) show the IPDs of all three parameters also with their $\mathrm{1}\mathit{\sigma }$ and $\mathrm{2}\mathit{\sigma }$ error bands. A box model inversion is used to calculate the hemispheric ${\mathrm{CH}}_{\mathrm{4}}$ emissions shown in the panels on the right. In panel (d) the emission strengths are shown for both hemisphere boxes together with the total ${\mathrm{CH}}_{\mathrm{4}}$ emissions. Panels (e) and (f) show the mean isotopic signatures of the emissions for both stable isotopes of ${\mathrm{CH}}_{\mathrm{4}}$. All the error bands represent $\mathrm{1}\mathit{\sigma }$ and $\mathrm{2}\mathit{\sigma }$ uncertainties. Note the inverted $y$ axes in panels (b), (c), (e), (f), (h), and (i). Red refers to Northern Hemisphere records and blue to Southern Hemisphere records throughout the paper.

[Figure omitted. See PDF]

Since the variations for both experiments were performed based on a reasonable range (observed changes in IPDs), the span for the calculated values of the emission parameters also provides a range for the changes that are needed to alter the measured tropospheric IPDs in the Holocene. Note, however, that the variation scenarios (e.g. a full decoupling of atmospheric concentration and the isotopic composition of ${\mathrm{CH}}_{\mathrm{4}}$) do not represent changes we expect to observe in nature and are presented for illustrative reasons only.

Results

The homogenised tropospheric data (this study complemented with published [${\mathrm{CH}}_{\mathrm{4}}$] data from Chappellaz et al., 1993, 1997; Blunier et al., 1995; Brook et al., 1996; Brook, 2009; Flückiger et al., 2002; Schilt et al., 2010; Mitchell et al., 2013; see Sect. 2.3, Fig. 1, and Table 1) are shown in Fig. 3 in the left-hand panels. The smoothing splines, which are all calculated with the same cut-off period of 3000 years, are shown including their $\mathrm{1}\mathit{\sigma }$ and $\mathrm{2}\mathit{\sigma }$ uncertainty bands. In Fig.3a the composite of the time-synchronised and offset-corrected ${\mathrm{CH}}_{\mathrm{4}}$ data is shown over the Holocene. The diamond-shaped symbols represent the NGRIP and TALDICE data of this study measured during the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ measurement campaign. Overall, our data compilation after the above-mentioned corrections supports previous trends: after high concentrations in the Preboreal (ca. 11.5–10 ka BP) the data show a decline of the atmospheric [${\mathrm{CH}}_{\mathrm{4}}$] over the first half of the Holocene followed by a reversal of this trend and an accelerating increase from the mid-Holocene to the late Holocene. The preindustrial level is comparable to the [${\mathrm{CH}}_{\mathrm{4}}$] during the Preboreal warm period. The splines emphasise this long-term evolution nicely but, deliberately, do not reproduce the fine structure in the data (see for example Mitchell et al., 2013). The difference between the [${\mathrm{CH}}_{\mathrm{4}}$] values in the Northern Hemisphere and Southern Hemisphere (the IPD${}_{{\mathrm{CH}}_{\mathrm{4}}}$) is shown in Fig. 3g. The concentration measured in ice cores from Greenland (red) is on average 45 ppb higher than the values derived from Antarctic ice cores (blue), which agrees with the result of Chappellaz et al. (1997).

In Fig. 3b the $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ values measured in ice from the NGRIP and the EDML ice cores are shown. As additionally illustrated in Fig. 3h we find a pronounced difference of $-\mathrm{16.3}$ ‰ between northern and southern $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ values, which is rather stable over the whole period of investigation. Limited by the coarser resolution of the EDML $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ data, the smoothed curves and the error bands illustrate the probable evolution of the true tropospheric signals only on millennial timescales and longer. Compared to the large [${\mathrm{CH}}_{\mathrm{4}}$] changes, the $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ signal over the Holocene is relatively small and similar in the two data sets. The values become isotopically lighter in the first half of the Holocene and heavier again in the second half. Note that our NGRIP $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ data differ significantly from the data set from the GISP2 ice core published by Sowers (2010) (see Fig. 4). First, there is a difference of the mean level between the two data sets of 10 to 15 ‰, which can be attributed to an inter-laboratory difference (Bock et al., 2014; Umezawa et al., 2018). Second, the GISP2 data suggested a shift of about 20 ‰ towards heavier values from 5 to 1 ka BP, which is not confirmed by our data. This disagreement, however, may be largely attributed to the much larger scattering of the GISP2 data set, which did not yet allow the quantification of an unambiguous trend in $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$.

Comparison of stable ${\mathrm{CH}}_{\mathrm{4}}$ isotope data. This study's stable ${\mathrm{CH}}_{\mathrm{4}}$ isotope values of the NGRIP ice cores are shown together with the data published by Sowers (2010) using the GISP2 ice core. The $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ data sets (a) differ in the absolute level and the variance, whereas the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ data (b) are in good agreement. Note that all data are corrected for gravitational fractionation in the firn column. All data of this study are free of krypton interference, whereas the GISP2 $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ record has been subsequently corrected for it (Schmitt et al., 2013). Note the inverted $y$ axes.

[Figure omitted. See PDF]

In Fig. 3c the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ measured on ice from the NGRIP and TALDICE ice cores are shown. Both data sets show a gradual trend towards lighter values over the first half of the Holocene and almost constant values over the second half. The point of inflection is well defined at about 5 ka BP. As shown in Fig. 3f, at about the same time the IPD${}_{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}$ changes quickly from $-$0.35 ‰ to $-\mathrm{0.63}$ ‰, which is rather constant before and after the point of inflection. Our NGRIP $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ record is well in line with the GISP2 data by Sowers (2010) (see Fig. 4), which have been subsequently corrected for the analytical krypton interference (Schmitt et al., 2013).

Our inversion model calculates the emissions that led to the measured tropospheric signals shown in the panels on the left in Fig. 3. In Fig. 3d the strength of the integrated emissions (in Tg ${\mathrm{CH}}_{\mathrm{4}}$ per year) in each of the two troposphere boxes are shown. Additionally, the total ${\mathrm{CH}}_{\mathrm{4}}$ global emissions (sum of north and south) are shown in grey. Note that the uncertainty ($\mathrm{1}\mathit{\sigma }$ and $\mathrm{2}\mathit{\sigma }$) in emission fluxes is very small and much smaller than the long-term changes in the emission fluxes over the Holocene indicated by the spline average. This uncertainty reflects only the error in emission fluxes for changes on millennial timescales that are resolved in our spline. As indicated in high-resolution data (for example Mitchell et al., 2013) the centennial variability in atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ is high and the uncertainties of reconstructed emission fluxes on such shorter timescales would be significantly higher as well.

The following observations can be made from the inversion results: the ${\mathrm{CH}}_{\mathrm{4}}$ emissions into the northern box are more than 2 times higher than the emissions into the southern box. After a decrease of roughly 30 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$ in the early Holocene the northern emissions remain relatively constant between 8 and 2 ka BP before increasing by 15 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$. In the Southern Hemisphere the emissions suggest a slight increase in the first 2000 years of the record followed by a significant 22 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$ decrease until 4–5 ka BP after which emissions increase by 25 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$ in the second half of the Holocene.

The weighted averaged $\mathit{\delta }$D signature of all the integrated ${\mathrm{CH}}_{\mathrm{4}}$ emissions is shown in Fig. 3e for each hemisphere. With values of about $-\mathrm{295}$ ‰ the northern emission $\mathit{\delta }$D signature is on average 30 ‰ lighter than the Southern Hemisphere emissions. Before 7.5 ka BP and around 4.5 ka BP the $\mathit{\delta }$D-$E_x$ signal shows long-term variations around the Holocene mean that are in antiphase between the Northern Hemisphere and Southern Hemisphere, whereas the global mean $\mathit{\delta }$D-E (grey line) remains essentially constant. Based on these two observations and the knowledge from our experiments with artificial time series (see Sect. 2.4.2) showing that a change in the isotopic signature in one hemisphere has a large influence of different sign on the emission isotope signal in the other hemisphere (Fig. 2), we attribute the $\mathit{\delta }$D-$E_x$ fluctuations to the measurement uncertainty of the individual data points of the low-resolution $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ EDML record which affect the spline reconstruction. The absence of a long-term trend in $\mathit{\delta }$D-$E_x$ over the Holocene despite the significant changes in the atmospheric [${\mathrm{CH}}_{\mathrm{4}}$] is remarkable and a strong constraint on the average Holocene ${\mathrm{CH}}_{\mathrm{4}}$ budget. However, the large uncertainty in $\mathit{\delta }$D-$E_x$ does not allow us to make robust conclusions about millennial variations in the hydrogen isotopic signature of ${\mathrm{CH}}_{\mathrm{4}}$ emissions.

The $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ values of the emissions shown in Fig. 3f indicate a significant shift of the northern emissions to heavier values over the Holocene, which is especially pronounced in the time interval before 4 ka BP. In the south, $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-$E_{\mathrm{s}}$ also shows a slight negative trend over the first half of the Holocene, which, however, is still well within the 1$\mathit{\sigma }$ inversion uncertainty. A stronger positive trend, however, occurs during the second half of the Holocene, which, however, is still within the $\mathrm{2}\mathit{\sigma }$ uncertainty of the reconstruction. The global mean of the isotopic emission signatures reflects the coherent trends in the measured atmospheric values, with a strong negative trend of $-\mathrm{1.8}$ ‰ over the first 5000 years of the Holocene.

Sensitivity studies

The size ratio of the two boxes of the box model inversion is defined by the annual mean latitude of the ITCZ ${\mathit{\phi }}_{\mathrm{ITCZ}}$. In fact, there is evidence from different proxies of a southward migration of the ITCZ during the Holocene (Haug et al., 2001; Zhao and Harrison, 2012; McGee et al., 2014). Since it is not possible to deduce robust numbers for the position and the movement of the global mean ITCZ from local studies, ${\mathit{\phi }}_{\mathrm{ITCZ}}$ has been kept constant over time at a reasonable Holocene value in our standard inversion (described in Sect. 2.4). However, sensitivity runs have been performed to quantify the impact of this assumption on the conclusions drawn from our inversion results.

Sensitivity study “SF${}_{\mathrm{6}}$ calibration”. In panels (a–c) the data splines with the $\mathrm{1}\mathit{\sigma }$ uncertainty bands of the measured tropospheric values are shown. In panels (d–f) the calculated strengths and the isotopic signature of the emissions are shown for the three different ${\mathit{\phi }}_{\mathrm{ITCZ}}$ and $\mathit{\theta }$ set-ups to investigate the sensitivity of the inversion to the SF${}_{\mathrm{6}}$ calibration. Note the inverted $y$ axes in panels (b), (c), (e), and (f).

[Figure omitted. See PDF]

As described in Sect. 2.4.1, SF${}_{\mathrm{6}}$ calibration has been used to determine the interpolar exchange time $\mathit{\theta }$ associated with different values of ${\mathit{\phi }}_{\mathrm{ITCZ}}$. For the first sensitivity study these two parameters have always been fitted in parallel to ensure that the different model configurations all fit the SF${}_{\mathrm{6}}$ constraint in an optimum way. In Fig. 5 the splines of the Holocene data (left panels) and the corresponding emission values derived by the inversion (right panels) are shown. In the sensitivity runs, the inversion equations have been solved using other (${\mathit{\phi }}_{\mathrm{ITCZ}}$, $\mathit{\epsilon })$ couples with ${\mathit{\phi }}_{\mathrm{ITCZ}}$ moving from 0.0 to 10.0${}^{\circ }$ N. In all six emission parameters, these variations only lead to small changes relative to the total signal variability. For the isotopic signatures of the emissions, all curves lie well within the $\mathrm{1}\mathit{\sigma }$ error of the 5.0${}^{\circ }$ N results.

Sensitivity study “ITCZ position”. In panels (a–c) the data splines with the $\mathrm{1}\mathit{\sigma }$ uncertainty error bands of the measured tropospheric values are shown. In panels (d–f) the calculated strengths and the isotopic signature of the emissions are shown for the three different ${\mathit{\phi }}_{\mathrm{ITCZ}}$ set-ups and the best-guess value for $\mathit{\theta }$ to investigate the sensitivity of the inversion to a shift of the ITCZ. Note the inverted $y$ axes in panels (b), (c), (e), and (f).

[Figure omitted. See PDF]

A change of the mean annual position of the ITCZ also has an influence on the box sizes. This has a direct impact on the north–south partitioning of the individual ${\mathrm{CH}}_{\mathrm{4}}$ sinks and thus on all six inversion parameters. Therefore, another sensitivity run was carried out, in which only ${\mathit{\phi }}_{\mathrm{ITCZ}}$ was varied and $\mathit{\theta }$ kept constant at the best guess value. As shown in Fig. 6, the sensitivity of the inversion to a change of ${\mathit{\phi }}_{\mathrm{ITCZ}}$ is very low if $\mathit{\theta }$ is kept constant. Note that the southward migration of the ITCZ over the Holocene is thought to be smaller than 1${}^{\circ }$ (McGee et al., 2014). The results of this second set of experiments in which only one parameter was changed (i.e. the constraints of the SF${}_{\mathrm{6}}$ calibration were not fulfilled) also shows that on average variations in ${\mathit{\phi }}_{\mathrm{ITCZ}}$ and $\mathit{\theta }$ contribute a similar amount to the deviation from our standard run (${\mathit{\phi }}_{\mathrm{ITCZ}}=\mathrm{5}$${}^{\circ }$ N) (see Fig. 6).

Discussion

The main goal of this study is to profit from the improved ${\mathrm{CH}}_{\mathrm{4}}$ multi-parameter data set using our inversion model in order to evaluate the different scenarios and to constrain the global ${\mathrm{CH}}_{\mathrm{4}}$ budget over the Holocene. To this end, we consider the early, mid-Holocene, and late Holocene separately.

Early Holocene (11–8 ka BP)

The significant early drop in the northern emission strength together with the slight shift towards isotopically lighter carbon emissions in the Northern Hemisphere from 11 to 8 ka BP is challenging to account for. One potential solution to account for the observed concentration shift would be a general weakening of the boreal ${\mathrm{CH}}_{\mathrm{4}}$ sources as a result of the decreasing northern summer insolation (Berger and Loutre, 1991). However, this should lead to a trend in both isotope emission signatures towards heavier values in the Northern Hemisphere, as in general the high-latitude ${\mathrm{CH}}_{\mathrm{4}}$ emissions from wetlands are isotopically light in both isotopes (Whiticar and Schaefer, 2007; Walter et al., 2008). Accordingly, this conflicts with the stable isotope results from our inversion. For $\mathit{\delta }$D-$E_{\mathrm{n}}$ the inversion might be affected by one single measurement point (EDML $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$) at 8.19 ka BP, forcing the emission signature curves of the two hemispheres apart and indicating that $\mathit{\delta }$D-$E_{\mathrm{n}}$ becomes lighter. Nevertheless, even if this single data point is removed, the inversion signal would be constant over this period but would not show a trend towards heavier values as expected from a decline in boreal emissions.

Based on ${}^{\mathrm{14}}$C dating, Walter et al. (2007) estimated the ${\mathrm{CH}}_{\mathrm{4}}$ emissions from thermokarst lakes to decrease from 26 to 5 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$ in the period from 11 to 8 ka BP in line with the decreasing trend in Northern Hemisphere emissions in our inversion. To match the result of the isotope inversion, another isotopically heavy northern source (e.g. biomass burning) has to be reduced simultaneously to compensate for the lack of isotopically light thermokarst ${\mathrm{CH}}_{\mathrm{4}}$ emissions. However, charcoal records from sediments do not indicate a decline of wild fires in the Northern Hemisphere but rather point to the opposite during the Holocene (Daniau et al., 2012).

Part of the answer of this puzzle may lie in the combined temporal evolution of the strengths and isotopic signature of the boreal wetlands. Due to land ice and permafrost retreat in the early Holocene, minerotrophic fens, characterised by relatively high ${\mathrm{CH}}_{\mathrm{4}}$ emissions which are not as strongly depleted in ${}^{\mathrm{13}}$C, turn into ombrotrophic bogs over time with more ${}^{\mathrm{13}}$C-depleted but also much weaker ${\mathrm{CH}}_{\mathrm{4}}$ emissions (Ding et al., 2005; Sowers, 2010; Hornibrook, 2013; Yu et al., 2013). This reduction of ${\mathrm{CH}}_{\mathrm{4}}$ emissions in the course of conversion of fens to bogs requires a very depleted ${}^{\mathrm{13}}$C source signature of the bogs to quantitatively explain our inversion results and may therefore not suffice to explain the ${\mathit{\delta }}^{\mathrm{13}}$C-$E_{\mathrm{n}}$ shift towards lighter values in the early Holocene. Furthermore, an enrichment trend in $\mathit{\delta }$D-$E_{\mathrm{n}}$ (of about 5 ‰) would be expected if the $\mathit{\delta }$D-depleted high-latitude emissions became relatively less important, while – if at all – we see a 20 ‰ depletion in $\mathit{\delta }$D-$E_n$ from 10 to 8 ka BP. Also the isotopic signature of the water used in methanogenesis has a direct imprint on the hydrogen isotopic signature of the emitted ${\mathrm{CH}}_{\mathrm{4}}$ (Whiticar and Schaefer, 2007). Speleothem $\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}$ records confirm an isotopic shift towards lighter values in the meteoric water in the summer monsoon regions of the Northern Hemisphere from 11 to 8 ka BP (Wang et al., 2014) and an inverse relationship in the Southern Hemisphere, in line with our inversion. However, the speleothem records also suggest opposite trends over the rest of the Holocene, which are not seen in our inversion results. In addition, isotopically light meltwater from the retreating Northern Hemisphere ice sheets (for example in proglacial lakes) may contribute to a depletion in the hydrogen isotopic signature of boreal ${\mathrm{CH}}_{\mathrm{4}}$ source in the early Holocene.

Mid-Holocene (8–5.5 ka BP)

Our inversion results show that the decrease in atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ concentrations over this interval is mainly caused by declining Southern Hemisphere emissions. This decrease in the Southern Hemisphere ${\mathrm{CH}}_{\mathrm{4}}$ emission from 8 to 5.5 ka BP – during the time of the strongest increase in austral summer insolation (Berger and Loutre, 1991) – appears counterintuitive but is in line with the model results by Singarayer et al. (2011), in which this is caused by changes in tropical rainfall location. This is also the period when the ITCZ is thought to have started moving southward as a response to the change in the orbital forcing (Haug et al., 2001; Zhao and Harrison, 2012; McGee et al., 2014), having an influence on the hemispheric ${\mathrm{CH}}_{\mathrm{4}}$ source distributions. ${\mathrm{CH}}_{\mathrm{4}}$ emissions from regions that were attributed to the Southern Hemisphere box in the early Holocene may now contribute ${\mathrm{CH}}_{\mathrm{4}}$ to the Northern Hemisphere due to the shift of the ITCZ. Note that such a change has an influence not only on the spatial emission distribution but also on other parameters (such as $\mathit{\theta })$, which are relevant for the inversion model, in which the box sizes were kept constant over time. Our sensitivity studies (see Sect. 3.1) show, however, that a movement of the ITCZ in the model has only a small effect on the inversion results.

As discussed above, the $\mathit{\delta }$D signal of the ${\mathrm{CH}}_{\mathrm{4}}$ emissions is mostly constant in both hemispheres. This is expected if one hemisphere box gains emissions from tropical wetlands, characterised by stable isotope signatures close to the integrated source mix, at the cost of the other box. The rather constant charcoal wildfire reconstructions for the Northern Hemisphere in this time window (Daniau et al., 2012) do not support a decline in biomass burning ${\mathrm{CH}}_{\mathrm{4}}$ emissions being responsible for the 1.5 ‰ depletion in $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-$E_{\mathrm{n}}$ between 7 and 5.5 ka BP. The $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ shift of the emissions could be related to a shift in the ratio of C3 to C4 plants as suggested by Sowers (2010), while a further evolvement of northern fens to bogs is incompatible with the relatively constant emissions in the Northern Hemisphere. Under Holocene climate conditions, with high availability of humidity and ${\mathrm{CO}}_{\mathrm{2}}$, C${}_{\mathrm{3}}$ plants successfully compete against the relatively ${}^{\mathrm{13}}$C-enriched C${}_{\mathrm{4}}$ plants and thereby alter the carbon isotopic signature of the precursor material of methanogenesis (Farquhar et al., 1989). This process was also suggested by Möller et al. (2013) to explain the close correlation of the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$-${\mathrm{CH}}_{\mathrm{4}}$ signal with atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and the decoupling from the [${\mathrm{CH}}_{\mathrm{4}}$] in the glacial period.

Late Holocene (5.5–1 ka BP)

The early anthropogenic influence hypothesis (Ruddiman and Thomson, 2001) calls upon increasing emissions from rice paddies mainly in eastern Asia. In our box model such emissions would be located in the Northern Hemisphere box and therefore lead to an increase in northern emissions after 5 ka BP. However, the box model inversion shows no significant increase in $E_{\mathrm{n}}$ before 2 ka BP. In contrast, at 5.5 ka BP, at the time when all six presented tropospheric records imply a change in the global ${\mathrm{CH}}_{\mathrm{4}}$ cycle, the inversion indicates an amplification of Southern Hemisphere ${\mathrm{CH}}_{\mathrm{4}}$ emissions. This clearly supports the idea of a natural Holocene wetland ${\mathrm{CH}}_{\mathrm{4}}$ emission increase as suggested by Singarayer et al. (2011). In contrast, a southern anthropogenic source can be excluded since the Southern Hemisphere was barely populated in the mid-Holocene (Kaplan et al., 2011). The observed shift in $\mathit{\delta }$D-$E_{\mathrm{s}}$ and especially ${\mathit{\delta }}^{\mathrm{13}}$C-$E_{\mathrm{s}}$ towards heavier values in this interval cannot be explained simply by increased tropical wetland emissions in the Amazon region but also requires an enhancement of a Southern Hemisphere ${\mathrm{CH}}_{\mathrm{4}}$ source strongly enriched in ${}^{\mathrm{13}}$C. Using a simple mass balance approach ($E_{\mathrm{total}}{\cdot }^i{\mathit{\delta }}_{\mathrm{total}}=E_{\mathrm{initial}}{\cdot }^i{\mathit{\delta }}_{\mathrm{initial}}+E_{\mathrm{additional}}{\cdot }^i{\mathit{\delta }}_{\mathrm{additional}}$) allows us to determine a $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ signature of the additional ${\mathrm{CH}}_{\mathrm{4}}$ emissions of 19 Tg ${\mathrm{CH}}_{\mathrm{4}}$ yr${}^{-\mathrm{1}}$ to be about $-\mathrm{45.5}$ ‰, which is isotopically heavier than tropical wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions ($-\mathrm{56.8}$ ‰ according to Whiticar and Schaefer, 2007). Note, however, that this quantitative approach is based on the assumption of unchanged carbon isotopic signatures of the initial ${\mathrm{CH}}_{\mathrm{4}}$ emissions and therefore needs to be taken with care. Increasing emissions from tropical and subtropical wild fires in the Southern Hemisphere, which are also documented in charcoal records (Daniau et al., 2012), could have become increasingly important, e.g. due to enhanced fuel production or higher ${\mathrm{CH}}_{\mathrm{4}}$ emission efficiency (more smoldering fires), both caused by wetter climatic conditions, and would readily explain the joint information derived from inversion results on Southern Hemisphere ${\mathrm{CH}}_{\mathrm{4}}$ emissions and their stable isotopic signatures. Other isotopically heavy sources such as geological emissions and marine ${\mathrm{CH}}_{\mathrm{4}}$ hydrates are unlikely to dramatically gain importance in the mid-Holocene and only in one hemisphere and can therefore be ruled out as important players driving the observed ${\mathrm{CH}}_{\mathrm{4}}$ changes.

Conclusions

Our new records of [${\mathrm{CH}}_{\mathrm{4}}$], $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$, and ${\mathit{\delta }}^{\mathrm{13}}$C-${\mathrm{CH}}_{\mathrm{4}}$ measured on samples from polar ice cores from both hemispheres provide valuable insights into the Holocene ${\mathrm{CH}}_{\mathrm{4}}$ cycle. A significant IPD for $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ and $\mathit{\delta }$D, as already documented for [${\mathrm{CH}}_{\mathrm{4}}$], exists over the entire period of investigation. While the [${\mathrm{CH}}_{\mathrm{4}}$] data confirm the well-known Holocene evolution, the two stable isotope records provide additional insights. This shows the value of our multi-isotope approach in order to avoid drawing wrong conclusions. A two-box model approach has been used to translate the measured signals into emissions. The calculated hemispheric ${\mathrm{CH}}_{\mathrm{4}}$ emission signatures, characterised by different trends in the two hemispheres for ${\mathit{\delta }}^{\mathrm{13}}$C and no statistically significant variation in $\mathit{\delta }$D, show a decoupling of the emission signatures of the two analysed isotopes. Therefore, we attribute the long-term changes of the atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ isotopes over the first half of the Holocene to shifts in the isotopic source signatures of individual ${\mathrm{CH}}_{\mathrm{4}}$ sources rather than to changes in the global ${\mathrm{CH}}_{\mathrm{4}}$ source mix.

In the early Holocene (11–8 ka BP) we associate the decline of Northern Hemisphere emissions and the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ trend towards lower values with the natural evolution of high-latitude wetlands from fens (high ${\mathrm{CH}}_{\mathrm{4}}$ emissions) to bogs (lower but strongly ${}^{\mathrm{13}}$C-depleted ${\mathrm{CH}}_{\mathrm{4}}$ emissions). An alternative explanation calling upon a significant decrease in ${\mathrm{CH}}_{\mathrm{4}}$ emissions from thermokarst lakes would require compensation by other processes to fulfil the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ constraint, e.g. a substantial decline of wild fire activity in the Northern Hemisphere. In the mid-Holocene (8–5.5 ka BP) the inversion shows a decline of the Southern Hemisphere ${\mathrm{CH}}_{\mathrm{4}}$ emissions and an ongoing shift of the $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ emission signatures towards lighter values, which might be related to the southward migration of the ITCZ and the change of the C3-to-C4 plant ratio altering the carbon isotopic composition of the precursor material for ${\mathrm{CH}}_{\mathrm{4}}$ production. During the second half of the Holocene (5.5–1 ka BP), which shows a major increase in ${\mathrm{CH}}_{\mathrm{4}}$ emissions accompanied by a shift towards higher $\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}$ signatures in the Southern Hemisphere, our results favour the ${\mathrm{CH}}_{\mathrm{4}}$ rise to be caused by natural ${\mathrm{CH}}_{\mathrm{4}}$ emissions (both wetlands and wildfires) in the southern tropics rather than by early rice agriculture in east Asia.

The ${\mathrm{CH}}_{\mathrm{4}}$ stable isotope data provide a powerful constraint on the Holocene ${\mathrm{CH}}_{\mathrm{4}}$ system (e.g. to benchmark ${\mathrm{CH}}_{\mathrm{4}}$ emission models); however, the relatively large uncertainty ranges in the calculated emission parameters still limit our ability to draw more robust conclusions on temporal variations. The inversion would benefit from an increased temporal resolution in the Southern Hemisphere $\mathit{\delta }$D-${\mathrm{CH}}_{\mathrm{4}}$ record. However, since the concentration data have a large influence on all six inversion parameters, we see the largest potential to better constrain the Holocene ${\mathrm{CH}}_{\mathrm{4}}$ emissions by actually improving the [${\mathrm{CH}}_{\mathrm{4}}$] records. Better IPD${}_{{\mathrm{CH}}_{\mathrm{4}}}$ data measured on the same analytical system during the same measurement campaign in high temporal resolution, which ensures an accurate synchronisation (such as the data by Mitchell et al. (2013) for the last 2800 years), would bring us a large step forward in further constraining the Holocene ${\mathrm{CH}}_{\mathrm{4}}$ cycle. This is especially true for the Greenland ice cores within the brittle ice zone, e.g. GISP2 records for depths between 650 and 1400 m, equivalent to the time window 2.7–8.1 ka BP when individual samples subject to modern air intrusion may compromise the Greenland ${\mathrm{CH}}_{\mathrm{4}}$ record. The new EGRIP ice core currently drilled in northeast Greenland has much lower accumulation; hence a brittle ice zone located at older ages may help to improve the ${\mathrm{CH}}_{\mathrm{4}}$ record in this time interval in the future. Alternatively, a high-accumulation ice core from Greenland where ice older than 2.7 ka BP is found below the brittle ice zone would circumvent this issue and allow for uncompromised ${\mathrm{CH}}_{\mathrm{4}}$ measurements over the Holocene.

Data presented in this study will be available at https://www.ncdc.noaa.gov/paleo/study/25350 (last access: 13 October 2018) (NOAA, 2018).

The authors declare that they have no conflict of interest.

Acknowledgements

The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme FP7/2007-2013 ERC grant 226172 (ERC Advanced Grant Modern Approaches to Temperature Reconstructions in Polar Ice Cores (MATRICs)) and the Swiss National Science Foundation. This work is a contribution to the NorthGRIP ice core project, which is directed and organised by the Department of Geophysics at the Niels Bohr Institute for Astronomy, Physics and Geophysics, University of Copenhagen. It is supported by funding agencies in Denmark (SNF), Belgium (FNRS-CFB), France (IFRTP and NSU/CNRS), Germany (AWI), Iceland (RannIs), Japan (MEXT), Sweden (SPRS), Switzerland (SNF), and the United States (NSF). Furthermore, this work is a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Science Foundation–European Commission scientific program funded by the European Union and national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland, and the United Kingdom. The main logistic support was provided by Institut Polaire Français Paul-Emile Victor (IPEV) and PNRA (at Dome C) and AWI (at Dronning Maud Land). The TALos Dome Ice CorE (TALDICE) project, a joint European program led by Italy, is funded by national contributions from Italy, France, Germany, Switzerland, and the United Kingdom. The main logistical support was provided by Programma Nazionale di Ricerche in Antartide (PNRA) at Talos Dome. This work is EPICA publication no. 309 and TALDICE publication no. 51. Edited by: Kirsten Thonicke Reviewed by: Hinrich Schaefer and one anonymous referee