Towards reconstructing the Arctic atmospheric methane history over the 20th century: measurement and modelling results for the North Greenland Ice Core Project firn

Systematic measurements of atmospheric methane (CH₄) mole fractions at the northern high latitudes only began in the early 1980s. Although CH₄ measurements from Greenland ice cores consistently covered the period before ∼ 1900, no reliable observational record is available for the intermediate period. We newly report a data set of trace gases from the air trapped in firn (an intermediate stage between snow and glacial ice formation) collected at the NGRIP (North Greenland Ice Core Project) site in 2001. We also use a set of published firn air data at the NEEM (North Greenland Eemian ice Drilling) site. The two Arctic firn air data sets are analysed with a firn air transport model, which translates historical variations to depth profiles of trace gases in firn. We examine a variety of possible firn diffusivity profiles, using a suite of measured trace gases, and reconstruct the CH₄ mole fraction by an iterative dating method. Although the reconstructions of the Arctic CH₄ mole fraction before the mid-1970s still has large uncertainties (> 30 ppb – parts per billion), we find a relatively narrow range of atmospheric CH₄ history that is consistent with both depth profiles of NGRIP and NEEM. The atmospheric CH₄ history inferred by this study is more consistent with the atmospheric CH₄ scenario prepared for the NEEM firn modelling than that for the CMIP6 (Climate Model Intercomparison Project Phase 6) experiments. Our study shows that the atmospheric CH₄ scenario used for the NEEM firn modelling is considered to be the current best choice for the Arctic CH₄ history, but it should not be used to tune firn air transport models until being verified by further measurements from sources such as the Arctic ice cores. Given the current difficulty in reconstructing the CH₄ history with low uncertainty from the firn air data sets from Greenland, future sampling and measurements of ice cores at a high-accumulation site may be the only way to accurately reconstruct the atmospheric CH₄ trend over the 20th century.

1 Introduction

Methane (CH₄) is an important atmospheric greenhouse gas emitted from both natural and anthropogenic sources. Despite making great efforts to understand its global budget, emission estimates of individual sources still have large quantitative uncertainties (e.g. Saunois et al., 2020; Chandra et al., 2021). Anthropogenic activities have enhanced CH₄ emissions globally and more than doubled the abundance of atmospheric CH₄ over the industrial era (e.g. Etheridge et al., 1998). The CH₄ emission histories have been estimated based on human activity statistics combined with emission factors (Stern and Kaufmann, 1996; van Aardenne et al., 2001). Such historical emission inventories have been examined by atmospheric chemistry transport modelling (Houweling et al., 2000; Monteil et al., 2011; Ghosh et al., 2015), in combination with the records of atmospheric CH₄ mole fraction reconstructed from polar ice cores (Blunier et al., 1993; Nakazawa et al., 1993; Etheridge et al., 1998; MacFarling Meure et al., 2006; Sapart et al., 2012) and air extracted from porous snow layers at the top of ice sheets (firn; Francey et al., 1999; Buizert et al., 2012; Sapart et al., 2013).

A large fraction of natural and anthropogenic CH₄ sources resides in the Northern Hemisphere, and thus, the atmospheric CH₄ trend of the Northern Hemisphere can provide important information on the evolution of anthropogenic CH₄ emissions and the variations in natural CH₄ emission in response to climatic variability. The interhemispheric gradient of CH₄ mole fraction is also key to the allocation of CH₄ emissions between both hemispheres (e.g. Dlugokencky et al., 2003; Ghosh et al., 2015; Chandra et al., 2021).

