Reconstruction of the carbon isotopic composition of methane over the last 50 yr based on firn air measurements at 11 polar sites

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 reconstructing the evolution of its sources and sinks in the past. This is important to improve predictions of atmospheric CH4 mixing ratios in the future under the influence of a changing climate. We present an attempt to reconcile methane carbon isotope records from 11 firn sites from both Greenland and Antarctica to reconstruct a consistent 13C(CH4) history over the last 50 yr. In the firn, the atmospheric signal is altered mainly by diffusion and grav itation. These processes are taken into account by firn transport models. We show that isotope reconstructions from individual sites are not always mutually consistent among the different sites. Therefore we apply for the first time a multisite isotope inversion to reconstruct an atmospheric isotope history that is constrained by all individual sites, generating a multisite “best-estimate” scenario. This scenario is compared to ice core data, atmospheric air archive results and direct atmospheric monitoring data.


Introduction
Methane (CH 4 ) is a strong greenhouse gas and plays an important role in atmospheric chemistry. Its atmospheric mixing ratios rapidly increased since 1800 in response to anthropogenic emissions (e.g. MacFarling Meure et al., 2006;IPCC, 2007). Understanding the factors responsible for changes in atmospheric CH 4 is not straightforward, 20 because the CH 4 atmospheric mixing ratio depends on the strength of its numerous sources 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).
Stable isotope analysis represents an excellent tool to investigate changes of individual CH 4 sources and sinks in the atmosphere, because each type of emis- 25 sion/removal process is associated with a characteristic isotope signature (Quay et al.,
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 5 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 mechanisms: molecular diffusion and gravitational settling (Craig et al., 1988;Schwander 10 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 those properties (e.g. Schwander et al., 1993;Buizert et al., 2011, and reference therein) 15 can calculate the evolution of trace gas concentrations in the atmosphere based on depth-dependent physical parameters of the firn (density, porosity, and tortuosity) and site-dependent surface conditions (temperature, precipitation and pressure).
Previous studies have shown the potential of reconstructing the history of the CH 4 budget using mixing and stable isotope ratios measurements from firn air (Etheridge 20 et al., 1998;Francey et al., 1999;Bräunlich et al., 2001;Sowers et al., 2005). Those 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 (noted thereafter: δ 13 C(CH 4 )) during the 20th century, which was ascribed to an increase of anthropogenic enriched (fossil and pyrogenic) CH 4 sources during this period. 25 Here we combine previously published firn air δ 13 C(CH 4 ) results from 5 sites and new results from 6 sites from both Greenland and Antarctica. Our goal is to reconstruct a multisite "best-estimate" atmospheric history for both hemispheres. 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., 2011;Wang et al., 2011). We then discuss the differences among the 11 sites and use a multi-site inversion to retrieve the best-estimate reconstruction of δ 13 C(CH 4 ) from firn air over the last 50 yr. Finally, the math/mismatch of these firn air results with atmospheric data and ice core results is discussed.

2 Firn air sampling and isotope analysis
Firn air samples discussed in this paper come from 11 sites 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 10 intervals of a few meters. 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  (2011), Sapart et al. (2011) and Sperlich et al. (2012). 20 To verify the mutual consistency among datasets from different laboratories, IMAU measured samples from NGRIP (NGR) and Dome C (DC), which were previously analyzed by LGGE (Bräunlich et al., 2001). The average δ 13 C(CH 4 ) difference between IMAU and LGGE is 0.14 ± 0.07 ‰ for which IMAU data have been corrected. The same correction is applied to NEEM 2008 firn samples from the European hole (NM-EU-08), 25 which were measured at both IMAU and CIC in this study. It has been shown previously that there are no systematic differences between LGGE and CSIRO (Bernard, 2004).

9591
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 samples from the same boreholes usually lead to the conclusion that possible sampling artifacts such as fractionation during the firn extraction or leaks of the bladders are negligible on 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 struc-10 ture 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

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 mathematics described in Rommelaere et al. (1997) and Witrant et al. (2011). The physical model uses atmospheric scenarios (best-guess historic evolution of trace gas concentrations) to calcu-5 late 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 10 diffusivity is not affected by the micro-structure of the firn, i.e. pore spreading, size and shape, but is related to molecular diffusion occurring at the macro-structure scale, i.e. the density fluctuations (Fabre et al., 2000). 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., 15 1997). Here we used a multi-gas constrained inverse method (Witrant et al., 2011) to calculate the firn diffusivity 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 20 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 function approach introduced by Rommelaere et al. (1997), which has been recently extended to isotopic ratios (Wang et al., 2011;25 Martinerie et al., 2012) to calculate the probability to have air of a certain age at a certain depth.

