Towards reconstructing the Arctic atmospheric methane history over the 20th century: measurement and modeling results for the NGRIP firn

Systematic measurements of atmospheric methane (CH4) mole fractions at the northern high latitudes only began in the early 1980s, and whilst CH4 measurements from Greenland ice cores covered the period before ~1900, no reliable 15 observational record is available for the intermediate period. In this study, we reconstruct the atmospheric CH4 for that period, when the mole fraction started to increase rapidly. We use a set of trace gas data measured from firn (an intermediate stage between snow and glacial ice formation) air samples collected at the NGRIP (North Greenland Ice Core Project) site in 2001, in combination with a firn air transport model whose performance is validated by using a set of published firn air data at the NEEM (North Greenland Eemian ice Drilling) site. We examine a variety of possible firn diffusivity profiles using a suite of 20 measured trace gases, and reconstruct the CH4 mole fraction by an iterative dating method. We find that, given the currently available firn air data sets from Greenland, reliable reconstruction of the Arctic CH4 mole fraction is possible only back to the mid 1970s. For the earlier period, it is difficult to identify the atmospheric CH4 history that consistently reproduce the depth profiles of CH4 in firn at both NGRIP and NEEM sites. Therefore, the currently proposed Arctic CH4 history should still be considered preliminary and uncertain, and should not be treated as the known history for constraining firn-air transport models. 25