Figure 1 Atmospheric CH₄ mole fraction data covering the last 300 years. Symbols in brown are from the Southern Hemisphere; crosses and open circles are from the Law Dome ice core and firn, respectively (Etheridge et al., 1998; MacFarling Meure et al., 2006; Rubino et al., 2019). Open squares are from Kennaook / Cape Grim (CGO; MacFarling Meure et al., 2006), and dots are from the South Pole (SPO; ftp://aftp.cmdl.noaa.gov/, last access: 10 December 2021). Other coloured data (except for yellow) are from the Northern Hemisphere. Closed green squares, triangles, and diamonds are from the EUROCORE (Blunier et al., 1993), EUROCORE and GISP2 (Etheridge et al., 1998), and Site J (Nakazawa et al., 1993) ice cores. Red dots are the NEEM-S1 ice core (Rhodes et al., 2013), blue dots are from Alert, Canada (ALT), and Barrow, Alaska (BRW; ftp://aftp.cmdl.noaa.gov/), and light blue and yellow lines are the atmospheric scenarios prepared for the NEEM firn modelling (Buizert et al., 2012) and the CMIP6 experiments (Meinshausen et al., 2017), respectively. For the latter, Arctic and Antarctic scenarios with higher and lower mole fractions are shown, respectively.

Systematic measurements of atmospheric CH₄ mole fraction began in the 1980s. In Fig. 1, CH₄ measurements at two Arctic sites, i.e. Barrow, Alaska (BRW), and Alert, Canada (ALT), whose records start in 1983 and 1985, respectively, are shown (blue dots; 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 CH₄ data only from around 1980, which means that some reconstruction methodology is required to infer atmospheric CH₄ mole fraction variations before that time. For this purpose, air extracted from ice cores and firn layers has been measured. Figure 1 also presents CH₄ mole fractions analysed in ice cores from Greenland (closed green 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 agreement with each other until ∼ 1900, after which the number of data and the consistency among the records are poor. Continuous measurements from North Greenland Eemian Ice Drilling (NEEM)-S1 ice core (Rhodes et al., 2013) presented CH₄ mole fractions before ∼ 1960 (red in Fig. 1), but their data are notably higher than other ice core data after ∼ 1850 (Blunier et al., 1993; Etheridge et al., 1998; Nakazawa et al., 1993). Therefore, the inconsistency among the different data sets indicates that the period between ∼ 1900 and ∼ 1980 is not covered reliably either by direct observations or ice core reconstructions.

In comparison to these data, higher-resolution CH₄ data are available from Antarctica (brown in Fig. 1). The comprehensive Law Dome ice core and firn data set (Etheridge et al., 1998; MacFarling Meure et al., 2006; Rubino et al., 2019) almost continuously covers the last 300 years and is well connected to the direct measurements at the South Pole (SPO; data provided by NOAA/ESRL/GML) and Kennaook / Cape Grim (CGO; Etheridge et al., 1998; MacFarling Meure et al., 2006). Such consistency and continuity among the data sets suggest that the Antarctic CH₄ data can serve as a good reference to represent the global atmospheric CH₄ trend over the past few centuries.

Despite the limited observational data as described above, two synthetic data sets for the Arctic historical CH₄ mole fractions are currently available (see Sect. 3.2 for details). Buizert et al. (2012) prepared Arctic historical trends of mole fractions of atmospheric trace gases, including CH₄ (light blue in Fig. 1), which in turn was used to constrain the gas diffusivity profile in firn at the NEEM site. Considered as the most likely atmospheric CH₄ trend for the northern high latitudes, this scenario was treated as a known history, by which the diffusivity profiles in firn were tuned in firn air transport models (Witrant et al., 2012; Trudinger et al., 2013). The other is the composite data set prepared for use in the CMIP6 experiments (yellow in Fig. 1; Meinshausen et al., 2017). It is seen that their scenario for the northernmost latitude (higher mole fraction; yellow line in Fig. 1) follows the NEEM-S1 ice core data set (red) and is inconsistent with the Buizert et al. (2012) scenario (light blue).

For constraining the global CH₄ budget, isotope ratios of CH₄ could provide useful information for relative contributions of different types of CH₄ sources. Sapart et al. (2013) examined the reconstruction of stable carbon isotope ratio (δ¹³C) of atmospheric CH₄ using firn air measurements from both the Northern and Southern hemispheres. They concluded that, with the available firn measurements and understanding of firn air transport, it is difficult to consistently reconstruct the past trend of δ¹³C of CH₄ because of multiple reasons, including uncertainty in the atmospheric CH₄ mole fraction scenario. Among many important and uncertain factors, the accurate reconstruction of the atmospheric CH₄ mole fraction is particularly important because the trend in the mole fraction can lead to a significant signal in the modelled δ¹³C profile in firn due to the difference in the molecular diffusion coefficient, even in the absence of a temporal trend in atmospheric δ¹³C. Reducing the uncertainty of the historical trend of the CH₄ mole fraction is therefore also important for utilising isotope data for a better understanding of the historical changes of different CH₄ emission categories.

In this study, we present a set of mole fractions of CH₄ and other trace gases in firn air samples collected at the NGRIP (North Greenland Ice Core Project) site. Using the available atmospheric scenarios, we simulate the depth profiles of trace gases in the NGRIP firn with our firn air transport model and those in the NEEM firn reported previously (Buizert et al., 2012). It will be shown in series of firn modelling in this study that, compared to other trace gases, CH₄ is uniquely underconstrained in its Arctic history. We examine a variety of modelling cases for different diffusivity profiles and reconstruct the Arctic atmospheric CH₄ over the late 20th century using an iterative dating approach (Trudinger et al., 2002). The reconstructed CH₄ trends from both NGRIP and NEEM firn air samples are evaluated by comparison to the atmospheric CH₄ scenarios. The uncertainty of the Arctic atmospheric CH₄ history for use in firn air modelling is discussed.

2 Experimental method

Firn air was sampled at the Greenland site NGRIP (75.10^∘ N, 42.32^∘ W; 2959 m a.m.s.l. – above mean sea level) in May–June 2001. Mean accumulation, surface density, temperature, and pressure are 179 kg m⁻² yr⁻¹, 300 kg m⁻³, 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 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. (2021a), only brief descriptions of the relevant data presented in this study are given here. CH₄ mole fractions of the firn air samples were measured using a gas chromatograph (Agilent 6890, Agilent Technologies, Inc.) equipped with a gas chromatograph flame ionisation detector (GC-FID) at Tohoku University (TU), with a reproducibility of 2 ppb (parts per billion; Umezawa et al., 2014). The CH₄ mole fractions were determined against our working standard gases that were calibrated on the TU1987 CH₄ scale (Aoki et al., 1992; Umezawa et al., 2014; Fujita et al., 2018). The difference between the TU1987 CH₄ scale and the WMO (World Meteorological Organization) CH₄ mole fraction scale (on which the NEEM CH₄ data were measured) is estimated to be ∼ 0.5 ppb at the current atmospheric CH₄ levels (Fujita et al., 2018). Oyabu et al. (2020) reported that ice core data analysed 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 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 CO₂ and SF₆ mole fractions, respectively, by using a nondispersive infrared gas analyser (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 (parts per million) for CO₂ and 0.09 ppt (parts per trillion) for SF₆, and mole fractions of both gases are reported on the TU2010 CO₂ and TU2002 SF₆ scales, respectively (Sugawara et al., 2018). Although not presented in this study, the firn air samples were also previously analysed for nitrous oxide (N₂O) and its isotope ratios (Ishijima et al., 2007).

As part of this study, the NGRIP firn air samples were newly analysed for selected halocarbons (CFC-11, CFC-12, CFC-113, and CH₃CCl₃) 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 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 Carboxen 1000/Tenax TA at −100 ^∘C. The trap was then heated to 180 ^∘C to inject the trapped gases onto a PoraBOND Q separation column for subsequent analysis on MS. Mole fractions of individual halocarbons are determined against a working standard gas (compressed dry air) that was calibrated against synthetic standards on the NIES-08 (National Institute for Environmental Studies) scales.

We also use a suite of trace gas measurement data from the NEEM site (77.45^∘ N, 51.06^∘ W; 2450 m a.m.s.l.). The firn air samples were collected in July 2008. The details of firn air sampling and gas measurements have been described by Buizert et al. (2012). The depth profile data of all the above trace gases (CH₄, CO₂, SF₆, CFC-11, CFC-12, CFC-113, and CH₃CCl₃) from the NEEM firn are used in this study together with HFC-134a and ¹⁴CO₂ data that are available from NEEM but not for NGRIP.

3 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 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 the following three zones: a convective zone (CZ), a diffusive zone (DZ), and a lock-in zone (LIZ; 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 is estimated to be below 2 m at NGRIP (Kawamura et al., 2006) and 4.5 m at NEEM (Buizert et al., 2012). In DZ, which is sufficiently isolated from the surface turbulence, the movement of air is governed by molecular diffusion. Gravitational enrichment according to the barometric equation (i.e. linear increases of δ¹⁵N of N₂ and δ¹⁸O of O₂) occurs with depth and stops at the top of LIZ (Sowers et al., 1992; Schwander et al., 1993; Kawamura et al., 2006). The top of LIZ (lock-in depth) is at a depth of 63 m, which is coincidently the same at NGRIP (Kawamura et al., 2006) and NEEM (Buizert et al., 2012). 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 was successfully made at 77.71 m at NGRIP (Kawamura et al., 2021a) and 77.75 m at NEEM (Buizert et al., 2012), and total pore closure is considered in the deeper layers in our modelling.

3.1 Modelling firn air transport

We use a one-dimensional diffusion model that has been used for the reconstruction of isotope ratios of CO₂ and N₂O (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 (Schwander et al., 1993) and a bubble-trapping process (Rommelaere et al., 1997). Air movement in the firn is driven by molecular diffusion and a gravitational effect. In other words, a trace gas flux (F) in firn is expressed by the following:

$$\begin{array}{}\text{(1)}& F=-D\left\{s{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{c}{s}}\right)-{\displaystyle \frac{mgc}{RT}}\right\},\end{array}$$

where D is the effective diffusivity of a trace gas molecule (m² s⁻¹), the variables s, c, and T are open porosity (unitless), trace gas molar concentration (mol m⁻³), and firn temperature (K), respectively, and the constants m, g, and R are the mass number of the trace gas (kg mol⁻¹), the acceleration of gravity (m s⁻²), and the gas constant (J mol⁻¹ K⁻¹), 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 the following:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial c}{\partial t}}+{\displaystyle \frac{\partial \left(vc\right)}{\partial z}}+{\displaystyle \frac{\partial F}{\partial z}}+rc=\mathrm{0},\end{array}$$

for the open pore space and, in the following, by:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial c_{\mathrm{b}}}{\partial t}}+{\displaystyle \frac{\partial \left(v_{\mathrm{f}}{c}_{\mathrm{b}}\right)}{\partial z}}-rc=\mathrm{0},\end{array}$$

for bubbles. Here c and c_b are trace gas molar concentrations in the open pore space and bubbles, respectively. The vertical speed of air (m yr⁻¹) in the open pore space v is distinguished from that of firn itself v_f. The vertical speed of firn v_f(z) is simply given by dividing the accumulation rate by the firn density under the assumption of the steady-state densification of firn. At the transition zone where the open pore air is gradually trapped into bubbles, mass conservation is given by using a bubble-trapping rate r (s⁻¹), 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 s_c in the NGRIP firn was calculated by the empirical equation, given by Schawander (1989), as follows:

$$\begin{array}{}\text{(4)}& s_{\mathrm{c}}=s_{\mathrm{t}}\exp \left\{\mathrm{75}\left({\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{close}}}}-\mathrm{1}\right)\right\},\end{array}$$

