Can the carbon isotopic composition of methane be reconstructed from multi-site ﬁrn air measurements?

. Methane is a strong greenhouse gas and large uncertainties exist concerning the future evolution of its atmospheric abundance. Analyzing methane atmospheric mixing and stable isotope ratios in air trapped in polar ice sheets helps in reconstructing the evolution of its sources and sinks in the past. This is important to improve predictions of atmospheric CH 4 mixing ratios in the future under the inﬂuence of a changing climate. The aim of this study is to assess whether past atmospheric δ 13 C(CH 4 ) variations can be reliably reconstructed from ﬁrn air measurements. Isotope reconstructions obtained with a state of the art ﬁrn model from different individual sites show unexpectedly large discrepancies and are mutually inconsistent. We show that small changes in the diffusivity proﬁles at individual sites lead to strong differences in the ﬁrn fractionation, which can explain a large part of these discrepancies. Using slightly modiﬁed diffusivities for some sites, and neglecting samples for which the ﬁrn fractionation signals are strongest, a combined multi-site can be which with are lower than been from Southern Hemisphere (SH) archived air samples and high-accumulation ice core data. We conclude that with the current datasets and understanding of ﬁrn air transport, a high preci-sion reconstruction of δ 13 C of CH 4 from ﬁrn air samples is not possible, because reconstructed atmospheric trends over the last 50 yr of 0.3–1.5 ‰ are of the same magnitude as in-herent uncertainties in the method, which are the ﬁrn fractionation correction (up to ∼ 2 ‰ at individual sites), the Kr isobaric interference (up to ∼ 0.8 ‰, system dependent), inter-laboratory calibration offsets ( ∼ 0.2 ‰) and uncertainties in past CH 4 levels ( ∼ 0.5 ‰).

because the CH 4 mixing ratio depends on the strength of its numerous sources (e.g. wetlands, ruminants, rice paddies, biomass burning, natural gas) and on the rate at which it is removed from the atmosphere by its different sinks (oxidation by tropospheric OH, stratospheric loss and soil uptake).
Air enclosed in polar firn and ice represents an archive of past atmospheric composition (e.g. Craig and Chou, 1982;Schwander et al., 1993;Sowers et al., 2005). The snow accumulates at the surface of the ice sheets and gradually densifies which leads to a slow trapping of the air in bubbles in the interstitial spaces. Firn represents the porous layer from the ice sheet surface down to the bubble close-off depth (40-120 m), where all the air is occluded in the ice, i.e. totally isolated from the atmosphere above. Trace gases moving downward with air in the firn column are affected by several processes: molecular diffusion and gravitational settling (Craig et al., 1988;Schwander et al., 1989;Sowers et al., 1989;Spahni et al., 2003), thermal diffusion and turbulent transport (Severinghaus et al., 2001), convection in the upper firn (Sowers et al., 1992) and advection caused by the firn sinking and the air being trapped in bubbles (Stauffer et al., 1985;Schwander et al., 1993). Firn transport models including some, but not all of these processes (e.g. Schwander et al., 1993;Rommelaere et al., 1997;Buizert et al., 2012;Witrant et al., 2012 and references therein) can be used to reconstruct the evolution of trace gas concentrations in the atmosphere based on depth-dependent physical parameters of the firn (density, porosity, and tortuosity), vertical trace gas profiles in firn and site-dependent surface conditions (temperature, precipitation and pressure).
Previous studies have shown the potential of reconstructing past variations in the CH 4 budget using mixing and stable isotope ratio measurements from firn air (Etheridge et al., 1998;Francey et al., 1999;Bräunlich et al., 2001;Sowers et al., 2005). These studies reported CH 4 stable isotope data of firn air from Antarctica and suggested a rise in the stable carbon isotope of atmospheric CH 4 δ 13 C(CH 4 ) during the 20th century, which was ascribed to an increase of anthropogenic 13 C-enriched (fossil and pyrogenic) CH 4 sources during this period.
Here we combine previously published firn air δ 13 C(CH 4 ) data from 5 sites and new data from 6 sites from both Greenland and Antarctica. Our goal is to reconstruct the δ 13 C(CH 4 ) atmospheric history of both hemispheres using a multi-site inversion. To this end, we first reconstruct the atmospheric isotope histories from measurements at individual sites, using the new LGGE-GIPSA firn model (Rommelaere et al., 1997;Witrant et al., 2012;Wang et al., 2012). Considering the relatively long lifetime of CH 4 compared to the interhemispheric exchange time and the fact that both hemispheres are well mixed on multi-annual timescales, we assume that in each hemisphere all reconstructions from the polar sites should agree. We discuss the differences among the 11 firn records and the possible processes responsible. We then perform a multi-site inversion to reconstruct the δ 13 C(CH 4 ) history from firn air over the last 50 yr. Finally, the consistency of these firn air results with atmospheric data and ice core results as well as the uncertainties related to the multi-site reconstruction are discussed.