9593
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5 Model results 5.1 Single-site δ 13 C scenarios In order to investigate the consistency among the different sites, the inverse scenario model has been applied to the isotope data from each site individually and single site δ 13 C(CH 4 ) atmospheric trends have been determined (Fig. 2). The center panels in 5 Fig. 2 show the comparison of the modeled δ 13 C(CH 4 ) firn profiles using the scenarios shown in the left hand panels and the measurements. 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 setup 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 15 the long time scales 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 likely imply that horizontal diffusion is no longer negligible, because gases may have to travel horizontally around the impermeable layers to reach the deeper firn. Low permeability 20 layers and horizontal transport are not represented in our 1-D firn model. In addition, the model assumes 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 multisite scenario reconstruction.
Whereas the inversion technically adequately fits the measurements at individual 25 sites, an important result of these single site inversions is that the atmospheric histories reconstructed from individual sites show a large variability, as seen in the left panels of Fig. 2. Differences could be expected between the δ 13 C(CH 4 ) trends at NH 9594 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). On the other hand the scenario based on SPO-01 shows large discrepancies (more than 2 ‰) in comparison with other SH sites. DC and SPO-95 based scenario are flatter than the previously published scenarios for these sites (Bräunlich et al., 2001;Sowers et al., 2005), whereas 5 our DML based scenario is consistent with the steepest slope scenarios in Bräunlich et al. (2001). On the other hand, 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 single species based diffusivity tuning, and a Monte-Carlo technique for scenario reconstruction. Here we used the Witrant et al. (2011) model, which has a supposedly 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 firn model intercomparison study, Buizert et al. (2011), conclude that diffusive fractionation of isotopes is insufficiently constrained by firn models and that using different 15 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 difference of molecular diffusion coefficient in air between the two isotopologues 12 CH 4 and 13 CH 4 produces large changes of δ 13 C(CH 4 ) in the firn column even in ab-20 sence of δ 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 runs in forward mode (Trudinger et al., 1997). An atmospheric CH 4 scenario based on NOAA-ESRL atmospheric network and ice core data (Etheridge et al., 1998;MacFarling Meure et al., 2006), assuming a con- 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 atmospheric history.
Results are shown in Fig. 3 as the difference between δ 13 C(CH 4 ) at the surface and δ 13 C(CH 4 ) at each sample depth level, and compared with the same difference in the 5 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 and 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 15 show rather 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 South Pole (1995 and, the model produces different estimates of isotopic fractionation in firn and the possible reason for this is debated in 20 the next section.

Sensitivity to small changes in diffusivity
The discrepancy between the reconstructed histories at 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 diffusivity (e.g. Buizert et al., 2011), which 25 is constructed for each site and drilling operation from an inversion approach using trace gas mixing ratios with relatively well-known 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 (Witrant et al., 2011). Sowers et al. (2005) only used the CO 2 atmospheric history and firn air profile to deduce the diffusivity profile used for both SPO-95 and 5 SPO-01. The CH 4 depth profile only qualified a posteriori the CO 2 -constrained diffusivity profile.
To investigate the sensitivity of δ 13 C(CH 4 ) to small differences in the diffusivity profile, atmospheric scenarios were simulated for SPO-01 using the better constrained diffusivity profile of SPO-95 (Fig. 4). The rationale behind this approach is that the dif-10 fusivity 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 location where only few species are available. Figure 4 shows that this slight change in diffusiv-15 ity has apparently a large effect on the δ 13 C(CH 4 ) reconstruction and removes much of the discrepancy that was apparent for the SPO-01 samples in the original inversion.
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 (Witrant et al., 2011). 20 A similar comparison was made for the two different NEEM holes drilled in 2008 and 2009. Although the single site isotope reconstructions from these sites are not as different as for SPO (Fig. 2), we ran the inversion of the 2009 hole with the diffusivity profile reconstructed from the 2008 hole. Diffusivity for the NM-09 hole was constrained with 2 species (CO 2 and CH 4 ) and diffusivity for NM-08 by 9 species (CO 2 , CH 4 , SF 6 , 25 CH 3 CCl 3 , CFC-11, CFC-12, CFC-13, HFC-134a, 14 CO 2 ). The atmospheric δ 13 C(CH 4 ) history for the 2009 reconstruction with the (better constrained) 2008 diffusivity is in good agreement with the original NM-EU-08 results (see Fig. 4, upper panel). Using the diffusivity of NM-US-08, constrained by three species (CO 2 , CH 4 , SF 6 ) to estimate 9597 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the isotopic fractionation in firn of NM-EU-08 leads to a reduced fractionation below ∼60 m ( Fig. 4: purple long dashed-line) than using the NM-09 diffusivity ( Fig. 4: purple long dashed-line). This means that even though the diffusivity profiles of NM-US-08 and NM-09 which were constrained by almost the same species, show significant differences. These diffusivity differences might be due to lateral inhomogenities in the firn or 5 by sampling or analytical artifacts, e.g.: sample contaminations or bladder outgassing (Buizert et al., 2011). Additionally, we compared the effect of diffusivity for two Antarctic sites with high (DE08) and low (DML) snow accumulation rates (Table 1). Indeed, high accumulation sites are expected to be less affected by isotopic fractionation in firn, due to a thinner 10 diffusive column height and due to a stronger advection/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., 2011) and this leads to a small difference in δ 13 C(CH 4 ) (∼0.2 ‰) in the calculated isotopic fractionation at 80 m depth ( Fig. 4: yellow full and dashed lines). This indicates that for high accumulation sites, changing the 15 diffusivity has a limited effect on isotopic fractionation. A very small firn fractionation estimate is obtained at DML whereas 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 singlesite 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. 20 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 single-site reconstructions 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). As expected, the simulations with this modified diffusivity profile lead to a considerably stronger isotopic fractionation along the DML firn column (Fig. 4: black dashed-line).
On the other hand, the 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 5 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, turns into an advantage.
In summary, these results show that model estimates of isotopic fractionation in firn are strongly dependent on the diffusivity profile used. Although constrained in our case 10 with multiple gas histories, the firn diffusivity is not always sufficiently constrained for an accurate estimate isotopic 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 15 single sites leads to sometime inconsistent solutions. We thus try in the following to extract a "best estimate" 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. 3). In the SH multi-site in-20 version, the SPO-95 diffusivity profile is used for the SPO-01 site. The deepest firn data points provide constraints further back in time compared with 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). As the correction for isotopic fractionation in firn is most uncertain for the deepest samples, where 25 strong differences between individual firn air models have been reported (Buizert et al., 2011), the deep firn data points showing strong deviations with the modeled isotopic 9599 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | fractionation were excluded (Fig. 3). In the NH, as NM-09 diffusivity is constrained with only two gases (CO 2 and CH 4 ) instead of nine gases for NEEM-EU, a larger number of data are eliminated for NEEM-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  S4). This restricts the duration 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 multisite inversion, each firn profile is given equal weight (red lines in the left panels of Fig. 5). In a second inversion (green lines in left panels of Fig. 5), 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 supposedly affected by larger experimental uncertainties. It leads to slightly flatter δ 13 C(CH 4 ) scenarios in 15 the latter multi-site inversion, but the results in Fig. 5 show that the differences are well within the uncertainty envelopes of the original multisite inversion. The middle panels in Fig. 5 show how the modeled firn air δ 13 C(CH 4 ) profiles, based on the multisite 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 0.05-0.3 ‰. The multisite δ 13 C(CH 4 ) reconstructions (Fig. 5) show an increase in δ 13 C(CH 4 ) over the last decades, but the magnitude of the increase (0.7 ‰ over 50 yr) is about 50-60 % smaller than reported previously (Francey et al., 1999;Bräunlich et al., 2001;Sowers et al., 2005;Ferretti et al., 2005).
As a further check of consistency of the individual firn air measurements with the multisite δ 13 C(CH 4 ) scenario, we convert the measurements from each depth into temporal isotope values, using the mean age of the samples and taking into account the isotopic fractionation in firn. This approach was also used in e.g. Francey et al. (1999) and Ferretti et al. (2005). Figure 6 shows that the data from the individual sites fall within the uncertainty range of the multisite reconstruction. Note that for the readability of the figure, the age spread is not represented in Fig. 6, but note that no discrete age exist 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 multisite reconstruction, whereas the data from VOS 5 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.
The multisite reconstruction is constrained by δ 13 C(CH 4 ) data from the ten sites presented in Table 1 (i.e. excluding DI). For the NH sites, the "best estimate" scenario is constrained by only a few data points before 1990. Furthermore, the firn profiles at the three sites that enter the multisite 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 much less constrained than the SH scenario, where data from sites with very different firn 15 characteristics and sampling dates are available. The lack of independent constraints and the relatively small number of data points may explain the higher variability of the NH δ 13 C(CH 4 ) reconstruction.
Although the multisite scenarios for the SH and NH agree within the uncertainty envelopes, the shape of the best estimate slightly different scenarios for the two hemi-20 spheres. 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 interhemispheric mixing time of 1 yr, large variations on short timescales are unlikely. In any case, given the larger errors of the multisite reconstructions it is not possible to extract robust information about the IPG. 25 Alternatively, in order to combine all available independent constraints, in the last step we attempt to reconstruct a combined global δ 13 C(CH 4 ) history, assuming that the IPG of δ 13 C(CH 4 ) has not changed over the period of the reconstruction. To facilitate this, the NH data were shifted by the magnitude of an assumed constant IPG, and then 9601 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a multisite inversion was performed using all data from both hemispheres together. This exercise was performed with values for the IPG between 0.2 and 0.5 ‰. The scenario with the lowest RMS difference is obtained for an IPG of 0.3 ‰. The results from this "global multisite" inversion are depicted in green dashed-lines in Fig. 7 and they may be compared to the hemispheric multisite reconstructions (red lines).