where s_t is the total porosity, and ρ and ρ_close are the density of firn and that at the close-off depth (kg m⁻³), respectively. For the NEEM firn, the closed porosity s_c was calculated as per Buizert et al. (2012), who followed the parameterisation by Goujon et al. (2003), as follows:

$$\begin{array}{}\text{(5)}& s_{\mathrm{c}}=\mathrm{0.37}{s}_{\mathrm{t}}{\left({\displaystyle \frac{s_{\mathrm{t}}}{s_{\mathrm{close}}}}\right)}^{-\mathrm{7.6}},\end{array}$$

where s_close is the mean close-off porosity (9.708 × 10⁻²), according to Buizert et al. (2012).

3.2 Atmospheric scenarios

To simulate depth profiles of trace gases in firn, atmospheric histories of the target gases are required. In this study, we used atmospheric histories prepared by the NEEM firn air modelling (Buizert et al., 2012) and by the CMIP6 experiments (Meinshausen et al., 2017) for all the trace gases presented in this study for the NGRIP and NEEM firn (CH₄, CO₂, SF₆, CFC-11, CFC-12, CFC-113, HFC-134a, CH₃CCl₃, and ¹⁴CO₂). Since Meinshausen et al. (2017) provide latitudinally gridded data sets, their historical data for the northernmost latitude (82.5^∘ N) are used. Note that the ¹⁴CO₂ history for CMIP6 is available from another study (Graven et al., 2017), which is also used in this study. The ¹⁴CO₂ data by Graven et al. (2017) are available in Δ¹⁴CO₂ for three zonal bands (Northern Hemisphere, the tropics, and Southern Hemisphere) and were converted to ¹⁴CO₂ mole fraction, as in Buizert et al. (2012). These atmospheric scenarios (hereafter referred to as the BZ and CMIP6 scenarios) are compared in Fig. 2. The BZ (light blue) and CMIP6 (blue) Arctic scenarios show compatible historical trends in many trace gases. There are, however, slight differences in some trace gases between the two scenarios, e.g. SF₆ and CFC-12, but we later show that these differences do not cause significant biases in reproducing their depth profiles at the NGRIP and NEEM firn sites.

Figure 2 Atmospheric scenarios of various trace gases used in this study. The Arctic histories prepared for the NEEM firn air modelling (Buizert et al., 2012) are in light blue. The Arctic and Antarctic histories for CMIP6 (Meinshausen et al., 2017) are in blue and red, respectively. For ¹⁴CO₂, the CMIP6 historical data are from Graven et al. (2017). The interpolar differences (IPDs) calculated from the Arctic histories with respect to the Antarctic CMIP6 histories are also shown (right axes).

In contrast, the Arctic CH₄ histories by the two studies differ considerably with maximum difference of ∼ 85 ppb around 1910 (Fig. 2b; see also Fig. 1). The two CH₄ histories show similar trends after 1960, but before this, the data sets diverge further into the past. This disagreement is also clear in the interpolar difference (IPD), which is calculated relative to the CMIP6 histories for Antarctica (right axes). Note that the CMIP6 scenario for the Antarctic latitude was constructed based on the Law Dome ice core data (see the agreement in Fig. 1). While the BZ scenario shows a gradual increase in IPD, the CMIP6 scenario indicates almost constant values of IPD at ∼ 130 ppb over the 20th century. This difference stems from the different methodologies that were used to produce the respective scenarios. As described by Buizert et al. (2012), the BZ CH₄ scenario was constructed by adding the presumed IPD to the Antarctic history (the Law Dome data), where the IPD was assumed to be proportionally correlated with the growth rate of CH₄. This seems to be a 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; Ghosh et al., 2015; Chandra et al., 2021). Meinshausen et al. (2017) compiled historical measurement records from the worldwide networks and Antarctic/Greenland ice core and firn samples and constructed latitudinally gridded data sets of various greenhouse gases. For the historical trend of CH₄, they relied on the data set from the NEEM-S1 ice core by Rhodes et al. (2013) to produce the atmospheric histories for the Northern Hemisphere. They used the 5-yearly averaged values with outliers removed to represent lower bounds of the raw data points, as shown in Fig. 1. Note that the Rhodes et al. (2013) data set was not available when the BZ scenario was constructed. Thus, CH₄ is unique because currently it has two diverging synthetic Arctic histories that are only loosely constrained by observational data.