Firn air sampling and isotope analysis
Firn air samples discussed in this paper come from 11 boreholes from both Greenland and Antarctica and were extracted between 1993 and 2009 ( Table 1). The air extraction was carried out using different set-ups based on the principle described by Schwander et al. (1993). Air was pumped out of the firn after a stepwise drilling of the borehole in intervals of a few metres. At each depth, the hole was sealed by an inflatable rubber bladder through which a tube was introduced to pump out firn air to the surface into stainless steel, glass or aluminum containers.
The containers were analyzed for δ 13 C(CH 4 ) in 4 different laboratories: Institute for Atmospheric and Marine research Utrecht (IMAU), Laboratoire de Glaciologie et Géophysique de l'Environnement (LGGE), the Commonwealth Scientific and Industrial Research Organisation (CSIRO) and the Centre for Ice and Climate in Copenhagen (CIC) ( Table  1). The different analytical systems are described in Francey et al. (1999);Bräunlich et al. (2001), Sowers et al. (2005), Röckmann (2011), Sapart et al. (2011). The offsets between measurements carried out in the different laboratories are discussed in Sect. 7 and in the Supplement.
The error bars assigned to individual data points in this paper generally represent the standard deviation of several measurements on the same firn air sample. Measurements of other trace gases on similar samples as those measured for δ 13 C(CH 4 ) lead to the conclusion that possible sampling artefacts such as fractionation during the firn extraction or leaks of the bladders are negligible at the investigated sites.

Data
The δ 13 C(CH 4 ) firn air results are separately presented for the Northern Hemisphere (NH) sites and the Southern Hemisphere (SH) sites in Fig. 1. At each location, δ 13 C(CH 4 ) shows a stable or slowly decreasing trend with depth in the upper firn and stronger isotope depletions in the bubble close-off zone (deeper firn). The firn structure and thickness depend on the snow accumulation rate and on temperature at each site (Table 1) hence these parameters cause differences in the depth profiles observed among the sites (Fig. 1). For Table 1. Information on the eleven firn air pumping sites used in this paper. The NH sites are Devon Island (DI), North GRIP (NGR), NEEM-2009 (NM-09), NEEM-EU hole (NM-EU-08). The SH sites are DE08-2 (DE08), Berkner Island (BI), Vostok (VOS), Dronning Maud Land (DML) and Dome Concordia (DC) (Bräunlich et al., 2001), South Pole 1995 and South Pole 2001 (SPO-01) (Sowers et al., 2005).
(1) Institute where the measurements used in this study were performed. (2) Date of the firn sampling campaigns. (3) Snow accumulation rate. (4) Annual mean temperature. (5) Z lowest is the depth of the lowest measurement for δ 13 C. (6) Estimated length of the δ 13 C scenario.  The SH sites are DE08-2 (DE08), Berkner Island (BI), Vostok (VOS), Dronning Maud Land (DML) and Dome Concordia (DC) (Bräunlich et al., 2001), South Pole 1995 and South Pole 2001 (SPO-01) (Sowers et al., 2005). The numbers associated with the names on the legend are the dates (year) of the firn sampling for each site. the NH sites, the differences are small between the Greenland sites (NEEM and NGR), because the surface conditions (temperature, snow accumulation, etc.) are rather similar at these sites. The Devon Island (DI) dataset shows a different pattern, as surface conditions differ significantly from the Greenland sites. Numerous melt layers have been detected in the DI firn column (e.g. Clark et al., 2007;Martinerie et al., 2009), which strongly reduce molecular diffusivity with depth. The SH data show larger differences between the sites, because the meteorological/glaciological variables (temperature, snow accumulation, snow morphology, etc.) in Antarctica differ strongly from site to site.

SITES
Firn air was sampled in different years (Table 1), but the upper firn data, corrected for seasonality (see Supplement) for the NH sites, do not show a systematic trend as a function of the sampling date. The similarity of near-surface δ 13 C(CH 4 ) values at all sites ( Fig. 1) already suggests no large time trend in δ 13 C(CH 4 ) between the oldest (1993) and youngest (2009) sampling campaigns, in qualitative agreement with direct atmospheric measurements Monteil et al., 2011).

Modeling trace gas transport in firn
To convert vertical δ 13 C(CH 4 ) depth profiles in the firn into temporal atmospheric isotope scenarios, we use a modelling approach based on the mathematical framework described in Rommelaere et al. (1997) and Witrant et al. (2012). The physical model uses a historic evolution of CH 4 atmospheric concentrations to calculate vertical profiles of concentrations in firn. In the case of isotopic ratios, two simulations are performed for the major and minor isotopologues, 12 CH 4 and 13 CH 4 , respectively. Besides atmospheric scenarios, the site temperature, accumulation rate, depth of the convective layer in the upper firn, and the close-off depth are used to constrain the model together with profiles of firn density and effective diffusivity. The effective diffusivity is determined as a function of depth for each firn drilling site by modeling the concentrations in firn for trace gases with well-known recent histories (Rommelaere et al., 1997;Trudinger et al., 1997). Here we used a multi-gas constrained inverse method  to optimise the effective firn diffusivity profile at each site. It is important to note that the diffusivity is not constrained equally well at all sites. This is firstly due to the fact that a different number of species has been used for diffusivity reconstruction for each site, and secondly because firn data were obtained with different depth resolutions (less than 1 m to ∼ 10 m sampling intervals).
Once the diffusivity is obtained, another version of the model based on the forward model physics was used to reconstruct atmospheric time trend scenarios from firn concentration data. We use the Green's function approach introduced in firn modeling by Rommelaere et al. (1997), which has been recently extended to isotopic ratios (Wang et al., 2012) to calculate the probability of having air of a certain age at a certain depth.