6 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. For the NH, no direct atmospheric data are available before 2000. For the SH, δ 13 C(CH 4 ) 10 was determined from archived air samples collected at Cape Grim, Australia from the late 1970's (Francey et al., 1999). Moreover, ice core data from Law Dome, Antarctica (Ferretti et al., 2005) overlap in time with the period of our multisite reconstruction. For the NH, no δ 13 C(CH 4 ) ice core data are available for the reconstructed period, because the age of the youngest air from Greenland ice cores is older than 50 yr. We 15 thus show in Fig. 7 a comparison with five data points from EUROCORE (Sapart et al., 2012) covering the period 1800-1950. It shows that the δ 13 C(CH 4 ) reconstruction from the firn gas samples agrees well with the direct measurements over the last decade (light grey lines), which provides evidence for a negligible δ 13 C(CH 4 ) trend over this period. The Cape Grim archive data from the SH (dashed grey line in Fig. 7, lower 20 panel) show a steeper temporal trend and fall outside the low end of our reconstruction before ∼1975. Here, the fact that inter-laboratory calibration offsets can considerably hamper direct comparison of the datasets has to be taken into account. The Cape Grim and Law Dome data were measured in the same lab, and have been reported to be on the same scale as the PSU measurements (Mischler et al., 2009). An offset 25 of 0.28 ‰ between the IMAU and PSU scales has been identified in an intercalibration exercise (Sapart et al., 2011). If the Cape Grim data are corrected upwards for this difference based on this indirect evidence, they become in acceptable agreement with the lower envelope of our best estimate scenario and the direct atmospheric data. Whereas an absolute offset can thus tentatively be attributed to uncertainties in the calibration scales, the Cape Grim data also show a somewhat stronger δ 13 C(CH 4 ) trend than the multisite firn air scenario. One should keep in mind that our reconstructed 5 scenarios are sensitive to 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 leads 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) 10 from Law Dome, which are also included in Fig. 7, lower panel. A shift of 0.28 ‰ would bring these results closer to the scenario range reconstructed from the firn air data, but the Law Dome data still suggest stronger isotope changes in the past than our multisite reconstruction (this was also inferred from the single site scenario reconstruction based on DE08 firn data only shown in Fig. 2). Ice core data are also affected 15 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 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:   20 red stars). The corrected firn air data with the CSIRO model (orange dots) and with the LGEE-GIPSA model with modified diffusivity (red stars) show quite similar results which implies 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 fractionation-corrected δ 13 C(CH 4 ) shown in Fig. 7 demonstrate 25 that the 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. 7, upper panel). However, here the 9603 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | δ 13 C(CH 4 ) data corrected for isotopic fractionation are rather at the high end of the firn air reconstruction. The ice core EUROCORE data are possibly enriched 0.25 ‰ with respect to the firn air data based on comparison of ambient and ice air samples (Sapart et al., 2011), and a corresponding correction would obviously bring the data in better agreement with the firn reconstruction. The firn/ice data comparison for the Northern 5 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.