A large fraction of natural and anthropogenic CH4 sources resides in the northern hemisphere, and thus the atmospheric CH4 trend of the northern hemisphere can provide important information on the evolution of anthropogenic CH4 emissions as well as the variations of natural CH4 emission in response to climatic variability. The interhemispheric gradient of CH4 mole fraction is also key to the allocation of CH4 emissions between both hemispheres (e.g. Dlugokencky et al., 2003;Chandra et al., 2021). 40 Systematic measurements of atmospheric CH4 mole fraction began in the 1980s. In Figure 1, CH4 measurements at two Arctic sites Barrow, Alaska (BRW) and Alert, Canada (ALT) are shown (orange dots), whose records start in 1983 and 1985, respectively (data provided by National Oceanic and Atmospheric Administration/Earth System Research Laboratory/Global Monitoring Laboratory, NOAA/ESRL/GML). Although some sparse data from the 1970s are available (e.g. Rice et al., 2016), such "direct" measurements provide CH4 data only since around 1980, which means that some reconstruction methodology is 55 required to infer atmospheric CH4 mole fraction variations before that. For this purpose, air extracted from ice cores and firn layers has been measured. Figure 1 also presents CH4 mole fractions analyzed in ice cores from Greenland (open blue symbols), such as Eurocore (Blunier et al., 1993), Eurocore and GISP2 (Greenland Ice Sheet Project 2) (Etheridge et al., 1998) and Site J (Nakazawa et al., 1993). These data show fairly good agreements with each other until ~1900, after which the number of data and the consistency among the records are poor. Continuous measurements from NEEM-S1 ice core (Rhodes et al., 2013) 60 presented CH4 mole fractions before ~1960 (red in Figure 1), but their data are notably higher than the ice core data after ~1850. Therefore, the data gap between ~1900 and ~1980 is not covered reliably either by direct observations or ice core reconstructions.
In comparison to these data, higher-resolution CH4 data are available from Antarctica (grey in Figure 1). The comprehensive Law Dome ice core and firn data set (Etheridge et al., 1998;MacFarling Meure et al., 2006) almost continuously covers the 65 last 200 years and are well connected to the direct measurements at South Pole (SPO, data provided by NOAA/ESRL/GML) and Cape Grim (CGO) (Etheridge et al., 1998;MacFarling Meure et al., 2006). Such consistency and continuity among the datasets suggests that the Antarctic CH4 data can serve as a good reference to represent the global atmospheric CH4 trend over the past centuries. Buizert et al., (2012) used the Antarctic CH4 trend with presumed interpolar difference (IPD) to propose the CH4 history over Greenland (light blue in Figure 1), which in turn was used to constrain the gas diffusivity profile in firn at 70 the North Greenland Eemian ice Drilling (NEEM) site. In their study, the Arctic CH4 scenario was constructed by adding IPD to the Law Dome data, where the IPD was assumed to be proportionally correlated with the growth rate of CH4. This seems reasonable assumption, given that the IPD and growth rate are both largely subject to changes in emissions from the northern hemisphere (Dlugokencky et al., 2003;Chandra et al., 2021). The Arctic CH4 scenario by Buizert et al. (2012) is consistent with the direct measurements that started in the 1980s, and ice core reconstructions before ~1900 (Figure 1). It is thus 75 considered as the most likely synthetic atmospheric CH4 trend for the northern high latitudes, and is treated as a "known" history, by which the diffusivity profiles in firn are tuned in firn-air transport models Trudinger et al., 2013). It is however noted that this scenario has not been validated against independent estimates for the data gap period (from ~1900 to ~1980).
In this study, we present a set of mole fractions of CH4 and other trace gases in firn air samples collected at the North Greenland 80 Ice Core Project (NGRIP) site. Using the atmospheric scenarios by Buizert et al. (2012), we simulate the depth profiles of trace gases in the NGRIP firn with our firn-air transport model. We examine a variety of modeling cases for different diffusivity profiles and reconstruct the Arctic atmospheric CH4 over the late 20th century using the iterative dating approach (Trudinger https://doi.org/10.5194/acp-2021-736 Preprint. et al., 2002). The reconstructed CH4 trend is evaluated by comparison to the original atmospheric CH4 scenario and to the NEEM firn air data. Uncertainty of the Arctic atmospheric CH4 history for use in firn air modeling is discussed. 85

Experimental method
Firn air was sampled at the Greenland site NGRIP (75°10′N, 42°32′W, 2959 m AMSL) in May-June 2001. Accumulation, surface density, mean temperature and pressure are 179 kg m −2 yr −1 , 300 kg m −3 , 241 K and 680 hPa, respectively. Details of the firn and firn-air sampling have been described elsewhere (Kawamura et al., 2006;Ishijima et al., 2007). At the NGRIP site, two shallow holes (EU and Japanese holes) were drilled (Kawamura et al., 2006;Landais et al., 2006), and the present data 90 are from the firn air samples collected from the Japanese hole. The total number of air-sampling depths is 24.
Since the technical details are reported in Kawamura et al. (2021), only brief descriptions of relevant data presented in this study are given here. CH4 mole fractions of the firn air samples were measured using a gas chromatograph (Agilent 6890, Agilent Technologies Inc.) equipped with a flame ionization detector (GC-FID) at Tohoku University (TU), with a reproducibility of 2 ppb (Umezawa et al., 2014). The CH4 mole fractions were determined against our working standard gases 95 that were calibrated on the TU1987 CH4 scale (Aoki et al., 1992;Umezawa et al., 2014;Fujita et al., 2018). The difference between the TU1987 CH4 scale and the WMO CH4 mole fraction scale (on which the NEEM CH4 data were measured) is estimated to be ~0.5 ppb at the current atmospheric CH4 levels (Fujita et al., 2018). Oyabu et al. (2020) reported that ice core data analyzed on the TU1987 and WMO scales showed good agreement within analytical uncertainties, indicating consistency of both scales, including for the lower mole fractions (e.g. ~700 ppb). It is therefore likely that the difference between both 100 scales is well below the variations of interest in this study, and thus no correction is applied for use of the NGRIP and NEEM firn data.
The firn air samples were measured for CO2 and SF6 mole fractions respectively by using a nondispersive infrared gas analyzer (NDIR) and a gas chromatograph equipped with an electron capture detector (GC-ECD) at TU. The measurement reproducibility is estimated to be 0.02 ppm for CO2 and 0.09 ppt for SF6 and mole fractions of both gases are reported on the 105 TU2010 CO2 and TU2002 SF6 scales, respectively (Sugawara et al., 2018).
The NGRIP firn air samples were also analyzed for selected halocarbons (CFC-11, CFC-12, CFC-113 and CH3CCl3) on the Vacuum Preconcentration and Refocusing-Gas Chromatography-Mass Spectrometry (VPR-GCMS) system, which was developed based on the work by Saito et al. (2006). An aliquot of the sample was transferred into an evacuated canister of ~0.3 L at around ambient pressure (~100 kPa) and the inner pressure of the canister was recorded. The air is extracted by a vacuum 110 pump through a preconcentration trap filled with HayeSep D cooled to -135° C using a Stirling cooler. The preconcentration trap was heated to -70°C to release major atmospheric constituents and then up to 100°C to transfer the trapped compounds to a cryofocusing trap containing Carboxene 1000/Tenax TA at −100°C. The trap was then heated to 180°C to inject the are determined against a working standard gas (compressed dry air) that was calibrated against synthetic standards on the 115 NIES-08 scales.

Firn air transport model
Since the gas diffusivity in firn layers is significantly lower than in the atmosphere, the movements of atmospheric constituents are driven mostly by molecular diffusion according to their vertical mole fraction gradients under the influence of gravity. In general, lighter air components (or isotopologues) diffuse faster under their mole fraction gradients, while heavier components 120 accumulate in the deeper layers due to the gravitational effect. Hence the depth profiles of trace gas mole fractions in firn are determined by the atmospheric histories transferred towards depth in the firn by the molecular diffusion driven by the mole fraction gradient and gravity. At the bottom of firn, the air is trapped as bubbles in the ice sheet, which creates slow downward motion of firn air.
The firn column can be divided into three zones: a convective zone (CZ), a diffusive zone (DZ) and a lock-in zone (LIZ) 125 (Sowers et al., 1992;Kawamura et al., 2006;Buizert et al., 2012). In CZ, primarily driven by surface winds and fluctuations of atmospheric pressure, air is mixed with the overlying atmosphere (Sowers et al., 1992;Kawamura et al., 2006). The CZ thickness at NGRIP is estimated to be below 2 m (Kawamura et al., 2006). In DZ, which is sufficiently isolated from the surface turbulence, movement of air is governed by molecular diffusion. Gravitational enrichment according to the barometric equation (i.e. linear increases of δ 15 N of N2 and δ 18 O of O2) occurs with depth, which stops at the top of LIZ (Sowers et al., 130 1992;Schwander et al., 1993;Kawamura et al., 2006). The top of LIZ (lock-in depth) at NGRIP is at depth 63 m (Kawamura et al., 2006). In LIZ, advection with the enclosing ice matrix dominates the transport of air, and air parcels are gradually isolated as bubbles. Traditionally, it was supposed that high-density impermeable ice layers stop diffusivity in LIZ completely, however, recent studies demonstrated finite diffusivity in LIZ (Severinghaus et al., 2010;Buizert et al., 2012;Trudinger et al., 2013). The deepest air sampling at NGRIP was successfully made at 77.71 m, and total pore closure is considered in the deeper 135 layers in our modeling.

Modeling firn air transport
We use a one-dimensional diffusion model that has been used for the reconstruction of isotope ratios of CO2 and N2O (Sugawara et al., 2003;Ishijima et al., 2007). The model is conceptually similar to that developed by Trudinger et al. (1997); it is based on a theoretical formation of diffusion ) and a bubble trapping process (Rommelaere et al., 140 1997). Air movement in the firn is driven by molecular diffusion and a gravitational effect. Namely, a trace gas flux in firn is expressed by where D is the effective diffusivity of a trace gas molecule, the variables s, c, and T are open porosity, trace gas molar concentration, and firn temperature, respectively, and the constants m, g, and R are the mass number of the trace gas, the 145 acceleration of gravity, and the gas constant, respectively. Vertical advection flux of the trace gas, caused by air trapping at the close-off zone and downward bulk motion of firn, is expressed by using the equation given by Rommelaere et al. (1997).