Single-site δ 13 C scenarios
In order to investigate the consistency among the different sites, the inverse model is first applied to the isotope data from each site individually and single site δ 13 C(CH 4 ) atmospheric trends are shown in Fig. 2a and d. A regularization term is used to control the smoothness of the scenario and find a unique solution (see Wang et al., 2012 and Supplement). Figure 2b and e show the comparison of the measurements with modeled δ 13 C(CH 4 ) firn profiles using the scenarios shown in panels a and d. The differences between the modeled and measured firn air datasets (panels c and f) are in the majority of the cases similar to or smaller than the assigned measurement uncertainties. This shows that mathematically the inversion set-up performs well, i.e. the reconstructed atmospheric δ 13 C(CH 4 ) scenarios produce firn air isotope profiles, which are in agreement with the measurements.
Strong deviations are observed for two data points only, the DI value at 54 m depth, and the SPO-01 value at 53 m depth. We consider the anomalous data point at SPO-01 as an outlier, because it occurs in the upper firn, where fast diffusion combined with the long timescales of isotope variations is inconsistent with strong isotope variability. Removing this point does not significantly affect the reconstructed scenario. The numerous melt layers (up to 7 cm of thickness) identified in the DI firn column may imply that horizontal diffusion is important, because gases may have to travel laterally around the impermeable layers to reach the deeper firn. Low permeability layers and horizontal transport cannot be represented in our 1-D firn model and the use of a 3-D model goes beyond the state of the art of firn modeling, because to date too few constraints exist to describe 3-D transport in firn. In addition, the model assumes the steady state of the firn structure, whereas at DI the diffusivity/depth profile probably changed with time due to the unequal distribution of melt layers with depth. DI will thus not be used in the final multi-site scenario reconstruction.
Whereas the inversion technically adequately fits the measurements at individual sites, an important result of these single site inversions is that the atmospheric histories re-constructed from individual sites differ strongly among each other, as seen in Fig. 2a and d. Small differences could occur between the δ 13 C(CH 4 ) trends at NH and SH high latitudes, due to the longer time required for isotopic ratios than mixing ratios to adapt to a new source configuration (Tans, 1997) or due to changes in the location of the different CH 4 sources. On the other hand, the scenario reconstructed with the SPO-01 data shows very large discrepancies (more than 2 ‰) compared to other SH sites. The DC and SPO-95 based scenarios are flatter than the previously published scenarios for these sites (Bräunlich et al., 2001;Sowers et al., 2005), whereas our DML based scenario is consistent with the steepest slope scenarios in Bräunlich et al. (2001). Moreover, our SPO-01 scenario is steeper than the scenario range in Sowers et al. (2005). These two studies used the Rommelaere et al. (1997) model with singlespecies-based diffusivity tuning, and a Monte-Carlo technique for scenario reconstruction. Here we used the Witrant et al. (2012) model, which uses similar firn physics, but has a more accurate multi-species-based diffusivity tuning. The differences between our δ 13 C(CH 4 ) single-site scenarios and the published scenarios for δ 13 C(CH 4 ) probably result from the use of different diffusivity profiles and modeling approach. Indeed in a recent firn model intercomparison study, Buizert et al. (2012), conclude that diffusive fractionation of isotopes is insufficiently constrained by firn models and that using different firn air models can result in significantly different reconstructions. We further investigate this issue in Sects. 5.3 and 5.4.

Effect of firn fractionation
The presence of a trend in CH 4 atmospheric mixing ratio implies a transport of the gas within the firn column which alters the isotopic signature of this gas. This is due to the fact that the difference of the molecular diffusion coefficient in air between the two isotopologues 12 CH 4 and 13 CH 4 produces large signals of δ 13 C(CH 4 ) in the firn column, even in the absence of a δ 13 C(CH 4 ) trend in the atmosphere. Therefore it is important to evaluate the amplitude of the combined effects of isotopic fractionation in firn due to molecular diffusion and gravitational settling using a firn model run in a forward mode (Trudinger et al., 1997). An atmospheric CH 4 scenario was built based on the NOAA-ESRL atmospheric network and ice core data (Etheridge et al., 1998;MacFarling Meure et al., 2006). To assess the role of firn fractionation only, this scenario was used as input, assuming a constant atmospheric δ 13 C(CH 4 ) history, to model 12 CH 4 (Buizert et al., 2012 for the NH and Witrant et al., 2012 for the SH), and 13 CH 4 in the firn column. The calculated 12 CH 4 and 13 CH 4 mixing ratios in firn air are then recombined to yield a δ 13 C(CH 4 ) depth profile which is only due to isotopic fractionation in firn and not to changes in the δ 13 C(CH 4 ) atmospheric history.
Results are shown in Fig. 3 as the difference between modeled δ 13 C(CH 4 ) at the surface and at each model depth level, and compared with the same difference in the firn data (δ 13 C(CH 4 ) − δ 13 C(CH 4 ) surface ). Gravitational fractionation produces a slight upward trend in δ 13 C(CH 4 ) in the diffusive part of the firn column (see Supplement, Fig. S1), whereas molecular diffusion in combination with the CH 4 mixing ratio gradient along the firn profile produces a clear trend in δ 13 C(CH 4 ) towards more depleted values in the deep firn, the amplitude being site-dependent (Fig. 3).
At most sites, isotopic fractionation in firn alone can already explain a large fraction of the trends in the observed δ 13 C(CH 4 ) depth profiles. This implies that only a small positive trend in atmospheric δ 13 C(CH 4 ) is needed to fit the measurements (Fig. 3). For DI, DML, and SPO-01, however, the difference between the constant δ 13 C(CH 4 ) scenario and the data is large. Therefore, the scenarios reconstructed from these sites show large atmospheric trends in comparison with the other sites as shown in Fig. 2.
Although firn data do not suggest strong differences in the firn physical properties between the two drill holes at the South Pole (1995 and 2001), the model produces very different estimates of isotopic fractionation in firn and the possible reason for this is discussed in the next section.