This study assumes that the atmospheric scenarios for trace gases other than CH₄ are known with sufficient accuracy. The scenarios of the individual trace gases have inherent uncertainties, but the comparisons of the two available scenarios (BZ and CMIP6) indicate that the data sources for other gases do not show inconsistent variations, as seen in CH₄. It should be, however, noted that, except for CO₂, many trace gases lack observational data for the early 20th century. Thus, both scenarios are, to a large extent, based on the same data sources.

3.3 Effective diffusivity

We follow the previous firn air studies in which the effective diffusivity in firn is optimised with an iterative method so as to minimise the difference between the simulated and observed depth profiles of CO₂ (Sugawara et al., 2003; Ishijima et al., 2007). In these previous studies, an initial guess of the depth profile of effective diffusivity for CO₂, D_init(z), was calculated as follows:

$$\begin{array}{}\text{(6)}& D_{\mathrm{init}}\left(z\right)=D_{\mathrm{0}}{\left({\displaystyle \frac{T}{\mathrm{253}}}\right)}^{\mathrm{1.85}}\left({\displaystyle \frac{\mathrm{1013}}{p}}\right)\left\{\mathrm{1.7}s\left(z\right)-\mathrm{0.2}\right\},\end{array}$$

where s(z) and D₀ represent the open porosity at a depth z and the diffusion coefficient of CO₂ at 253 K and 1013 hPa, respectively. D₀ was set to 1.247 × 10⁻⁵ (m² s⁻¹) according to Trudinger et al. (1997). p is the mean atmospheric pressure (Pa). The bulk density (kg m⁻¹) was determined by measuring the dimension and weight of cylindrically cut firn core samples (Kawamura et al., 2006). The effective diffusivity of CO₂, 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 due to the scaling factors. In this study, the effective diffusivity profile prepared for the NGRIP firn by Ishijima et al. (2007) is referred to as the initial diffusivity, and it was modified to improve the reproducibility of our newly measured trace gas profiles. For simulating trace gas profiles for the NEEM firn, we began with the effective diffusivity profiles available from Buizert et al. (2012). Those effective diffusivity profiles, which were originally optimised for individual firn air transport models that participated in that study, were modified and used for simulating the various trace gas profiles reported for the NEEM firn. The various diffusivity profiles were constructed by modifying the original profiles at a certain range of depths in a stepwise empirical manner; the depth range of the diffusivity profile key to improve reproducibility of trace gas depth profiles was first diagnosed and then the diffusivity was perturbed up and down in the depth range to the degree in which the corresponding simulated profiles do not deviate substantially. Although 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.

Figure 3 The 100 sets of effective diffusivity profiles of CO₂ in the NGRIP (a, b, where panel a is on a linear scale and panel b is on a log scale) and NEEM firn (c, d). The initial diffusivity profile (Ishijima et al., 2007) is shown by a black dotted line (NGRIP only) and modified diffusivity profiles are in colours. The diffusivity profiles whose corresponding mole fraction profiles have root mean square deviation (RMSD) values of < 1.0 are coloured red and the others light blue.

We eventually prepared 100 different sets of diffusivity profiles so as to cover a considerable range of diffusivity. Each set of the diffusivity profiles was evaluated based on the root mean square deviation (RMSD) between the model and data, according to Buizert et al. (2012). All sets of effective diffusivity profiles for the NGRIP and NEEM firn sites are shown in Fig. 3 (top and bottom panels, respectively). The different colours of the diffusivity profiles will be explained later. Those diffusivity profiles were evaluated against the observed trace gas profiles, which were regarded as constraints. Note that CH₄ was not used in the evaluation, and the atmospheric scenarios of other trace gases are assumed to be known with sufficient accuracy to infer a range of acceptable diffusivity profiles that reproduce the depth profiles of the firn air composition.

3.4 Performance of the firn air transport model

To validate our firn air transport model, we began by simulating 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 the BZ scenarios as per their model intercomparison. The simulated depth profiles of the nine trace gases are presented in Fig. 4 and compared with those by other models presented in Buizert et al. (2012). The results confirm that the performance of our model is comparable to those by other groups. As a measure of the model performance, Buizert et al. (2012) compared RMSD, 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 for the NEEM EU borehole. This RMSD value was achieved with an effective diffusivity profile that was prepared by modifying the profile originally optimised for the CIC (Centre for Ice and Climate) model at a certain range of depths. Note that the RMSD value here was calculated to include CH₄, as per Buizert et al. (2012), but, as described in Sect. 3.5, CH₄ is excluded in calculation of RMSD in the following sections.