Conservation of the trace gas is given by
for the open pore space and 150 zone where the open pore air is gradually trapped into bubbles, mass conservation is given by using a bubble trapping rate r 155 (s −1 ), which simply means that a portion of the trace gas molar concentration in the open pore space (rc) is added to bubbles.
The bubble trapping rate is given as a function of the open porosity, the total porosity, and the vertical speed of firn itself (Rommelaere et al., 1997). The total porosity was calculated from the firn density data. At the transition zone, the total porosity should be divided into the open and closed porosity. The closed porosity sc was calculated by the empirical equation given by Schawander (1989). 160

Effective diffusivity
We follow the previous firn-air studies in which the effective diffusivity in firn is optimized with an iterative method so as to minimize the difference between the simulated and observed depth profiles of CO2 (Sugawara et al., 2003;Ishijima et al., 2007). An initial guess of the effective diffusivity for CO2, Dinit(z), was calculated by: where s(z) and D0 represent the open porosity at a depth z and the diffusion coefficient of CO2 at 253 K and 1013 hPa, respectively. D0 was set to 1.247×10 −5 (m 2 s −1 ) according to Trudinger et al. (1997). The bulk density was determined by measuring the dimension and weight of cylindrically cut firn core samples (Kawamura et al., 2006). The effective diffusivity of CO2 thus obtained was converted to those of other trace gases by multiplying by scaling factors from Buizert et al. (2012). Therefore, the depth profile pattern of the effective diffusivity is identical among all gases, but the magnitude is gas-dependent 170 due to the scaling factors. In this study, the effective diffusivity profile prepared for the NGRIP firn by Ishijima et al. (2007) is modified to improve the reproducibility of our newly measured trace gas profiles (section 4).