Sensitivity to small changes in diffusivity
The discrepancy between the reconstructed histories at the South Pole from the two firn air campaigns in 1995 and 2001 is puzzling and we investigated possible causes. An important input parameter in firn models is the effective diffusivity profile (e.g. Buizert et al., 2012), which is constructed for each site and drilling operation from an inversion approach using trace gas mixing ratios with relatively well-known 35 Figure 4: δ 13 C(CH4) firn fractionation calculated with regular (full lines) and modified diffusivities (dashed lines). a: NH sites, NM--EU--08 (purple dots) run with regular diffusivity (full line), with NM--09 diffusivity (short--dashed lines) and with NM--US--08 (long--dashed lines) diffusivity. b: SH sites, SPO--01 (light blue), DML (black) and DEO8 (orange) run with regular diffusivities (full lines) and with modified diffusivity: DEO8 run with DEO8 "slightly modified" diffusivity (see text, orange dashed--line), DML run with δ 13 C(CH4) constrained diffusivity (black dashed--line) and SPO--01 run with SPO--95 (light blue dashed--line). atmospheric histories, as explained in Sect. 4. For SPO-95, the diffusivity is constrained by a combination of 6 species (CO 2 , CH 4 , SF 6 , CH 3 CCl 3 , CFC-11, CFC-12) whereas data for only 3 species (CO 2 , CH 4 , SF 6 ,) are available for SPO-01, hence these diffusivities are not equally well constrained . In their analysis Sowers et al. (2005) only used one species, CO 2 , to deduce the diffusivity profile used for both SPO-95 and SPO-01. The CH 4 depth profile only provided an a posteriori test of the diffusivities derived from CO 2 .
To investigate the sensitivity of δ 13 C(CH 4 ) to small differences in the diffusivity profile, firn fractionation was simulated for SPO-01 using the better constrained diffusivity profile of SPO-95 (Fig. 4b). The rationale behind this approach is that the diffusivity profiles at nearby drilling locations (with density profiles and close-off depths nearly similar) are expected to be comparable even if the firn core is drilled in a different year. A diffusivity profile from a nearby firn drilling site where more trace gas profiles are available may then be more realistic than the profile from a certain lo-cation where only few species are available. Figure 4b shows that this slight change in diffusivity apparently has a large effect on the δ 13 C(CH 4 ) firn fractionation. This change also removes much of the discrepancy that was apparent between SPO-01 and SPO-95. It is important to note that the model runs with the modified diffusivity still produce acceptable agreement with the CO 2 , CH 4 and SF 6 mixing ratios (the tracers used for determining the diffusivity profile) measured at SPO-01 .
A similar comparison was made for the two different NEEM boreholes drilled in 2008 and 2009. Although the single site isotope reconstructions from these sites are not as different as for SPO (Fig. 2), we calculated firn fractionation of the 2009 borehole with the diffusivity profile reconstructed from the 2008 hole. Diffusivity for the NM-09 borehole was constrained by 2 species (CO 2 and CH 4 ) and diffusivity for NM-08 by 9 species (CO 2 , CH 4 , SF 6 , CH 3 CCl 3 , CFC-11, CFC-12, CFC-13, HFC-134a, 14 CO 2 ). As an additional test we used the diffusivity of NM-US-08, but constrained by three species (CO 2 , CH 4 , SF 6 ) to estimate the isotopic fractionation in firn of NM-EU-08. This leads to a smaller fractionation below ∼ 60 m than using the NM-09 diffusivity ( Fig. 4a: purple long-dashed line). This means that even though the diffusivity profiles of NM-US-08 and NM-09 were constrained by almost the same species, they do show significant differences. These diffusivity differences might be due to lateral inhomogeneities in the firn. Due to the high sensitivity of δ 13 (CH 4 ) to diffusivity, the uncertainties on diffusivity profiles constrained with only a few gases (Trudinger et al., 2012) may play a major role.
Additionally, we compared the effect of diffusivity for two Antarctic sites with high (DE08) and low (DML) snow accumulation rates (Table 1). High accumulation sites are expected to be less affected by isotopic fractionation in firn due to a stronger advection and faster gas trapping in bubbles. For DEO8, we used two different diffusivity profiles (constrained with/without the CO 2 outlier at 80 m depth, Witrant et al., 2012) and this leads to a small difference in δ 13 C(CH 4 ) (∼ 0.2 ‰) in the calculated isotopic fractionation at 80 m depth (Fig. 4b: yellow full and dashed lines). This indicates that for high accumulation sites, changing the diffusivity has a limited effect on isotopic fractionation.
A very small firn fractionation estimate is obtained at DML, despite the fact that it belongs to the low accumulation rate category. This results in a very strong change in atmospheric δ 13 C(CH 4 ) with time in the single-site (DML) reconstruction, which is not in agreement with the isotope histories at most other sites (Fig. 2). Here, we do not have an alternative diffusivity profile from a similar site. We thus attempted to use a data based δ 13 C(CH 4 ) atmospheric trend as an additional constraint in the construction of the diffusivity profile. We used a δ 13 C(CH 4 ) scenario based on our multi-site reconstruction and ice core data from Ferretti et al. (2005), in addition to six other gases species (CO 2 , CH 4 , SF 6 , CH 3 CCl 3 , CFC-12, CFC-113) to re-tune the diffusivity Atmos. Chem. Phys., 13, 6993-7005, 2013 www.atmos-chem-phys.net/13/6993/2013/ at DML. As expected, the simulations with this modified diffusivity profile lead to a considerably stronger isotopic fractionation along the DML firn column (Fig. 4b: black dashed line). It is important to note that all other gases are still fitted well with this modified diffusivity profile (not shown). The large uncertainty of the isotopic fractionation in the deep DML firn may be due to an insufficient sampling resolution below 70 m depth (2 depth levels), where steep concentration gradients are observed for all species. This example of using δ 13 C(CH 4 ) as a constraint confirms the strong sensitivity of δ 13 C(CH 4 ) to the firn diffusivity profiles. This implies that -once a reliable isotope history has been establishedδ 13 C(CH 4 ) may be a valuable constraint to obtain realistic firn diffusivity profiles. In this case, the "weakness" of δ 13 C(CH 4 ), namely that the isotopic fractionation in firn is of the same order of magnitude or even larger than the atmospheric changes, would turn into an advantage.
In summary, these results show that model estimates of isotopic fractionation in firn are strongly dependent on the diffusivity profile used, similar to the findings of Buizert et al. (2012). Although constrained in our case with multiple gas histories, the firn diffusivity is not always sufficiently constrained for an accurate estimate of the firn fractionation for δ 13 C(CH 4 ), at least in the deep firn, and thus for atmospheric reconstruction of δ 13 C(CH 4 ) from a single site.