Figure 4 Modelled depth profiles of various compounds in the NEEM firn. Closed circles in black are the observed data with uncertainty estimated by Buizert et al. (2012). Our model results are shown with black solid lines, and the results from other models are in colours. Different models are labelled by institution, according to Buizert et al. (2012), as follows: CIC (Centre for Ice and Climate), CSIRO (Commonwealth Scientific and Industrial Research Organisation), INSTAAR (Institute of Arctic and Alpine Research), LGGE (Laboratoire de Glaciologie et Géophysique de l'Environnement), OSU (Oregon State University), and SIO (Scripps Institution of Oceanography).

3.5 Reconstruction of CH₄ history

Our modelling procedure for reconstructing the atmospheric histories of CH₄ was as follows:

1. To represent atmospheric trends of different trace gases in the Arctic region, we began by employing the two sets of atmospheric scenarios (Sect. 3.2). The firn transport model calculates depth profiles of the various trace gases at the NGRIP and NEEM firn sites using the large set of modified (N=100) effective diffusivities described in Sect. 3.3. The simulation case with each diffusivity profile was evaluated based on RMSD. As we aim to estimate the historical atmospheric trend of CH₄, the RMSD-based evaluations were made using all of the available trace gas data, excluding CH₄. In the RMSD calculation, we used the measurement uncertainties of 0.2 ppm for CO₂, 0.2 ppt for SF₆, 1.1 ppt for CFC-11, 3.3 ppt for CFC-12, 0.6 ppt for CFC-113, and 3.2 ppt for CH₃CCl₃ for the NGRIP firn. For the NEEM firn, we employed the uncertainties provided by Buizert et al. (2012). It should be noted that the present study does not follow the uncertainty estimation as done by Buizert et al. (2012). They indicated that uncertainties in the atmospheric scenarios and measurement uncertainties are the two largest contributors to the total uncertainties for individual data points for the NEEM firn. We consider that the uncertainties in the atmospheric scenarios are appreciably examined through comparisons of series of simulations using the two independent scenarios.

2. We ran the model with 100 different sets of diffusivity profiles to calculate the depth profile of CH₄. Note that, based on the earlier step, we know the diffusivity profiles that generate reasonable firn air profiles for the trace gases other than CH₄. Every diffusivity profile was used in combination with the firn air CH₄ data for reconstructing an atmospheric CH₄ history. We employed an iterative dating approach (Trudinger et al., 2002), where the initial atmospheric scenario (the BZ scenario) was modified to improve model reproducibility of the CH₄ depth profile (see below). The corrected atmospheric CH₄ scenarios were then compared to the original scenarios (BZ and CMIP6) for further discussion.

The iterative dating for CH₄ was performed as follows:

• i.

  The depth profile of CH₄ was calculated with the initial atmospheric CH₄ scenario.

• ii.

  The modelled CH₄ mole fraction, calculated in step (i), was compared to the input atmospheric CH₄ scenario, and an effective age at each sampling depth was determined as the time when the modelled CH₄ agreed with a value in the atmospheric CH₄ scenario. It is noted that the smoothing spline curve applied to the BZ CH₄ scenario was used for the calculation of the effective age, as the input scenario with seasonal variation (Fig. 2) would not allow the effective age to be uniquely determined.

• iii.

  A new atmospheric CH₄ scenario was constructed by assigning the observed CH₄ mole fraction, at each sample depth, to the effective age determined in step (ii). The observed CH₄ versus the effective age data set was interpolated by a smoothing spline function, and it is considered as a revised atmospheric CH₄ scenario.

• iv.

  The depth profile of CH₄ was again calculated with the revised atmospheric CH₄ scenario constructed in step (iii).

• v.

  The above steps (ii–iv) 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.

4 Result

4.1 Initial simulations: NGRIP

Figure 5 presents the simulation results with the initial diffusivity in comparison to the observed profiles for the six trace gases, excluding CH₄ (CO₂, SF₆, CFC-11, CFC-12, CFC-113, and CH₃CCl₃) for the NGRIP firn. It is again noted that the initial diffusivity profile was tuned only for the depth profile of CO₂ (Ishijima et al., 2007). As seen in this figure, measured profiles of these trace gases (except CH₃CCl₃) show gradual decreases with depth in the DZ and sharp decreases in the LIZ. In contrast, CH₃CCl₃ increases with depth in the DZ and sharply decreases in the LIZ. The difference in the depth profile pattern among species is due to their different historical atmospheric trends. It is known that, since the mid-20th century, the atmospheric mole fractions of the five trace gases (CO₂, SF₆, CFC-11, CFC-12, and CFC-113) have increased either monotonically or shown peak/slowed increase in the early 1990s (Sturrock et al., 2002; Martinerie et al., 2009). In contrast, CH₃CCl₃ increased until the early 1990s and has rapidly decreased since then (Sturrock et al., 2002; Rigby et al., 2017), which is also observed in Fig. 2. Our simulation reproduces the observed depth profiles of these six trace gases in the NGRIP firn fairly accurately by using the BZ scenarios.

Figure 5 Depth profiles of CO₂, SF₆, CFC-11, CFC-12, CFC-113, and CH₃CCl₃ in the NGRIP firn. The measurement and model results are shown by red circles and black solid lines, respectively (left axis). The measurement uncertainties are shown as vertical error bars, though in many cases they are smaller than the circle sizes. Black solid lines show the modelled profiles with the initial diffusivity and the BZ scenarios. The black dashed lines indicate the profiles calculated with the atmospheric scenarios for Antarctica (red lines in Fig. 2). The model–data differences are also shown (right axis). The vertical solid line in each panel indicates the upper depth of LIZ.

It is interesting to note that the depth profiles of CH₃CCl₃ at NGRIP and NEEM are remarkably different. Whereas the NEEM data show a relatively sharp CH₃CCl₃ peak in the LIZ (∼ 65 m; Fig. 4), 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 signal of the maximum atmospheric CH₃CCl₃ 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, 7 years later, at the NEEM site, such a signal was found deeper in the LIZ, where the age of air changes rapidly with depth in both deeper and shallower sides. We emphasise that, despite the differences in the depth profiles of CH₃CCl₃ at the two sites, our simulations reproduce the profiles measured at both sites well when using the same atmospheric CH₃CCl₃ scenario.

4.2 Sensitivity to diffusivity profile: NGRIP

Figure 5 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 modelling, e.g. a depth profile of diffusivity. To examine the impact of diffusivity modification on the simulated depth profiles and their agreements with the data, we examined the 100 sets of modified diffusivity profiles (Fig. 3). It was found that, to reduce the model–data difference in the LIZ (Fig. 5), the diffusivity needs to be increased in the shallower layers compared with 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 (Fig. 6). In this figure, the modelling results with the modified diffusivity profiles are shown in colours 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. The RMSD values are as little as 0.51 for a particular case. In Fig. 6 and associated figures, the model results with RMSD of < 1.0 are coloured red and other cases are coloured light blue.

Figure 6 Depth profiles of the six trace gases below 50 m depth in the NGRIP firn. Black circles indicate the measurements, and solid lines in colours are the model results with different diffusivity profiles and the BZ scenarios (left axis). Black dotted lines are modelled results with the initial diffusivity. Also shown are the model–data differences (right axis). See the text for the differences in the line colours.

4.3 Sensitivity to atmospheric scenario: NGRIP

In Fig. 7, modelled depth profiles with different sets of the atmospheric scenarios (BZ and CMIP6) are compared for trace gases other than CH₄. For simplicity, only the results with RMSD of < 1.0 are presented. This figure shows that the differences in the atmospheric histories (Fig. 2) produces relatively small differences in the depth profiles in the firn. There are small offsets due to the differences of the histories in some gases; the difference in the SF₆ history before 1980 (< 0.2 ppt) corresponds to the small (< 0.1 ppt) offsets below 65 m, significant differences in the histories of CFC-11 (< 10 ppt) and CFC-12 (< 20 ppt) for 1960–1990 resulted in the overall offsets (roughly < 5 and < 10 ppt, respectively) below 50 m, and the difference in the CH₃CCl₃ history (< 10 ppt) in the early 1990s produced the offsets (∼ 3 ppt) above 66 m. The smaller offsets in the calculated depth profiles than in the input atmospheric scenarios are due to the smoothing effect of diffusion in the firn layers. We calculated the difference between the modelled profiles with the two scenarios for individual diffusivity cases and found that those differences in the LIZ are within the measurement precisions for SF₆, CFC-113, and CH₃CCl₃. They are a bit larger for CO₂, CFC-11, and CFC-12, which are up to 5 times the respective measurement precisions.

Figure 7 Same as Fig. 6 but for a comparison of the modelled profiles in the NGRIP firn with different atmospheric scenarios of BZ and CMIP6, which are shown in light blue and green, respectively. Only the simulation results with RMSD values of < 1.0 are shown.

4.4 Simulations of CH₄: NGRIP

Modelled depth profiles of CH₄ are shown in Fig. 8. These calculations were made with the 100 diffusivity profiles and the atmospheric CH₄ scenarios of BZ (left) and CMIP6 (right). It is interesting to note that the characteristics of the model–data difference for CH₄ are different from those for the other six trace gases (Fig. 6). For the other trace gases, the initial simulation showed increased model–data differences around 70 m, and they were reduced with some modified diffusivity profiles. For CH₄, the model run with the initial diffusivity profile and the BZ scenario reproduces the observed CH₄ profile quite well down to ∼ 73 m but significantly overestimates by > 20 ppb for the lowest three depths (the black dotted line in the left panel). Using the modified diffusivity profiles that allowed better agreements for the other six trace gases (red lines), we find larger model–data CH₄ differences than in the initial simulation. These features are also seen for the simulation results with the CMIP6 scenario (right panel), but the overestimate in the LIZ is more pronounced because the CH₄ mole fractions for the early- to mid-20th century are higher in the CMIP6 than in the BZ scenario. In comparison to other gases (Sect. 4.3), the largest impact of the scenario difference occurs in CH₄. In the LIZ, the difference between depth profiles with the two scenarios reaches to 5–10 times the measurement precision, namely > 30 ppb in some diffusivity cases.

Figure 8 Same as Fig. 6 but for CH₄ in the NGRIP firn. Simulation results with the atmospheric scenarios of BZ and CMIP6 are shown in panels (a) and (b), respectively. The black dashed lines indicate the profiles calculated with the atmospheric scenarios for Antarctica (red line in Fig. 2).

4.5 Sensitivity to IPD: NGRIP

Reflecting geographical source/sink allocation and atmospheric lifetime, different trace gases show a different magnitude of IPD, as shown in Fig. 2. It would be possible that, for some gases, the northern to southern difference in mole fraction has little influence on depth profile in firn due to the relatively small interhemispheric gradient. To examine magnitude of the sensitivity of the depth profiles to the IPD for individual trace gases, we calculated the depth profiles that would be expected if the Antarctic atmospheric scenarios (red lines in Fig. 2) are given to force the firn model for the NGRIP firn. The results are shown by dashed lines in Fig. 5. This sensitivity experiment shows that IPD causes significant biases larger than the measurement precisions of the respective trace gases. We calculated the difference between the simulations for the BZ scenario (solid line) and the Antarctic scenario (dashed line) and found that sensitivities to the IPD for these six trace gases are no more than 20 times the respective measurement uncertainties. Such relative sensitivities of these suite of gases to the IPDs are much smaller than that of CH₄, which reaches 40 times the measurement uncertainty. As seen in Fig. 8, the calculated CH₄ profile with the Antarctic scenario for the NGRIP firn is aligned ∼ 120 ppb below the original simulation, showing the pronounced impact on CH₄.

4.6 Sensitivity to atmospheric scenario: NEEM

In addition to the above simulations for the NGRIP firn, we examined the impact of the difference of the atmospheric scenarios on the depth profiles for the NEEM firn. In Fig. 9, comparisons between simulations with the different scenarios for the nine trace gases, including CH₄, are presented. We found some common characteristics at the NGRIP and NEEM firn sites. While the difference between the simulations with the two scenarios are relatively small for most trace gases, the large difference between the CH₄ scenarios (i.e. higher CH₄ mole fraction in the CMIP6 scenario; see Fig. 2) results in an overestimation of CH₄ mole fraction in the LIZ (> 63 m) in the modelled profiles with the CMIP6 scenario. In the LIZ, the difference in the depth profiles with the two scenarios exceeds 30 ppb, which is comparable to the magnitude observed in the NGRIP firn (Fig. 8).

Figure 9 Same as Fig. 7 but for the NEEM firn. Only the simulation results with RMSD values of < 1.0 are shown.

4.7 Reconstruction by iterative dating

As shown in the earlier sections, CH₄ shows uniquely high uncertainty in the firn modelling compared to other trace gases. First, simulated depth profiles of CH₄ do not show improved reproducibility with effective diffusivity profiles that satisfy depth profiles of other gases (Sect. 4.2 and 4.4). Second, while other gases show limited sensitivities to the different atmospheric scenarios, inconsistency between the currently available two scenarios have particularly large impact on the simulated depth profile of CH₄ (Sect. 4.3 and 4.4). Third, sensitivity to IPD in the atmospheric scenarios is most pronounced for CH₄ compared to other trace gases (Sect. 4.5). These characteristics are common for both NGRIP and NEEM sites (Sect. 4.6).

Taking all these unique characteristics of CH₄ into account, we explore the possibility of reconstructing the Arctic CH₄ mole fractions by the iterative dating method. This approach has a fixed diffusivity profile (assumed to be correct) and aims to find an acceptable atmospheric history to reproduce the firn air depth profile. It is considered that the modified diffusivity profiles with the low RMSD values are adequately evaluated by the trace gases, except CH₄, and that the model–data mismatch in the CH₄ modelling is therefore attributable to uncertainty in the atmospheric CH₄ scenario. The historical atmospheric CH₄ variations obtained by the iterative dating method are presented in Fig. 10 for the 100 modelling cases of the modified diffusivity. The different simulation cases are coloured in the same manner as in the earlier figures according to the RMSD (Figs. 3, 6, and 8). Note that the reconstruction cases coloured in light blue are considered to be less likely due to poorer reproduction of the depth profiles (Figs. 6 and 8). The NGRIP reconstruction results (red; Fig. 10a) are in good agreement with the BZ scenario after around 1980. For the earlier period, however, the upper bounds of the reconstructions (line connecting far-left red circles) match with the BZ scenario, and the overall range of acceptable histories is below the BZ scenario (circles and shades in red). On the other hand, the NEEM reconstruction results show that the range of acceptable histories are distributed closely around the original BZ scenario or below it after around 1960, whereas the lower bounds of the reconstruction (line connecting far-right red circles) align with the BZ scenario before that (Fig. 10b).

Figure 10 Results of reconstructions by the iterative dating method for the model cases with the 100 different diffusivity profiles (circles coloured in the same manner as in Figs. 3, 6, and 8) from the (a) NGRIP and (b) NEEM firn air and (c) the ranges of reconstructions from both firn sites. The red shades indicate the range of CH₄ mole fraction trends reconstructed using the diffusivity profiles that better reproduce the trace gases, except for CH₄. Black and grey solid lines are the atmospheric CH₄ scenarios of BZ and CMIP6, respectively, and the black dotted line is the smoothed Antarctic CH₄ mole fraction (spline curve) from the Antarctic data shown in Fig. 1. The thick horizontal bars in red and blue in panel (c) correspond to red shades in panels (a) and (b). The light blue thin bars in panel (c) indicate the range of the reconstructions for cases in which ¹⁴CO₂ data are excluded for evaluations of the diffusivity profiles. The inset in panel (c) shows the magnified view of the NGRIP upper bounds (red) and the NEEM lower bounds (blue) in comparison to the two scenarios.

It should be stressed that the reconstructions for the time period before 1980 rely on CH₄ data from the lowest four or five depths of the NGRIP and NEEM firn air. The differences between the initial and corrected atmospheric CH₄ scenario from the three deepest data for the NGRIP firn are up to ∼ 100 ppb. As seen in Fig. 10, such a reduction in the Arctic atmospheric CH₄ scenario over the period would result in alignment with the atmospheric CH₄ trend in Antarctica inferred from the Law Dome and other data sets (black dotted line). It is again noted that all the reconstruction cases coloured in red used diffusivity profiles that yield relatively good model reproducibility with RMSD values of < 1.0. It is therefore seen that iterative-dating-based reconstruction from the NGRIP firn data suggests lower CH₄ mole fraction than the BZ or CMIP6 scenarios from the 1950s to 1970s in any case, albeit with large uncertainty. It is, however, noted that, as shown in Fig. 10c, this result is not fully compatible with that from the NEEM firn. In particular, the discrepancies between the reconstructions from both firn data diverge with time before 1970. When compared to these reconstructions, the BZ scenario (black line) tracks the overlapping ranges of both reconstructions, while the CMIP6 scenario passes above them before around 1955.

5 Discussion and conclusion

5.1 CH₄ history

In Sect. 4, we have shown that, with the two atmospheric scenarios (BZ and CMIP6), depth profiles of CO₂, SF₆, CFC-11, CFC-12, CFC-113, and CH₃CCl₃ in the NGRIP and NEEM firn are reproduced with sufficient accuracy by using a range of modified diffusivity profiles (Figs. 6, 7, and 9). This suggests that the atmospheric scenarios of these trace gases used in the modelling are consistent with the depth profiles observed at both firn sites in Greenland (NEEM and NGRIP). In contrast, the observed CH₄ profile in the NGRIP firn was not accurately reconciled using the two atmospheric scenarios and the diffusivity profiles that allow adequate reproducibility for the trace gases (except CH₄). This suggests either that the Arctic atmospheric scenarios of CH₄ are uncertain or that the diffusivity profile of the NGRIP firn is underconstrained. We explored the correction of the atmospheric scenario of CH₄ by using the iterative dating approach. This method improves agreement to the observed CH₄ depth profile, with an implicit assumption that the diffusivity profile in each case is correct. Although uncertainty due to the underconstrained diffusivity profile in the LIZ is large, this attempt for the NGRIP firn suggested that the CH₄ mole fractions over the period 1950–1980 could be ∼ 100 ppb lower than the original BZ scenario (Fig. 10). In contrast, the iterative dating reconstruction for the NEEM firn agrees with the BZ scenario. Although the spread of the reconstructions is large, particularly for the NGRIP firn, it was found that the BZ scenario passes within the ranges of the reconstructions from both firn data but that the CMIP6 scenario is notably higher than the reconstructions for the period before 1960.

Whereas uncertainties in the reconstructions from the individual firn sites are large, Fig. 10c could suggest a relatively narrow range of the CH₄ history that satisfies both NGRIP and NEEM reconstructions. Although the overlapping range of the two reconstruction is ∼ 90 ppb in the 1970s, it is as small as ∼ 30 ppb in the 1950s. This suggests that the combined NGRIP and NEEM firn data could provide a stronger constraint to the range of the CH₄ mole fraction, e.g. 1185–1215 and 1225–1260 ppb in 1950 and 1955, respectively (Fig. 10c inset). It is again noted that only the BZ scenario fall within these ranges for the period, suggesting that it is likely closer to the true atmospheric CH₄ history than the CMIP6 scenario for the mid-20th century.

The IPD of the atmospheric CH₄ mole fraction is important for a better understanding of the evolution of the global CH₄ budget. Given that the Antarctic ice core and firn measurements have provided relatively reliable CH₄ records over the 20th century, improved reconstructions from Greenland ice cores and firn air should better constrain the changes in the IPD. To sufficiently constrain the historical global CH₄ budget, the reconstruction for Greenland needs to be accurate within ∼ 10 ppb, corresponding to ∼ 30 Tg CH₄ yr⁻¹ global emission. Based on the current best firn CH₄ data from NGRIP and NEEM, we demonstrated that consistent reconstruction of the Arctic CH₄ mole fraction is achievable back to the 1950s, but the uncertainty of the reconstruction is still large (> 30 ppb) for the 1950s to 1970s.

5.2 Uncertainty of effective age

It is important to note that the reconstructions for the period before 1980 from the NGRIP firn were heavily influenced by the five deepest data in the LIZ (> 72 m). Figure 11a shows distributions of the effective age of CH₄ at depths below 55 m in the NGRIP firn, which is coloured 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 five lowest sampling depths at NGRIP have effective ages corresponding to the period from ∼ 1950 to the late 1970s. At those depths, even the acceptable diffusivities yield the spread of the effective age of > 5 years (black vertical bars). This shows that the reconstruction of the CH₄ mole fraction for the period is subject to much uncertainty in effective age. For the NEEM firn, the reconstructions before 1980 also rely on the five deepest data in the LIZ (Fig. 11b). The σ_age values at those depths in the NEEM firn range from 5 to 8 years (thick vertical bars), comparable to those in the NGRIP firn. The Antarctic atmospheric CH₄ record (see Figs. 1 and 10) indicates that the atmospheric increase rate of CH₄ was 10–15 ppb yr⁻¹ 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 CH₄ level. This is comparable to the IPD of CH₄ mole fraction during the 1950s to 1970s, which was assumed when Buizert et al. (2012) prepared the BZ CH₄ scenario.

Figure 11 Depth profiles of effective age of CH₄ calculated after the iterative dating calculations for depths below 55 m at the (a) NGRIP and (b) NEEM firn sites (left axis). The solid lines represent cases for different diffusivity profiles and are coloured in the same manner as in the earlier figures. Vertical bars in black indicate the spread of the effective age (σ_age; maximum minus minimum) at the individual sampling depths among the modelling cases in red (right axis). For the NEEM firn, the thin vertical bars correspond to the spread for simulation cases in which ¹⁴CO₂ data are excluded for evaluations of the diffusivity profiles.

It is also interesting to note that σ_age at the four deepest depths in the NEEM firn is almost constant, whereas it increases with depth in the NGRIP firn and exceeds 10 years at the two deepest depths. This indicates that effective age in the oldest firn air layers can be estimated with better accuracy in the NEEM firn than in the NGRIP firn, thereby providing the reconstructions with smaller uncertainties. As σ_age was calculated from the simulation cases with acceptable ranges of diffusivity, its magnitude reflects how tightly the diffusivity profile is constrained at each firn site. In other words, our simulations infer that the diffusivity profile in the NEEM firn can be better constrained than the NGRIP firn.

5.3 Constraint to diffusivity profile

In order to see the degree of constraint to the effective diffusivity from different trace gases, we calculated RMSD values for different combinations of the trace gases for the NEEM firn. The choice of the NEEM firn is due to the availability of a larger number of gas species. It was found that the ¹⁴CO₂ data provide strong constraints for narrowing the acceptable range of diffusivity profiles in the LIZ, as its history has a unique peak in the early 1960s due to nuclear weapons testing (Fig. 2). The evaluated RMSD, with the ¹⁴CO₂ data excluded and the historical CH₄ reconstruction from the simulation cases with RMSD < 1.0, is presented in Fig. 10c (thin horizontal bars in light blue). The corresponding spread of the effective age is also shown in Fig. 11b (thin vertical bars in black). These results show that, for the NEEM reconstruction, the uncertainty of the effective age would be doubled at the two deepest depths if it were not for the ¹⁴CO₂ data. It is also interesting that the range of the CH₄ reconstruction without ¹⁴CO₂ would deviate in the younger ages and would suggest a historical trend lower than the BZ scenario (Fig. 10c). In turn, the CH₄ trend reconstructed from the NGRIP firn might be different if the ¹⁴CO₂ measurement was available. The above contribution of the ¹⁴CO₂ data for the NEEM firn implies that the NGRIP reconstruction could have been closer to the NEEM reconstruction. It is also noted that constraints from halocarbon species to the diffusivity profile are relatively weak, as their mole fractions in the LIZ decreased to close to zero.

5.4 Outlook

Previous studies have used CH₄ as one of tracers that constrain diffusivity profiles in firn (Witrant et al., 2012; Trudinger et al., 2013). These studies have shown that CH₄ effectively contributed to constraining diffusivity for deep layers, i.e. LIZ in the NEEM firn. However, the use of CH₄ as an effective constraint is valid only when its atmospheric scenario given as input to the models was assumed to be correct. Figs. 6 and 8 show a larger spread in the model–data differences of CH₄ among the different diffusivity cases spread widely in the deepest layers of the LIZ of the NGRIP firn, in comparison to the other six gases. This fact highlights that subtle changes of the diffusivity profile in the LIZ have a large impact on simulating CH₄, and it indeed indicates that the species could serve as an effective diffusivity constraint if its atmospheric scenario was correctly given with low uncertainty. This study indicated that the two currently available Greenland firn data sets (NGRIP and NEEM) are more consistent with the BZ CH₄ scenario (Buizert et al., 2012) than the CMIP6 scenario (Meinshausen et al., 2017). Furthermore, it should be again pointed out that the CMIP6 scenario suggests an almost constant IPD of ∼ 130 ppb over the 20th century (Fig. 2). Such constant IPD is unlikely because CH₄ emissions are considered to have increased in the Northern Hemisphere for that period, which requires IPD to increase with time as discussed in the previous studies (Dlugokencky et al., 2003; Ghosh et al., 2015; Chandra et al., 2021). Accordingly, the BZ scenario appears to be a useful choice in firn modelling at Greenland sites, but it should be kept in mind that the use of CH₄ as a tuning tracer could lead to overfitting of the diffusivity profile.

A possibility to improve the reproducibility of the depth profile of CH₄ in the NGRIP firn could come from an additional constraint to the diffusivity profile, along with those currently made by the six trace gases (CO₂, SF₆, CFC-11, CFC-12, CFC-113, and CH₃CCl₃). In particular, it was indicated that ¹⁴CO₂, if available, would have strongly constrained the diffusivity profile and reduced the uncertainty of the historical CH₄ trend reconstructed from the NGRIP firn. This study showed that the currently available firn data from Greenland (NGRIP and NEEM) are in better agreement with the historical CH₄ scenario prepared for the NEEM firn modelling (Buizert et al., 2012) than that for the CMIP6 experiments (Meinshausen et al., 2017). Since the latter scenario relies on the NEEM-S1 ice core data (Rhodes et al., 2013), this study highlighted inconsistency between the ice core and two sets of firn data in Greenland. Given that reconstruction of the CH₄ history from the deepest firn layers is challenging (in terms of the diffusivity versus history problem, as shown in this study), future sampling and measurements of ice cores at a high-accumulation site in Greenland (where the age of the air occluded can be determined accurately) may be the only way to reconstruct the atmospheric CH₄ trend over the 20th century.