Performance of the firn air transport model
To validate our firn air transport model, we simulated depth profiles of various trace gases in the NEEM firn. Our model did not participate in the model intercomparison study using the NEEM data (Buizert et al., 2012). In this simulation, we employed 175 the atmospheric scenarios prepared by Buizert et al. (2012) as per their model intercomparison. The diffusivity profile optimized for the CIC (Center for Ice and Climate, Niels Bohr Institute, University of Copenhagen) model was tuned for our model. The simulated depth profiles of various compounds are presented in Figure 2 and compared with those by other models presented in Buizert et al. (2012). The results confirm that performance of our model is comparable to those by other groups.
As a measure of the model performance, Buizert et al. (2012) compared root mean square deviation (RMSD) between the 180 model and data, which ranged from 0.73 to 0.92 for the six participating models. Following the same approach, our model yields the RMSD value of 0.83 of for the NEEM EU borehole.

Modeling procedure
Our modeling procedure for reconstructing the atmospheric histories of trace gases was as follows: (1) To represent atmospheric trace gas trends in the Arctic region, we began by employing atmospheric scenarios prepared for the NEEM firn air modeling (Buizert et al., 2012) for all the trace gases presented in this study (CH4,CO2,SF6,. The firn transport model calculates depth profiles of these trace gases in the NGRIP firn with the initial effective diffusivity described in section 3.2. (2) We examined the different sets of diffusivity profiles so as to improve overall reproducibility of the modeled profiles of the six trace gases except CH4. In other words, we regarded these six trace gases as constraints to the diffusivity profile. Note that this step assumes that the atmospheric histories of the six gases are known with sufficient accuracy to infer a range of 200 acceptable diffusivity profiles that reproduce the depth profiles of the firn air composition. We prepared 100 different sets of modified diffusivity profiles, and the overall model-data agreement for each set was evaluated based on the RMSD of the six gases.
(3) We ran the model with the 100 different sets of diffusivity profiles to calculate the depth profile of CH4. Note that, based on the earlier steps, we know the diffusivity profiles that generate reasonable firn-air profiles for the six trace gases. Every 205 diffusivity profile was used in combination with the firn-air CH4 data for reconstructing an atmospheric CH4 history. We employed an iterative dating approach (Trudinger et al., 2002) where the initial atmospheric scenario was modified to improve model reproducibility of the CH4 depth profile (see below). The corrected atmospheric CH4 scenarios were then compared to the initial scenario for further discussion.
The iterative dating for CH4 was performed as follows: 210 (1) Depth profile of CH4 was calculated with the initial atmospheric CH4 scenario.
(2) Effective age at each sampling depth was calculated so that the modeled CH4 mole fraction, at that depth, agreed with a value in the atmospheric CH4 scenario.
(3) A new atmospheric CH4 scenario was constructed by assigning the observed CH4 mole fraction, at each depth, to the effective age. 215 (4) Depth profile of CH4 was again calculated with the revised atmospheric CH4 scenario.
(5) The above steps 2-4 were repeated until the model-data difference converged within an acceptable range (typically after a few iterations) (Trudinger et al., 2002;Ishijima et al., 2007). In this study, we made five iterations for each modified diffusivity case as we confirmed sufficient convergence of the result.  Figure 3 presents the initial simulation results in comparison to the observed profiles for the six trace gases (CO2, SF6, CFC-11, CFC-12, CFC-113 and CH3CCl3). As seen in this figure, measured profiles of these trace gases (except CH3CCl3) show gradual decreases with depth in the DZ and sharp decreases in the LIZ. In contrast, CH3CCl3 increases with depth in the DZ and sharply decreases in the LIZ. The difference of the depth profile pattern among species is due to their different historical 225 atmospheric trends. It is known that the atmospheric mole fractions of the five trace gases (CO2, SF6, CFC-11, CFC-12 and CFC-113) have increased monotonically since the mid 20th century (Sturrock et al., 2002;Martinerie et al., 2009). In contrast, CH3CCl3 has increased until the early 1990s and has decreased since then (Sturrock et al., 2002;Rigby et al., 2017). Our simulation reproduces the observed depth profiles of these six trace gases in the NGRIP firn fairly accurately with the atmospheric scenarios prepared for modeling the NEEM firn air (Buizert et al., 2012). 230 It is interesting to note that the depth profiles of CH3CCl3 at NGRIP and NEEM are remarkably different. Whereas the NEEM data show a relatively sharp CH3CCl3 peak in the LIZ (~65 m, Figure 2), NGRIP does not show such a narrow peak. This is due to the timing of firn air sampling i.e. 2001 for NGRIP and 2008 for NEEM. When the NGRIP firn air was sampled, the 240 signal of the maximum atmospheric CH3CCl3 in the early 1990s had only reached near the top of the LIZ at the site, thereby formulating the relatively gentle changes at the shallower depths. On the other hand, seven years later at the NEEM site, such signal was found deeper in the LIZ where the age of air changes rapidly with depth in both deeper and shallower sides. We emphasize that, despite the differences in the depth profiles of CH3CCl3 at the two sites, our simulations reproduce both profiles measured at both sites well, using the same atmospheric CH3CCl3 scenario. 245 Figure 3 shows that the model-data difference increases in the LIZ for all the trace gases. In particular, the model-data difference is pronounced as a dip around 70 m for all trace gases, implying that the mismatches may originate in a common factor in the modeling e.g. depth profile of diffusivity. To examine the impact of diffusivity modification on the simulated depth profiles and their agreements with the data, we prepared 100 sets of modified diffusivity profiles (Figure 4). The various diffusivity profiles were constructed by modifying the original profile at certain range of depths in a stepwise manner. Although 250 this simple method does not guarantee identification of a best-match profile, we are confident that an acceptable range of the diffusivity profile is satisfactorily constrained. To reduce the model-data difference in the LIZ (Figure 3), we found that the diffusivity needs to be increased in the shallower layers than the top LIZ i.e. 50-65 m. The diffusivity was also modified in the deeper layers (>65 m) and simulations were made accordingly. The simulated profiles for depths deeper than 50 m using the 100 diffusivity profiles are presented for the six trace gases 260 ( Figure 5). In this figure, the modeling results with the modified diffusivity profiles are shown in colors on the left axis, and the model-data differences of the respective cases are plotted on the right axis. It is clear that the model-data differences could be significantly reduced with some diffusivity cases. We calculated RMSD in the same manner as Buizert et al. (2012) with the measurement uncertainties of 0.2 ppm for CO2, 0.2 ppt for SF6, 1.1 ppt for CFC-11, 3.3 ppt for CFC-12, 0.6 ppt for CFC-113 and 3.2 ppt for CH3CCl3. The RMSD values are as little as 0.51 for a particular case. In Figure 4 and   Likewise, depth profiles of CH4 are shown in Figure 6. These calculations were made with the 100 diffusivity profiles and the atmospheric CH4 scenario for the NEEM firn-air modeling (Buizert et al., 2012). It is interesting to note that the characteristics 275 of the model-data difference for CH4 are different from those for the other six trace gases. For the other trace gases (Figure 5), the initial simulation showed increased model-data differences around 70 m, and they were reduced with some modified diffusivity profiles. For CH4, the model run with the original diffusivity profile reproduces the observed CH4 profile quite well down to ~73 m, but significantly overestimates by >20 ppb for the lowest three depths (the black solid line). Using the modified diffusivity profiles that allowed better agreements for the other six trace gases (orange and red lines), we actually find larger 280 model-data CH4 difference than in the initial simulation.

Reconstruction by iterative dating
We explore the possibility of reconstructing the Arctic CH4 mole fractions by the iterative dating method. This approach sets the diffusivity profile fixed (assumed to be correct) and aims to find an acceptable atmospheric history to reproduce the firnair depth profile. It is considered that the modified diffusivity profiles with the low RMSD values for the six trace gases are adequately evaluated and that the model-data mismatch in the CH4 modeling is therefore attributable to uncertainty in the 290 atmospheric CH4 scenario. The atmospheric CH4 scenarios obtained by the iterative dating method are presented in Figure 7 for the 100 modeling cases of the modified diffusivity. The different simulation cases are colored in the same manner as in the earlier figures (Figures 4-6). Note that the reconstruction cases colored in light blue are considered to be less likely, due to poorer reproduction of the depth profiles (Figures 5 and 6). The model results in red and orange are in good agreement with the initial atmospheric scenario by Buizert et al. (2012) after around 1980. For the earlier period, however, the upper bounds 295 of the reconstructions matches with the initial scenario, and the overall range of acceptable histories is below the initial atmospheric CH4 scenario (circles and shades in red and orange), suggesting that the CH4 mole fractions may have been lower than the initial modeling scenario.

305
It should be stressed that the reconstructions for the time period before 1980 rely on the CH4 data from the lowest three depths of the NGRIP firn air. The differences between the initial and corrected atmospheric CH4 scenario from these three deepest data are up to ~100 ppb. As seen in Figure 7, such reduction in the Arctic atmospheric CH4 scenario over the period would result in alignment with the atmospheric CH4 trend in Antarctica inferred from the Law Dome and other datasets (black dotted 310 line). It is again noted that all the reconstruction cases colored in orange and red used diffusivity profiles that yield relatively good model reproducibility with RMSD values of <1.0 and 0.75, respectively. It is therefore seen that iterative dating-based reconstruction from the NGRIP firn data suggests decreased CH4 mole fraction from the 1950s to 1970s in any case, albeit with large uncertainty.

Comparison to the NEEM profile
To further evaluate consistency of the atmospheric CH4 scenario reconstructed from the NGRIP firn, we ran the model for the NEEM firn using the corrected atmospheric CH4 histories obtained from the iterative dating approach. The modeling results with the initial and modified atmospheric CH4 histories are shown in Figure 8. With the initial scenario (black line), the modeldata difference was <25 ppb, with relatively large differences appearing in the LIZ (>63 m). However, when the atmospheric CH4 histories reconstructed from the NGRIP firn were employed, the reconstruction cases based on the improved reproducibility for the NGRIP firn (orange and red lines) resulted in larger model-data differences of >50 ppb in the LIZ.

Discussion and conclusion 330
In section 4, we have shown that, with the atmospheric scenarios prepared for the NEEM firn, depth profiles of CO2, SF6, CFC-11, CFC-12, CFC-113 and CH3CCl3 in the NGRIP firn are accurately reproduced by using the modified diffusivity profiles ( Figure 5). This suggests that the atmospheric scenarios of these trace gases used in the modeling are consistent with the depth profiles observed at both firn sites in Greenland (NEEM and NGRIP). In contrast, the observed CH4 profile in the NGRIP firn was not reconciled using the atmospheric scenario prepared for the NEEM firn and the diffusivity profiles that 335 allow adequate reproducibility for the above six trace gases. This suggests either that the Arctic atmospheric scenario of CH4 is uncertain or that the diffusivity profile of the NGRIP firn is underconstrained.
We explored the correction of the atmospheric scenario of CH4 by the iterative dating approach. This method improves agreement to the observed CH4 depth profile in the NGRIP firn, with an implicit assumption that the diffusivity profile in each case is correct. Although uncertainty due to the under-constrained diffusivity profile in the LIZ is large, this attempt suggested that the CH4 mole fractions over the period 1950-1980 could be decreased in comparison to the original scenario prepared for the NEEM modeling (Figure 7). However, the follow-up modeling with the corrected atmospheric CH4 histories by the iterative dating based on the NGRIP data did not show improved reproducibility for the CH4 profile at the NEEM site (Figure 8). It is therefore implied that neither the original CH4 scenario by Buizert et al. (2012) nor the corrected CH4 histories by this study, are well validated against the available two Greenland firn data sets (NGRIP and NEEM). The decrease of up to 100 ppb from 345 the original scenario over the period 1950-1980, as suggested by the iterative dating in this study, would make the CH4 mole fraction as low as that in Antarctica (Figure 7). We point out that such a nearly-zero IPD of the CH4 mole fraction is highly unlikely, given that major fraction of both natural and anthropogenic CH4 sources resides in the northern hemisphere in both preindustrial and industrial periods (Houweling et al., 2000;Ghosh et al., 2015;Saunois et al., 2020;Chandra et al., 2021).
It is important to note that the reconstructions for the period before 1980 were heaviliy influenced by the three deepest data in 350 the LIZ (> 74 m). Figure 9 shows distributions of the effective age of CH4 at depths below 55 m in the NGRIP firn, colored the same as in the earlier figures. In addition, the spread of the effective age (σage) at each sampling depth is shown on the right axis. This figure shows that firn air samples collected at the three lowest sampling depths have effective ages corresponding to the period from ~1950 to the mid 1970s. In addition, even the acceptable diffusivities yield the spread of the effective age of >5 years at those depths (red and orange). This shows that the reconstruction of the CH4 mole fraction for the period is 355 subject to much uncertainty in effective age. The Antarctic atmospheric CH4 record (see Figures 1 and 7) indicates that the atmospheric increase rate of CH4 was 10-15 ppb yr −1 over the period. The 5-year uncertainty in the age estimate for the NGRIP firn air samples could therefore be translated to an uncertainty of >50 ppb in the Arctic atmospheric CH4 level. This is comparable to the IPD of CH4 mole fraction during the 1950s to 1970s, which was assumed when Buizert et al. (2012) prepared the Arctic atmospheric CH4 scenario for the NEEM modeling. 360 The IPD of the atmospheric CH4 mole fraction is important for better understanding the evolution of the global CH4 budget.
Given that the Antarctic ice core and firn measurements have provided relatively reliable CH4 records over the 20th century, 370 improved reconstructions from Greenland ice cores and firn air should better constrain the changes in the IPD. To sufficiently constrain the historical global CH4 budget, the reconstruction for Greenland needs to be accurate within ~10 ppb, corresponding to ~30 Tg CH4 yr −1 global emission. Unfortunately, based on the firn CH4 data from NGRIP and NEEM, this study demonstrated that consistent reconstruction of the Arctic CH4 mole fraction is achievable only back to the mid 1970s, and that the uncertainty of reconstruction is still very large (around 50 ppb) for the 1950s to 1970s. 375 Previous studies have used CH4 as one of tracers that constrain diffusivity profiles in firn Trudinger et al., 2013). These studies have shown that CH4 effectively contributed to constraining diffusivity for deep layers i.e. LIZ in the NEEM firn. However, the use of CH4 as a strong constraint is valid only when its atmospheric scenario given as input to the models was assumed to be correct. Figures 5 and 6 show a larger spread in the model-data differences of CH4 among the different diffusivity cases spread widely in the deepest layers of the LIZ, in comparison to the other six gases. This fact 380 highlights that subtle changes of the diffusivity profile in the LIZ have a large impact on simulating CH4, and it indeed indicates that the species could serve as an effective diffusivity constraint, if its atmospheric scenario was correctly given with low uncertainty. Although the currently proposed Arctic CH4 history (Buizert et al., 2012) appears a useful choice in firn modeling at Greenland sites, it should be kept in mind that the use of CH4 as a tuning tracer could lead to overfitting of the diffusivity profile. 385 Sapart et al. (2013) examined the reconstruction of stable carbon isotope ratio (δ 13 C) of atmospheric CH4 using firn air measurements from both northern and southern hemispheres. They concluded that, with the available firn measurements and understanding of firn air transport, it is difficult to reliably reconstruct the past trend of δ 13 C of CH4 because of multiple reasons including uncertainty in the atmospheric CH4 mole fraction scenario. Although there are many important and uncertain factors, the accurate reconstruction of the atmospheric CH4 mole fraction is particularly important, because the trend in the mole 390 fraction can lead to significant signal in the modeled δ 13 C profile in firn due to the difference in the molecular diffusion coefficient, even in the absence of a temporal trend in atmospheric δ 13 C. This study revealed a large uncertainty in the Arctic CH4 mole fraction trend over the 20th century, which supports the conclusion of Sapart et al. (2013) on the difficulty of reconstruction of the δ 13 C of atmospheric CH4 in the northern hemisphere. The NGRIP firn air samples were also analyzed for δ 13 C and stable hydrogen isotope ratio (δD) of CH4 (Kawamura et al., 2021), but we regrettably report that reconstruction of 395 δ 13 C of CH4 is difficult despite our best modeling efforts (not shown).
A possibility to improve the reproducibility of the depth profile of CH4 in the NGRIP firn may come from additional constraint to the diffusivity profile along with those currently made by the six trace gases (CO2, SF6, CFC-11, CFC-12, CFC-113 and CH3CCl3). It could be performed by introducing additional effective tracers e.g. 14 CO2, although such measurements cannot be foreseen for the NGRIP samples. This study showed that, with the currently available firn data from Greenland (NGRIP 400 and NEEM), reliable reconstruction of the Arctic CH4 mole fraction history (independent of Antarctic records and assumed IPD) is possible back to the mid 1970s only. For the earlier period, consistent reconstruction with a small uncertainty of age estimate is currently difficult. Future sampling and measurements of ice cores at a high-accumulation site in Greenland (where age of air occluded can be determined accurately) may be the only way to reconstruct the atmospheric CH4 trend over the 20th century. 405

Data Availability
The composition data of the NGRIP firn air samples are available on the Arctic Data archive System (ADS) of National Institute of Polar Research (https://ads.nipr.ac.jp/dataset/A20210609-001). The NEEM firn air data are available in the supplementary file of Buizert et al. (2012). Our modeling data are available upon request (umezawa.taku@nies.go.jp).

Author Contribution 410
TU, SS, KK and IO discussed on design of the study. KK, SA and TN conducted firn air sampling at the NGRIP site. SS, KK, TU and TS analyzed the firn air samples for trace gases. SJA set up the measurement system for the halocarbons. TU and SS made firn air model simulations. TU analyzed the measurement/modeling data and prepared the manuscript with contributions from all co-authors.

Completing interests 415
The authors declare that they have no conflict of interest.