Multi-site δ 13 C reconstruction
We have shown that trying to reconstruct atmospheric scenarios of δ 13 C(CH 4 ) from single sites leads to inconsistent solutions. In the following we carry out a multi-site reconstruction by performing two inversions for all sites from the SH and NH, respectively (Fig. 5). In the NH multi-site inversion, the DI data are excluded as the firn column at this site is affected by melt layers, which likely caused the inconsistency with other sites as shown by the single site reconstructions (Fig. 2). In the SH multi-site inversion, the SPO-95 diffusivity, constrained with 6 species instead of only 3 for SPO-01 (see Sect. 5.3), is used for the SPO-01 firn. The deepest firn data points provide constraints further back in time than shallower data points in the bubble close-off zone, thus potentially extending the historical reconstruction of atmospheric δ 13 C(CH 4 ) back to the early 20th century (Sowers et al., 2005). However, the correction for isotopic fractionation in firn is most uncertain for the deepest samples, where strong differences between individual firn air models have been reported (Buizert et al., 2012). Therefore the deep firn data points showing strong deviations with the modeled isotopic fractionation (Fig. 3) were excluded. As the NM-09 diffusivity is constrained by only two gases (CO 2 and CH 4 ) instead of nine gases for NM-EU, a larger number of data are eliminated for NM-09. The excluded data points are shown in grey in Fig. 5. As this data elimination procedure is somewhat subjective, an alternative method was also used where all data below the lock-in depth (LIDgas in Table 2 of Witrant et al., 2012) were eliminated. The two methods lead to very similar results (see Supplement, Fig. S4). The exclusion of the oldest firn air samples restricts the length of the scenario constrained by the remaining firn data to the last ∼ 50 yr (mean ages at the last depth levels used are provided in Table 1), but it provides a more robust scenario.
In the preferred multi-site inversion, each firn profile is given equal weight (red lines in Fig. 5a and d). In a second inversion (green lines in Fig. 5a and d), the contribution of each site is weighted according to the inverse of the RMS (root mean squared) difference of the model results from single site inversions (Fig. 2) and the measurements. This weighting reduces the relative weight of sites that are affected by larger experimental uncertainties. It leads to slightly flatter δ 13 C(CH 4 ) scenarios in the latter multi-site inversion, but the results in Fig. 5 show that the differences are well within the uncertainty envelopes of the original multi-site inversion. The fact that the green and red reconstructions are similar indicates that the multi-site reconstructions are not affected by individual sites with obviously poor data quality. Due to the difficulty to precisely evaluate overall uncertainties, the equal-weight approach is the preferred multi-site reconstruction. Figure 5b and e show how the modeled firn air δ 13 C(CH 4 ) profiles, based on the multi-site atmospheric reconstructions, agree with the experimental data. The right hand panels show the model-measurement differences. Overall, the differences are of the order of the analytical uncertainty of 0.05-0.3 ‰.
As a further consistency check of the individual firn air measurements with the multi-site δ 13 C(CH 4 ) scenario, we place isotope measurements in firn on a timescale, using the mean age of the samples and taking into account the isotopic fractionation in firn. This simplified approach compared to inverse methods was also used in, for example, Francey et al. (1999) and Ferretti et al. (2005). It is less accurate in terms of gas age and trend determination than inverse scenario reconstruction but does not require a regularization term to ensure the uniqueness of the solution. Figure 6 shows that the data converted with this simplified approach fall within the uncertainty range of the inverse reconstruction, thus the two methods are overall consistent. Note that for the readability of the figure, the age spread is not represented in Fig. 6, but that no discrete age exists for a given firn air sample. The lower panel shows that the data from DE08 and SPO-01 are near the lower envelope of the multi-site reconstruction, whereas the data from VOS are for a large part of the record near the higher end. This reflects the results from the single site reconstructions, where the VOS data required only a weak δ 13 C(CH 4 ) atmospheric trend, whereas DE08 and SPO-01 required a strong δ 13 C(CH 4 ) atmospheric trend. Figure 5: Multi--site δ 13 C(CH4) trend reconstructions (a and d), the model fit of the firn data (b and e) and the difference between model and data (c and f). NH sites (a, b and c) and SH sites (d, e and f). Two multi--site inversions are shown in a) and d), one giving equal weight to each site (red lines) and one weighting each site by the RMS difference of the single site reconstruction (Fig. 2) to the data (green lines). The light grey line (left panel) represents air--archive results from Cape Grim (Francey et al., 1998) and continuous atmospheric monitoring data from the NOAA--ESRL network.