Discussion
The data evaluation and interpretation above highlights the importance of isotopic frac-10 tionation in polar firn for reconstructions of δ 13 C(CH 4 ). This is because isotopic fractionation effects are of the same order of magnitude or even larger than the atmospheric δ 13 C(CH 4 ) signal. When analyzed with the same inverse firn air model, single site reconstructions from 11 sites lead to strongly different atmospheric isotope histories. Investigation of the data used to derive diffusivity profiles indicate that in general diffu-15 sivity profiles constrained by a larger number of atmospheric species lead to stronger isotopic fractionation effects and thus smaller isotope trends in the model than at sites with diffusivity profiles constrained by less species. For the two firn air sampling campaigns performed at South Pole, our work shows that a single site inversion of the SPO-01 data agrees much better with reconstructions from other sites when the pre- 20 sumably better constrained diffusivity profile of SPO-95 is used. When the deepest firn data-points are neglected, due to assumed larger uncertainties in the isotopic fractionation affecting these samples, it is possible to build a multisite firn air δ 13 C(CH 4 ) reconstruction that leads to consistent results at all sites (when the diffusivity profile of SPO-95 is used for the SPO-01 data). This scenario is consistent with recent at-core data. Those differences maybe due to calibration issues, uncertainties on the CH 4 scenario and/or insufficiently constrained diffusivity. Another source of uncertainty is the firn model itself. Buizert et al., (2011) 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 fraction-5 ation in the lock-in-zone. This may bias our reconstructions slightly towards smaller required atmospheric trends. On the other hand, the correction for isotopic fractionation along the Law Dome firn column and in ice bubble show that our model is consistent with the CSIRO model and yields the same trend in δ 13 C(CH 4 ). Our best estimate scenario based on multi-site inversion can be considered as the state-of-the-art of con-10 clusions bearable from the existing δ 13 C(CH 4 ) firn air data. However, due to uncertainties in estimating the isotopic fractionation effects, it remains tentative to reconstruct δ 13 C(CH 4 ) during periods of strong CH 4 increase. The difficulties in constraining the atmospheric isotope history of δ 13 C(CH 4 ) may turn into an advantage when investigating the consistency between diffusivity profiles at dif-15 ferent firn air sites. The sensitivity of δ 13 C(CH 4 ) to diffusivity profiles provide a useful model constraint if the atmospheric record is independently known. But for the time being, this turns into a circular argument since the δ 13 C(CH 4 ) history since pre-industrial times is not uniquely defined yet. Our reconstructed δ 13 C(CH 4 ) trend between 1960 and 2009 is only half of the pre-20 viously published trends in δ 13 C(CH 4 ) (Francey et al., 1999;Bräunlich et al., 2001;Sowers et al., 2005). On the other hand, the measurements on ice core samples in the period before the strong increase in CH 4 started clearly implies that there must have been a strong increase in δ 13 C(CH 4 ) since pre-industrial times. Whereas the previous studies concluded that most of the enrichment occurred during the last 50 yr, the firn air 25 reconstruction would imply that there was also a strong increase already before 1960.