!" #"
Atmospheric scenarios Fit to firn data Differences model-data (equal-weight)  (d), one giving equal weight to each site (red lines) and one weighting each site by the RMS difference of the single site reconstruction (Fig. 2) to the data (green lines). The light grey line (left panel) represents air-archive results from Cape Grim (Francey et al., 1998) and continuous atmospheric monitoring data from the NOAA-ESRL network.
The multi-site reconstructions are constrained by δ 13 C(CH 4 ) data from 10 sites presented in Table 1 (i.e. excluding DI). For the NH sites, the multi-site scenario is constrained by only a few data points before 1990. Furthermore, the firn profiles at the three sites that enter the multi-site inversion (NGR, NM-EU-08 and NM-09) are alike, as shown by the similar vertical profiles in Fig. 1. Therefore, the NH scenario is less constrained than the SH scenario, where data from sites with very different firn characteristics and sampling dates are available. The lack of independent constraints and the relatively small number of data points may explain the higher uncertainty range of the NH δ 13 C(CH 4 ) reconstruction.
Due to its high snow accumulation rate, DE08 is less affected by firn fractionation than other sites (see Fig. 3). An atmospheric trend scenario, based on DE08 firn and ice data, was built and compared to firn data from the SH (see Sect. 6 of the Supplement and Fig. S5). For the two sites that were sampled relatively close in time to DE08, namely SPO-95 and VOS, the model results show clear differences to the measurements, also well above the LID (see also Fig. S5b). This suggests that firn fractionation may not be the only explanation for the discrepancies between our firn drilling sites (see Sect. 7).
Although the multi-site scenarios for the SH and NH agree within the uncertainty envelopes when we shift the SH data by a constant IPG of 0.5 ‰ (Fig. 7a), the shape of the best estimate scenario is slightly different for the two hemispheres. Due to the fact that the main anthropogenic CH 4 sources are in the NH, it cannot be excluded that the inter-polar gradient (IPG) of δ 13 C(CH 4 ) changed, but given the short inter-hemispheric mixing time of 1 yr, large variations on short timescales are unlikely. In any case, given the large uncertainty envelopes on the multi-site reconstructions it is not possible to extract robust information about the IPG.