9605
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Conclusions
We reconstructed δ 13 C(CH 4 ) atmospheric scenarios over the last 50 yr for both Northern and Southern hemispheres using a firn transport model taking as input δ 13 C(CH 4 ) firn air measurements at 11 sites. Our calculations reveal discrepancies between sites while running single-site inversion and show that the estimate of isotopic fractionation 5 between 12 CH 4 and 13 CH 4 along the firn column is highly dependent on how good the diffusivity profiles are constrained for 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 a best estimate of δ 13 C(CH 4 ) for each hemisphere from a multisite- 10 inversion, which fits all measurements within the errors of the reconstruction. This scenario shows an increase of about 0.7 ‰ between 1960 and 2000 followed by a period of relatively stable δ 13 C(CH 4 ) values since then, in agreement with direct atmospheric measurements. The scenario shows, however, remaining discrepancies with the deepest firn air data, and with published ice core data (Ferretti et al., 2005) for the same 15 period. Interlaboratory calibration offsets likely contribute to the differences with the ice core data, but uncertainties in the isotopic fractionation in firn air models is important as well. The reconstructed isotope trend between 1960 and 2009 is a factor of 2 smaller than previously published trends in δ 13 C(CH 4 ) (Francey et al., 1999;Bräunlich et al., 2001;Sowers et al., 2005). Aside from NEEM, firn air sampling programs were supported by: at DI, IPEV (France) and the Canadian Continental Polar Shelf Project; at BI, IPEV and the British Antarctic Survey; at VOS, IPEV; at SP: US-NSF; at DC, IPEV, ENEA (Italy) and the ESF/EC program EPICA; at 10 DML: British Antarctic Survey; at DE08; at NGRIP: Funding Agencies in Denmark (SNF), Belgium (FNRSCFB), France (IPEV and INSU/CNRS), Germany (AWI), Iceland (RannIs), Japan (MEXT), Sweden (SPRS), Switzerland (SNF) and the United States of America (NSF). We thank all participants to the field work at the 11 sites for drilling and firn air sampling operations.

9617
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 grey curve (left panel) represents air--archive from Cape Grim (Francey et al., 1998) and 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 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. Figure 7: Multi--site "best estimate" δ 13 C(CH4) scenario compare with atmospheric archive from Cape Grim (dashed grey line) and ice core data (dots). The red dashed--lines are the "best--estimate" scenario for each hemisphere. The green dashed--lines are the scenarios including all sites from both hemispheres together and using an inter--polar gradient of 0.3‰. a: NH raw δ 13 C(CH4) data (light blue dots) and δ 13 C(CH4) data corrected for firn fractionation from EUROCORE ice core (Sapart et al., under review). 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 "best estimate" δ 13 C(CH 4 ) scenario compare with atmospheric archive from Cape Grim (dashed grey line) and ice core data (dots). The red dashed-lines are the "bestestimate" scenario for each hemisphere. The green dashed-lines are the scenarios including all sites from both hemispheres together and using an inter-polar gradient of 0.3 ‰. (a) NH raw δ 13 C(CH 4 ) data (light blue dots) and δ 13 C(CH 4 ) data corrected for firn fractionation from 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).