Comparison of our firn reconstructions with atmospheric and ice core data
In Fig. 7, we compare our firn air reconstructions with direct and archive atmospheric measurements, and ice core data. Continuous atmospheric measurements of δ 13 C(CH 4 ) have been carried out over the last ten years in both hemispheres Monteil et al., 2011). For the NH, no direct atmospheric data are available before 2000. For the SH, δ 13 C(CH 4 ) was determined from archived air samples collected at Cape Grim, Australia from the late 1970s (Francey et al., 1999). Moreover, ice core data from Law Dome, Antarctica (Ferretti et al., 2005) overlap in time with the period of our multi-site reconstruction. For the NH, no δ 13 C(CH 4 ) ice core data are available for the reconstructed period, and we show in Fig. 7a a comparison with five data points from EUROCORE (Sapart et al., 2012) covering the period 1800-1950. Figure 7 shows that the direct measurements over the last decade (light grey lines) fall within the uncertainty envelopes of the δ 13 C(CH 4 ) reconstruction from the firn air samples. The Cape Grim archive data from the SH, which are not affected by firn air fractionation (dashed grey line in Fig. 7b), show a steeper temporal trend than our reconstruction. Interlaboratory calibration offsets can considerably hamper direct Atmos. Chem. Phys., 13, 6993-7005, 2013 www.atmos-chem-phys.net/13/6993/2013/ C. J. Sapart et al.: Carbon isotopic composition of methane 7001 37 Figure 6: Firn air δ 13 C(CH4) datasets corrected for isotopic fractionation in firn (diffusion and gravitational settling) as a function of age. A: NH sites, NM--EU--08 (purple), NM--09 (red), NGR (green). b: SH sites, DE08 (orange), BI (purple), SPO--95 (dark blue), SPO--01 (light blue), DML (black), DC (green), VOS (brown). The dashed grey line in b) represents air--archive data from Cape Grim (Francey et al., 1998) and the solid grey lines represent continuous atmospheric monitoring data from the NOAA--ESRL network. Note that the age spread of each data point is not presented here in order to keep the figure readable.  (Francey et al., 1998) and the solid grey lines represent continuous atmospheric monitoring data from the NOAA-ESRL network. Note that the age spread of each data point is not presented here in order to keep the figure readable. comparison of the datasets and an absolute offset could potentially be attributed to uncertainties in the calibration scales. However, the Cape Grim data also show a somewhat stronger δ 13 C(CH 4 ) trend than the multi-site firn air scenario. Experimental uncertainties are further discussed in the next section. One other possible parameter that can influence the model results is the uncertainty in the best-guess CH 4 input scenario. A sensitivity test (Supplement, Fig. S2) shows that using atmospheric CH 4 values at the upper end of their uncertainty envelope may lead to a slope nearly consistent with the one from the Cape Grim air archive.
The same interpretation holds for the firn and ice core data from Ferretti et al. (2005) from Law Dome (Fig. 7b). The Law Dome data suggest an isotope change of 1.5 ‰ over the past 50 yr, about twice as large as derived from our multi-site firn reconstruction, which is in agreement with the single site scenario reconstruction based on DE08 firn data only shown in Fig. 2. Ice core data are also affected by isotopic fractionation effects in the overlying firn, and the data published in Ferretti et al. (2005) were corrected for these effects using the CSIRO firn air model (Fig. 7: orange dots). To be consistent with our own reconstructions, the corresponding corrections 38 Figure 7: Multi--site δ 13 C(CH4) scenarios compared with direct atmospheric measurements (solid grey lines), atmospheric archive data from Cape Grim (dashed grey line) and ice core data (dots). The red continuous lines are the "best--estimate" scenarios for each hemisphere with uncertainties (short dashed lines). The red long--dashed--line is the SH scenario shifted by 0.5‰ a: NH raw δ 13 C(CH4) data (light blue dots) and δ 13 C(CH4) data corrected for firn fractionation from the EUROCORE ice core (Sapart et al. 2012). b: SH raw δ 13 C(CH4) firn and ice data at DEO8 (yellow dots) and the same data corrected for firn fractionation with the CSIRO model (orange dots) (Francey et al., 1999, Ferretti et al., 2005, with the LGGE--GIPSA model using regular diffusivity (black stars) and with the LGGE--GIPSA model but with "slightly modified diffusivity" as in Fig.  4  (red  stars). The grey lines represent continuous atmospheric monitoring data from the NOAA--ESRL network (NH and SH).

!"
#" Fig. 7. Multi-site δ 13 C(CH 4 ) scenarios compared with direct atmospheric measurements (solid grey lines), atmospheric archive data from Cape Grim (dashed grey line) and ice core data (dots). The red continuous lines are the "best-estimate" scenarios for each hemisphere with uncertainties (short-dashed lines). The red long-dashed line is the SH scenario shifted by 0.5 ‰ (a) NH raw δ 13 C(CH 4 ) data (light blue dots) and δ 13 C(CH 4 ) data corrected for firn fractionation from the EUROCORE ice core (Sapart et al., 2012). (b) SH raw δ 13 C(CH 4 ) firn and ice data at DEO8 (yellow dots) and the same data corrected for firn fractionation with the CSIRO model (orange dots) (Francey et al., 1999;Ferretti et al., 2005), with the LGGE-GIPSA model using regular diffusivity (black stars) and with the LGGE-GIPSA model but with "slightly modified diffusivity" as in Fig. 4 (red stars). The grey lines represent continuous atmospheric monitoring data from the NOAA-ESRL network (NH and SH).
were also calculated with the LGGE-GIPSA model using regular diffusivity (Fig. 7: black stars) and "slightly modified diffusivity" as explained in Sect. 5.3. (Fig. 7: red stars). The corrected firn air data with the CSIRO model (orange dots) and with the LGGE-GIPSA model with modified diffusivity (red stars) show similar results which suggests that the difference in firn air models is not responsible for the differences between the Law Dome records and our SH firn reconstruction. The differences between the raw and firn fractionationcorrected δ 13 C(CH 4 ) shown in Fig. 7 demonstrate that the DE08 young ice core data are also significantly affected by isotopic fractionation in the firn. This is expected since these samples contain air from a period when CH 4 mixing ratios in the atmosphere strongly increased.
The same is true for the only available ice core data from the NH site EUROCORE (Fig. 7a). However, here the δ 13 C(CH 4 ) data corrected for isotopic fractionation are rather Figure  8: NH Multi--site δ 13 C(CH4) trend scenarios based on Kr corrected firn data. (a) Note that here NGR data are data measured at IMAU with lower depth resolution than previously shown, because the NGRIP data from LGGE could not be corrected for the Kr interference. (b) the fit to the firn data and (c) model--data differences. The continuous grey curve (a) represents continuous atmospheric monitoring data from the NOAA--ESRL network. The dashed grey curve shows the preferred multi-site scenario based on non Kr--corrected data (red continuous line in Fig.5a).

Atmospheric scenarios
Fit to firn data Differences model-data (equal-weight) !" #" !"#$%&&!"#'(#$)#*"+(&&!,-&&&& $" Fig. 8. NH Multi-site δ 13 C(CH 4 ) trend scenarios based on Kr corrected firn data. (a) Note that here NGR data are data measured at IMAU with lower depth resolution than previously shown, because the NGRIP data from LGGE could not be corrected for the Kr interference. (b) the fit to the firn data and (c) model-data differences. The continuous grey curve (a) represents continuous atmospheric monitoring data from the NOAA-ESRL network. The dashed grey curve shows the preferred multi-site scenario based on non Kr-corrected data (red continuous line in Fig. 5a).
at the high end of the firn air reconstruction. The firn/ice data comparison for the Northern Hemisphere is limited by the absence of overlapping period but the large difference between the EUROCORE and Law Dome data at 1950 additionally implies a possible unidentified interlaboratory offset between these datasets.

Additional model and measurement uncertainties
As discussed above, the discrepancies between the different trend reconstructions may be due to calibration issues, uncertainties of the CH 4 scenario and/or insufficiently constrained diffusivity, but also to other features.
A major drawback while working with multiple datasets measured at different times and in different laboratories is the uncertainty associated with intercalibration differences between measurement systems, but also reference gas drifts over time which may affect datasets measured in the same laboratory. Sapart et al. (2011) suggested that the intercalibration should be carried out using at least 3 reference datapoints with different δ 13 C(CH 4 ) values to identify and quantify possible system offsets and scale contractions. This exercise could not be carried out here, because most of the datasets presented were measured long ago and some systems have been dismounted. The available information about intercalibration is presented in Table S1.
It has been discovered recently that Krypton (Kr) produces a strong interference for CH 4 isotope measurements in many of the currently used analytical systems (Schmitt et al., 2013). When the Kr/CH 4 ratio of the sample differs from the Kr/CH 4 of the reference air, which is generally the case for firn air over the last century and for ice air samples, this introduces a strong bias in the δ 13 C(CH 4 ) measurements. This system specific interference introduces a strong bias dependent on both the concentration of the air sample and the analytical set-up. Considering the CH 4 concentration trend (about 800 ppb) between the lock-in zone and the firn surface, the bias induced by the Kr interference could reach −0.8 ± 0.2 ‰ (Schmitt et al., 2013). For the IMAU system, it was possible to correct for the Krinterference (Schmitt et al., 2013). Using this correction, we carried out a δ 13 C(CH 4 ) reconstruction for the NM-EU-08, NM-09 and NGR datasets with Kr correction (Fig. 8). The temporal trend of the Kr-corrected scenario is about 0.5 ‰ larger than the non-Kr corrected scenario showing the significance of the Kr-correction on δ 13 C(CH 4 ) trend reconstruction. Unfortunately, the Kr-correction is not possible for the SH datasets, which limits the interpretation of our SH reconstruction in term of atmospheric trends. Buizert et al. (2012) showed that different firn air models produce isotopic fractionation effects of different magnitude. For δ 13 C(CO 2 ), the LGGE-GIPSA model produces the largest diffusive fractionation below 65 m depth at NM-EU-08, because a dispersive mixing term is not considered in the lock-in zone. Sensitivity tests (Sect. 7 of the Supplement) show as expected that including a dispersive mixing term reduces the firn fractionation and thus requires larger atmospheric trends. The effect of dispersive mixing for the Greenland sites is consistent with the estimate of Buizert et al. (2012). However, the discrepancies between the SH sites (see Sect. 5.4) cannot be easily reduced because they do not only affect the lock-in zone but also the diffusive zone. Adding one common dispersive parameter in the diffusivity profiles thus cannot solve the discrepancies between sites. Dispersive mixing tends to reduce firn fractionation; however, as detailed in Sect. 5.3, we suspect isotope fractionation to be underestimated at three sites (DI, DML, SPO-01) where the diffusivity profiles are less constrained than at the other sites. At the better constrained DE08 site, isotopic fractionation along the firn column and in ice bubbles from our model is consistent with the former CSIRO model and yields the same trend in δ 13 C(CH 4 ).
To our knowledge, the LGGE-GIPSA model is the only model available that allows multi-site inversions including a large number of sites. It remains to be investigated whether models with different representation of firn physics would yield more consistent results between available firn and young ice core data, but the experiment above indicates that this may not be the case. Experimental issues such as the Kr interference and inter-laboratory offsets, and the sensitivity of the isotope signals to small changes in diffusivity are presently limiting factors that lead to rather large uncertainties in reconstructing an atmospheric δ 13 C(CH 4 ) scenario from firn air samples.

Conclusion and perspectives
We reconstructed δ 13 C(CH 4 ) atmospheric scenarios over the last 50 yr for both northern and southern high latitudes using a firn transport model and taking as input δ 13 C(CH 4 ) firn air measurements from 10 firn sampling campaigns. Our calculations reveal discrepancies between sites and show that the modeled isotopic fractionation between 12 CH 4 and 13 CH 4 along the firn column is highly sensitive to slight variations in the diffusivity profiles at each site.
After plausible adjustments of diffusivity profiles at individual sites, and neglecting the deepest firn air samples (where the isotopic fractionation is strongest), it is possible to reconstruct an estimate of δ 13 C(CH 4 ) for each hemisphere from a multi-site inversion, which fits the selected measurements within reasonable uncertainties; however, discrepancies between the different datasets and with previously published results remain. At present, it is difficult to discriminate between different sources of uncertainty such as scale differences between laboratories, the recently discovered Krinterference affecting the δ 13 C(CH 4 ) trends, the representation of fractionation effects in the lock-in zone and the uncertainty of the 12 CH 4 scenario used. After the recent discovery of a Kr interference, analytical systems will need to be revised. This could eliminate an important source of discrepancies between datasets in the future.
Our analysis showed that results of firn air models for δ 13 C(CH 4 ) strongly depend on accurate determination of the diffusivity profiles. Sampling at high-accumulation rate sites, such as DEO8 is an advantage, because firn fractionation is less strong at such sites, hence the δ 13 C(CH 4 ) reconstruction will be more robust. Once a reliable isotope history has been constructed, δ 13 C(CH 4 ) data themselves may be used in the future to construct diffusivity profiles better than what can be achieved with currently used species.