Post-bubble close-off fractionation of gases in polar ﬁrn and ice cores: effects of accumulation rate on permeation through overloading pressure

. Gases in ice cores are invaluable archives of past environmental changes (e


Introduction
Atmospheric gases trapped in the firn layer (unconsolidated snow layer; ∼ 70 m at the Greenland summit and preserved in the underlying ice sheets provide precious and continuous records of the past atmosphere and environments (Petit et al., 1999;Spahni et al., 2005;Ahn and Brook, 2008;Kobashi et al., 2008a).However, to reconstruct the original environmental records, it is important to understand the processes of air trapping in the firn and how the air is retained in the ice until it is analyzed in laboratories.There are two well-known processes that change air composition before the air is trapped within bubbles in the firn.First, gravitational fractionation separates gases according to their mass differences and diffusive column height of the firn layer (Craig et al., 1988;Schwander, 1989).Second, a temperature gradient ( T ) between the top and bottom of the firn layer induces thermal fractionation generally pulling heavier gases toward the colder end (Severinghaus et al., 1998).In this study, we investigate a third process that occurs after the bubbles are closed (post-bubble close-off fractionation) and that preferentially affects gases with smaller molecular sizes (< 3.6 Å; for example, helium, neon, oxygen, and argon) but also gases with larger molecular sizes in smaller magnitudes (Ikeda-Fukazawa et al., 2005;Huber et al., 2006;Ikeda-Fukazawa and Kawamura, 2006;Severinghaus and Battle, 2006;Ahn et al., 2008).This fractionation continues deep in ice sheets smoothing signals (Ahn et al., 2008;Bereiter et al., 2014), and the process continues during/after coring (Ikeda-Fukazawa et al., 2005;Kobashi et al., 2008b;Suwa and Bender, 2008b;Bereiter et al., 2009;Vinther et al., 2009).
Clear evidence of the diffusive gas loss from ice cores through ice crystals has been observed in the oxygen content in ice cores as a depletion of oxygen relative to nitrogen (Bender et al., 1995;Ikeda-Fukazawa et al., 2005;Suwa and Bender, 2008b).Depletion of air content by ∼ 10 % was observed for the Camp Century ice core after storage for 35 years, although possible analytical differences between early and late measurements cannot be rejected (Vinther et al., 2009).The process is highly temperature dependent, and the gas loss is induced by the pressure gradients between the bubbles and the atmosphere (Ikeda-Fukazawa et al., 2005).In ice sheets, the concentration gradients at different depths drive the gas diffusion through ice crystals, which smooth climate signals (Bereiter et al., 2014).Firn air studies showed that smaller molecules such as helium, neon, oxygen, and argon preferentially leak out from the closed bubbles, leading to enrichments of these gases in open pores near the bubbleclose-off depth, which leads to depletions of these smaller gases in the closed bubbles (Huber et al., 2006;Severinghaus and Battle, 2006;Battle et al., 2011).However, the mechanisms creating δAr/N 2 or δO 2 /N 2 variations in the time domain (i.e., ice cores) are still poorly understood.
On a longer timescale (i.e, orbital), variations in δO 2 /N 2 closely follow local insolation changes (Bender, 2002;Kawamura et al., 2007;Suwa and Bender, 2008a;Landais et al., 2012).As a possible mechanism, it has been hypothesized that changes in local insolation affect physical properties of the snow at the surface and persist into the bubble close-off depth, controlling the δO 2 /N 2 fractionation (insolation hypothesis) (Bender, 2002;Fujita et al., 2009).In addition, air content in ice cores is also found to co-vary with δO 2 /N 2 on the orbital timescale, indicating common causes (Raynaud et al., 2007;Lipenkov et al., 2011).According to this hypothesis, the orbital signals in δO 2 /N 2 in ice cores are linked to the ice chronology rather than to the gas chronology, which differ by up to a few thousand years.Therefore, the precise understanding of the gas loss process in the firn is essential to determine how climate signals in the bubbles are placed between the ice ages and gas ages on the orbital timescale.
In this paper, encouraged by the observation of a significant negative correlation between δAr/N 2 and the accumulation rate over the past 6000 years in the GISP2 (Greenland Ice Sheet Project 2) ice core (Fig. 1), we investigated the processes of multidecadal to centennial δAr/N 2 variability in three ice cores (GISP2, NGRIP, and Dome Fuji) as well as gas loss processes during storage.δAr/N 2 variations are generally highly correlated with δO 2 /N 2 in ice cores, suggesting similar processes driving these variations (Bender et al., 1995).As δAr/N 2 is nearly constant in the atmosphere over the relevant period (Kobashi et al., 2010), it is better suited to assess the permeation processes in the firn and ice cores than δO 2 /N 2 that varied in the atmosphere by ∼ 1.5 ‰ during the glacial-interglacial cycles (Bender et al., 1995).In the following sections, we first describe the data and investigate the relationships between δAr/N 2 and changes in accumulation rates and surface temperatures.Then, the fractionation processes are examined by applying a permeation model to the ice cores and the firn under two processes: "pressure sensitive processes" (e.g, microbubbles) and "normal bubble process".Finally, we discuss our findings, draw conclusions and mention implications.
2 Data description δAr/N 2 data from three ice cores covering the past millennia (NGRIP, Dome Fuji, and GISP2) were used for the analyses.GISP2 and NGRIP (North Greenland Ice Core Project) data have been published earlier (Kobashi et al., 2008b(Kobashi et al., , 2015)), and Dome Fuji data are new.Importantly, storage histories of these cores (i.e., temperatures) are known and methods for measuring δAr/N 2 are all comparable.GISP2 and NGRIP ice cores were drilled from the Greenland ice sheet, and Dome Fuji was drilled from the Antarctic ice sheet (Table 1).For GISP2, sample resolution varies from 10 to 20 years with high-resolution analyses covering the past 1000 years (Kobashi et al., 2008b(Kobashi et al., , 2010) ) and around the 8.2 ka event (8100 ± 500 years before present Figure 1.δAr/N 2 vs. accumulation rates or air contents in GISP2 over the past 6000 years.Note that δAr/N 2 are corrected values for the post-coring fractionation (1.5 ‰ added).A spline with a 21-year cut-off period (blue line) was applied to the δAr/N 2 data.Two σ error bounds are shown, which were estimated by 1000 iterations of Monte Carlo simulations.Accumulation rates (m ice year −1 ) (black line) were filtered by 21-year RMs.Note that the y axis for the accumulation rate is reversed.δAr/N 2 vs. accumulation rates are significantly negatively correlated over the past 6000 years (r = −0.35,p = 0.03).δAr/N 2 and air contents are significantly positively correlated over the past 6000 years (r = 0.38, p < 0.001 after linear detrending).A slight shift of the air contents around 3700 BP is probably due to the analytical changes that occurred between two different periods of measurements (Kobashi et al., 2008b).The correlations between δAr/N 2 and air contents before and after 3700 BP are similar and significant (r = 0.30, p = 0.002 and r = 0.26, p = 0.008, respectively).Therefore, the δAr/N 2 variation can explain 7-14 % of the total variance of the air contents.
(BP), "present" is defined as 1950) (Kobashi et al., 2007).For NGRIP, sample resolution is about 10 years throughout the past 2100 years (Kobashi et al., 2015).Both GISP2 and NGRIP have similar annual average temperatures of approximately −30 • C (Table 1).However, the accumulation rate of NGRIP (∼ 0.19 m ice year −1 ) is 20 % less than that of GISP2 (0.24 m ice year −1 ) over the past 2100 years and, importantly, its variation (standard deviation after 21-year running means; RMs) is lower by 40 % than that of GISP2 (see later discussion).Dome Fuji has a radically different environment from Greenland with the current annual average air temperature of −54.3 • C and a mean accumulation rate of ∼ 0.03 m ice year −1 (Watanabe et al., 2003).
For the timescale of GISP2 and NGRIP ice ages, we used the GICC05 (Greenland Ice Core Chronology 2005; Vinther et al., 2006;Seierstad et al., 2014).To obtain gas ages, we applied a firn-densification, heat-diffusion model (Goujon et al., 2003) that calculates firn density structure, close-off depth, and delta age.The gas age uncertainties relative to ice age were estimated as ∼ 10 % of the estimated gas age-ice age difference (Goujon et al., 2003).To investigate the δAr/N 2 fractionation, we used reconstructed temperature records from argon and nitrogen isotopes in the trapped air within the GISP2 ice core for the past 4000 years (Kobashi et al., 2011) and within NGRIP for the past 2100 years (Kobashi et al., 2015) and layer-counted accumulation rate data for the Holocene (Alley et al., 1997;Cuffey and Clow, 1997;Gkinis et al., 2014).Dome Fuji data have neither precise temperatures nor accumulation rates over the past 2100 years.The annual resolution accumulation rate data were smoothed with 21-year RMs to correspond to gas diffusion and the bubble close-off process in the firn (Kobashi et al., 2015).A spline fit (Enting, 1987) was applied to gas data (e.g., δAr/N 2 ) with a 21-year cut-off period to be consistent with 21 RMs of other parameters and used for the following analyses to investigate signals longer than the decadal timescale.
GISP2 and NGRIP ice cores were analyzed for δAr/N 2 ∼ 14 years after coring, however, with different temperature histories.GISP2 (82.4-540 m) was drilled in summer 1991.After shipment, they were stored at −29 • C in a commercial freezer until they were moved to a freezer (−36  1993(G. Hargreaves, personal communication, 2015)).The ice samples were then cut and moved to the Scripps Institution of Oceanography, where δAr/N 2 was measured in 2005 (Kobashi et al., 2008b).On the other hand, the NGRIP2 ice cores (one of the two NGRIP ice cores; 64.6 m-445.2m) Table 1.Environmental parameters for GISP2, NGRIP and Dome Fuji.Temperatures for GISP2 and NGRIP are averages over the past 2100 years (Kobashi et al., 2015).Accumulation rates (Alley et al., 1997;Cuffey and Clow, 1997;Gkinis et al., 2014) for GISP2 and NGRIP are averages for the past 2100 years, and accumulation rate variations are calculated as standard deviations of accumulation rates in 21-year RMs.Annual average temperature and accumulation rate for Dome Fuji are from Watanabe et al. (2003).
Latitude Longitude Altitude Average Accumulation Accumulation m a.s.l.temperature rate rate variation ( • C) (m ice year −1 ) (m ice year −1 ) GISP2 72.59  (S. Kipfstuhl, personal communication, 2015).In fall 2011, we cut the ice samples and shipped them to a freezer at the National Institute of Polar Research at −30 • C until 2013 when we analyzed the ice cores (Kobashi et al., 2015).The ice cores from Dome Fuji were drilled in late 1995 and stored at −50 • C, with a short period (2.5 months) at < −25 • C during shipment from Antarctica to Japan (S. Fujita, personal communication, 2015).The ice core was analyzed in early 2014.
The conventional delta notation is used to express δAr/N 2 as follows: where the subscript "sample" indicates ice core values, and "standard" is the present atmospheric composition.For GISP2, mass 40 of argon and 29 of nitrogen, and for NGRIP and Dome Fuji, mass 40 of argon and 28 of nitrogen were used to calculate δAr/N 2 .All δAr/N 2 data presented in this study were corrected for gravitational and thermal fractionations in the firn using the conventional method based on δ 15 N (Bender et al., 1995;Severinghaus and Battle, 2006;Severinghaus et al., 2009) as follows: (2) The coefficient 11 arises because the mass difference of δAr/N 2 ( 40 Ar and 29 N 2 ) is 11 times larger than that of the nitrogen isotopes ( 29 N 2 and 28 N 2 ).This coefficient is replaced with 12 for the calculation of δAr/N 2corr for NGRIP and Dome Fuji because the mass difference between 40 Ar and 28 N 2 is 12.As the temperature sensitivities of δ 15 N and δAr/N 2 are slightly different, the correction is not perfect.However, the variability induced by the gas loss is much bigger than the uncertainties introduced by the differences of the thermal sensitivities.After these corrections, the δAr/N 2corr variations in the ice cores can be attributed only to gas loss.δAr/N 2corr of the GISP2 data using the mass 28 or 29 leads to negligible differences (an average difference is 0.4 × 10 −3 ‰ and the standard deviation is 0.94 × 10 −3 ‰) which are much smaller than the measurement uncertainty of δAr/N 2 (1σ < 0.7 ‰).We also note that the standard deviation (0.07 ‰) of δ 15 N × 11 in GISP2 is much smaller than standard deviation of raw δAr/N 2 (1.33 ‰) over the past 6000 years, indicating that the variations of δAr/N 2corr mostly originate from the raw δAr/N 2 and not from δ 15 N.For the sake of simplicity, we denote all the δAr/N 2corr as δAr/N 2 in later sections.The significance of correlations were calculated considering the autocorrelation of time series (Ito and Minobe, 2010;Kobashi et al., 2013).We consider > 95 % confidence as significant, unless otherwise noted.All error bounds in figures and texts are 2σ .

Post-coring fractionation
Before evaluating δAr/N 2 in ice cores for the changes that have occurred in the firn, it is necessary to consider the postcoring fractionation (Ikeda-Fukazawa et al., 2005).For this purpose, we applied a molecular diffusion model (permeation model) through ice crystals (Ikeda-Fukazawa et al., 2005).It has been applied to observed depletions of oxygen in the Dome Fuji and GISP2 ice cores in ∼ 10 ‰ with respect to nitrogen (Ikeda-Fukazawa et al., 2005;Suwa and Bender, 2008b).The model was also implemented with modifications for gas permeation processes in the firn (Severinghaus and Battle, 2006) and in ice cores (Bereiter et al., 2009).The gas permeation in ice cores is driven by the pressure gradients between two spaces isolated by ice walls (e.g., between bubbles or between bubbles and the atmosphere).The concentration (U m ; mol mol −1 ice ) of species m (i.e., nitrogen, oxygen, and argon) in bubbles in 1 mol of ice after a time t can be described as follows (Ikeda-Fukazawa et al., 2005): where U 0 m (mol mol −1 ice ) is the original concentration of species m. k m (m s −1 ) is the mass transfer coefficient, and it is equal to D m / l, where D m (m 2 s −1 ) is the diffusion coefficient of the species m and l (m) is the thickness of the surface layer of ice (Ikeda-Fukazawa et al., 2005).X m (mol mol −1 ice M Pa −1 ) is the solubility of species m in ice.P i and P a are the pressures in the bubbles and in the atmosphere, respectively.Z i m and Z s m are molar fractions of species m in the bubbles and in the atmosphere, respectively.S (m 2 ) and V (m 3 ) represent the surface area and the volume of an ice sample such that S/V can be understood as specific surface area (m −1 ), an important variable for the gas exchange between the atmosphere and the ice (Matzl and Schneebeli, 2006).
For Eq. ( 3), we assumed an initial air content of 6.53 × 10 −5 mol in 1 mol of ice (a typical air content in ice cores).U 0 m for each gas is calculated from the total air content multiplied by the atmospheric molar ratio of each gas.In this case, Z i m and Z s m are set to the atmospheric partial pressures for each molecule.Another factor that affects the gas loss is the specific surface area.The GISP2 ice core has a larger diameter (0.132 m) and longer length (1 m) during the storage than the NGRIP core (0.098 m diameter and 0.55 m length).Dome Fuji core has a diameter of 0.093 m and length of 0.50 m.Therefore, the specific surface areas (S/V ) were calculated to be 32.3, 44.5, and 47.0 m −1 for GISP2, NGRIP, and Dome Fuji, respectively.It is noted that these specific surface areas are approximations given that ice cores during storage often have different shapes from earlier sampling shapes and because we shaved the ice surface by ∼ 5 mm before the analyses (Kobashi et al., 2008b(Kobashi et al., , 2015)).However, we also note that shallow late Holocene ice cores often had nearly intact shapes (no sampling) at the time of our sampling from ice cores.
Diffusivity (D Ar ) and solubility (X Ar ) for argon in ice are less known than for nitrogen and oxygen.Therefore, we attempted to estimate two possible functions (Ar(I) and Ar(II)) for k Ar X Ar (i.e., D Ar / l • X Ar ) in relation to those for nitrogen and oxygen (Fig. 2).K N 2 X N 2 and k O 2 X O 2 in different temperatures can be estimated using Eqs.( 4) and (2) with l = 12 and 7 mm for nitrogen and oxygen in Ikeda-Fukazawa et al. (2005) for the Dome Fuji core (Fig. 2), which were consistent with various observations (Ikeda-Fukazawa et al., 2005;Severinghaus and Battle, 2006;Suwa and Bender, 2008b;Bereiter et al., 2009).
The first estimate of Ar(I) uses a diffusion coefficient (D Ar ; 4.0 × 10 −11 m 2 s −1 ) of argon at 270 K calculated from molecular dynamic simulations with those of nitrogen (D N 2 ; 2.1 × 10 −11 m 2 s −1 ) and oxygen (D O 2 ; 4.7 × 10 −11 m 2 s −1 ) (Ikeda-Fukazawa et al., 2004).Owing to the molecular-sizedependent fractionation, argon permeation occurs slower than oxygen but faster than nitrogen (Fig. 2), which cannot be explained by their mass differences (Huber et al., 2006;Severinghaus and Battle, 2006).Then, temperature-dependent k Ar and X Ar were estimated assuming that the geometrical relationship between D N 2 , D Ar , and D O 2 at 270 K from the molecular dynamics simulations holds for k Ar and X Ar at different temperatures as follows: Second, we estimated Ar(II) from an observation that δAr/N 2 in ice is often depleted about half of δO 2 /N 2 in ice cores (Bender et al., 1995).To satisfy this condition, k Ar X Ar can be written as Estimated k Ar X Ar for Ar(I) and Ar(II) are higher than k N 2 X N2 and increase with temperatures, resulting in a general depletion of δAr/N 2 in ice compared to the atmospheric composition, and the depletion is faster in warmer temperatures (Fig. 2).The use of Ar(I) induces faster depletion of δAr/N 2 than that of Ar(II) owing to faster permeation of argon.With the two estimates of k Ar X Ar , we explore the range of uncertainties associated with argon permeation.
In a pioneering study by Bender et al. (1995), δAr/N 2 in a shallow core of GISP2 was analyzed after 1 week, 3 months, and 7 months of drilling in 1989 to study the timedependent gas loss process (Fig. 3).As the data from three different periods are not significantly different, we treat Bender's data as the true Ar/N 2 values before coring.By comparing the data (Bender et al., 1995) with our data set analyzed 14 years after the coring (Kobashi et al., 2008b), we estimated the post-coring fractionation of δAr/N 2 in GISP2 to be −1.5 ± 0.6 ‰, a difference of the two data sets for common depths (124-214 m) (Fig. 3).Using this value, we derived an unknown parameter (i.e., bubble pressure) in Eq. ( 3).The bubble pressures are calculated as 0.6 ± 0.    (Lipenkov, 2000).
4 Post-bubble close-off fractionation in firn: empirical evidence

GISP2 δAr/N 2 variation over the Holocene
The δAr/N 2 record over the Holocene in the GISP2 ice core exhibits relatively constant values of around −3 ‰, except for a prominent rise of up to 10 ‰ around 7000 years BP (Fig. 4).The rise is located within the depths of the brittle zone (650-1400 m), where air in the bubbles changes to clathrate inducing anomalously high pressure (Gow et al., 1997).The dissociation pressure of nitrogen in the clathrate phase is higher than that of argon (or oxygen) so that nitrogen is enriched in the gas phase in relation to the clathrate (more stable state), resulting in a preferential leakage of nitrogen and thus argon (or oxygen) enrichments in these depths (Ikeda et al., 1999;Ikeda-Fukazawa et al., 2001;Kobashi et al., 2008b).As the dissociation of gases from the clathrate depends on various factors, δAr/N 2 in these depths are highly variable (Fig. 4).It is noted that δ 15 N and δ 40 Ar exhibit little influences from the anomalous δAr/N 2 fractionation, indicating that the processes are mass independent  (1995), and black data points with error bounds are from Kobashi et al. (2008b).We did not use shallower data of Bender et al. (1995) as they exhibit depletions similar to our shallow NGRIP data (Fig. 8); an anomalous value (−16.91 ‰ at 145.4 m) in the Bender data set was also excluded.Squares, diamonds, and triangles represent the data measured after 1 week, 3 months, and 7 months of coring, respectively (Bender et al., 1995).
The average difference between the Kobashi and Bender data sets is −1.51 ± 0.58 ‰, which we interpret as the post-coring fractionation for GISP2.
in first order (Huber et al., 2006;Severinghaus and Battle, 2006) (Fig. 4).Changes in the surface temperatures and accumulation rates are the dominant controlling factors for the state of firn layers (e.g., density profile, bubble close-off depth, and firn thickness) (Herron and Langway, 1980;Schwander et al., 1997;Goujon et al., 2003).Therefore, we investigated if changes in surface temperature or accumulation rate have any controls on the δAr/N 2 variations.We found a significant  (Kobashi et al., 2008b).The grey arrow indicates the brittle zone (Gow et al., 1997).
negative correlation (r = −0.35,p = 0.03) between δAr/N 2 on the gas-age scale and the accumulation rate for the past 6000 years, a time interval in which the abnormal δAr/N 2 fractionation is not observed (Figs. 1, 4).This negative correlation is opposite of what an earlier study (Severinghaus and Battle, 2006) suggested for the permeation fractionation in the firn (positive correlation).In addition, the significant correlation was found for δAr/N 2 on the "gas ages" scale rather than the "ice ages" that the insolation hypothesis predicts; an indication that new processes need to be considered for the gas loss processes in the firn.
The subset of the GISP2 data covering the past 4000 years provides a unique opportunity to investigate δAr/N 2 variations because precise temperature (Kobashi et al., 2011) and accumulation rate by layer counting (Alley et al., 1997;Cuffey and Clow, 1997) are available.Using these data, we applied a linear regression and lag analysis on δAr/N 2 .It is found that the surface temperature is positively correlated with δAr/N 2 on the gas ages (r = 0.47, p = 0.04; r = 0.28, p = 0.001 after linear detrending) with a 68-year lag (Fig. 5a), indicating that cooler (warmer) temperatures induced more (less) depletions in δAr/N 2 with a multidecadal lag.On the other hand, the accumulation rate is negatively correlated with δAr/N 2 on the gas ages (r = −0.47,p = 0.12; r = −0.26,p = 0.01 after linear detrending) with a 38-year lag (Fig. 5b), indicating that high (low) accumulation rates induced more (less) depletions in δAr/N 2 over the past 4000 years.We note that the surface temperature and accumulation rate have a negative but insignificant correlation (r = −0.32,p = 0.13; after linear detrending r = −0.11,p = 0.2) over the past 4000 years.
To estimate the relative contribution of the accumulation rate and the surface temperature changes on δAr/N 2 , we applied a multiple linear regression, which finds the best linear combination of variables (i.e., temperature and accumulation rate) for a response variable (i.e., δAr/N 2 ).Before the regression is applied, the temperature and accumulation records were shifted toward younger ages to account for the lags (38 and 68 years for accumulation rate and temperature, respectively), and δAr/N 2 is corrected for the post-coring fractionation (1.5 ‰ added).As ordinary least squares including the multiple linear regression underestimate the variance of target time series when the data are noisy (Von Storch et al., 2004), we used "variance matching" by linearly scaling regression coefficients according to the ratio between the variance of the target and model time series.Figure 5c shows the original and modeled results of δAr/N 2 over the past 4000 years.As expected, the model of the multiple linear regression captures the δAr/N 2 variations better than the individual variables do (Fig. 5a-c) with a correlation coefficient of r = 0.58, p = 0.09 (r = 0.36, p < 0.001 after linear detrending).For the centennial variations, the model captures nearly half of the total variance of the observed δAr/N 2 variations with a 95 % confidence (r = 0.71, p = 0.05 after linear detrending with 200-year RMs).The high and significant correlation between the model and observed δAr/N 2 indicates that changes in the surface temperature and accumulation rate play important roles in controlling the δAr/N 2 variations.From the multiple linear regression, δAr/N 2 on Figure 5. Observed and modeled δAr/N 2 , surface temperatures, and accumulation rates from the GISP2 ice core over the past 4000 years.Note that the observed δAr/N 2 is corrected for the post-coring fractionation (1.5 ‰ added).(a) δAr/N 2 and surface temperatures (Kobashi et al., 2011).Ages of the temperatures were adjusted for the lag (68 years).(b) δAr/N 2 and accumulation rates in 21-year RMs (Alley et al., 1997;Cuffey and Clow, 1997).Ages of the accumulation rates were adjusted for the lag (38 years) (c) Observed and modeled δAr/N 2 from multiple linear regression (see text).(d) Observed and modeled δAr/N 2 of multiple linear regression using δ 18 O ice as a temperature proxy (see text).
the gas ages in GISP2 can be expressed by temperature ( • C) and accumulation rate (m ice year −1 ) as a function of time after adjusting for the lags: and t, t temp , and t accm are time, lag for temperature, and lag for accumulation rate, respectively (all in years).
Next, we attempted to use oxygen isotopes of ice (δ 18 O ice ) as a temperature proxy for the same regression analyses of δAr/N 2 since we do not have the N 2 /Ar isotope-based temperature information before 4000 years BP.Temperature records derived from δ 18 O ice can be quite noisy but stacking several δ 18 O ice records can improve the derived temperature histories (White et al., 1997;Kobashi et al., 2011).Thus, we stacked three oxygen isotope records (GISP2, GRIP, and NGRIP) over the Holocene in the 20-year RMs (Stuiver et al., 1995;Vinther et al., 2006).The stacked record was calibrated to temperatures using the relation obtained from borehole temperature profiles (Cuffey and Clow, 1997).Using the regression coefficients obtained in Fig. 5c, a δAr/N 2 model was calculated from the oxygen-isotope-based temperature and the accumulation rate (Fig. 5d).We found that the correlation between the model and the observed δAr/N 2 performs not as well as the one with the temperature and accumulation rate records for the last 4000 years based on Ar/N 2 isotope values (Fig. 5c) but does slightly better than the correlations with the temperature or accumulation rate individually (Fig. 5a, b).
The δAr/N 2 regression model with the δ 18 O ice -based temperatures and accumulation rates can span the entire Holocene, including the periods when the observed δAr/N 2 are highly variable owing to the post-coring fractionation as discussed earlier.Except for the time interval around 7000 years BP, the model and observed δAr/N 2 exhibit rather constant values of −1 to −3 ‰ during the Holocene (Fig. 6).Interestingly, the model indicates that the constant δAr/N 2 during the early Holocene is the result of a cancellation between the effects of the accumulation rate and the temperature, both of which were rapidly rising in the early Holocene (Fig. 6).The observed δAr/N 2 variations remained higher or noisier from the early Holocene to ∼ 6000 BP than for the later period, which probably make it difficult to decipher the original multidecadal to centennial signals in δAr/N 2 (Fig. 6).

NGRIP and Dome Fuji δAr/N 2 variation over the past 2100 years
The δAr/N 2 data of NGRIP ice cores allows for a good comparison with the GISP2 data (Fig. 7).The averages of δAr/N 2 for the past 2100 years are −3.36 and −2.40 ‰ for NGRIP and GISP2, respectively (Fig. 7).The δAr/N 2 variability in NGRIP (1σ = 0.91 ‰) over the past 2100 years is 24 % smaller than that of GISP2 (1σ = 1.19 ‰) after correcting for the post-coring fractionation (Table 2), likely owing to the smaller variations of the accumulation rate at NGRIP than that of GISP2 (Fig. 7).The pooled standard deviations of replicated samples are 0.94 ‰ for NGRIP over the past 2100 years and 0.66 ‰ for GISP2 over the past 1000 years (replicates are available only for the past 1000 years in GISP2) (Kobashi et al., 2008b).The noisier data for NGRIP than for GISP2 should not be analytical as the mass spectrometer used for the NGRIP had better precision on δAr/N 2 than the one used for GISP2 (Kobashi et al., 2008b(Kobashi et al., , 2015)).δAr/N 2 for GISP2 and NGRIP are weakly but significantly correlated with a correlation coefficient of r = 0.24 and p = 0.02 (after linear detrending) for the past 1000 years of the high-resolution part of GISP2 but not for the deeper part, likely owing to the difference of sampling densities between the two periods (Kobashi et al., 2015).The surface temperatures at NGRIP were only weakly correlated with δAr/N 2 in the deeper part of NGRIP (r = 0.20, p = 0.06 after linear detrending) and were uncorrelated in the shallower part.
No significant correlations were found between δAr/N 2 and the accumulation rate for NGRIP, probably due to the lower variation of the accumulation rate at NGRIP than that of GISP2.It is consistent with the fact that the signal to noise ratio (SNR = variance of signals/variance of analytical er-rors = 1.2) for NGRIP is about one-fifth of that for GISP2 (6.1) estimating the NGRIP signals from Eq. ( 7).
From the relationship between δAr/N 2 and the temperature or accumulation rate of GISP2 in Eq. ( 7), we can calculate expected δAr/N 2 for NGRIP and Dome Fuji.Using the past 2100 years of temperatures and accumulation rates for NGRIP (Fig. 7a, b) and the current observation (Table 1) for Dome Fuji, the expected δAr/N 2 from Eq. ( 7) were calculated as 0.3 ± 1.3 and −6.4 ± 1.2 ‰, respectively.The value for NGRIP is significantly higher than the observed value of −3.3 ± 1.2 ‰ corrected for the post-coring fractionation (Table 2).For Dome Fuji, the value is similar to the observed −6.3 ± 0.8 ‰ corrected for the post-coring fractionation (Fig. 7, Table 2).This may indicate that the relationship of δAr/N 2 with the temperature and accumulation rate becomes non-linear when the firn thickness becomes thinner than that of GISP2 as δAr/N 2 is not expected to be positive without the existence of clathrate (see later discussion).
δAr/N 2 ice core data of NGRIP from the depth range 64.6-80 m exhibit some interesting features (Fig. 8).The depth from ∼ 60 to 78 m corresponds to the lock-in zone in NGRIP, where vertical mixing of gas is limited so that δ 15 N stays nearly constant in these depths (Huber et al., 2006;Kawamura et al., 2006).Therefore, the shallowest data at 64.6 m are located in the lock-in zone.Generally, gas data from the lock-in zone are not used owing to possible contamination (Aydin et al., 2010).However, a recent study (Mitchell et al., 2015) demonstrated that δ 15 N can be used to estimate the amount of ambient air contamination using ice samples in the lock-in zone, and the original methane concentration in the firn was reconstructed with a range of uncertainties.Therefore, we interpret the observed rapid decreases of δ 15 N and δ 40 Ar toward shallower depths in the lock-inzone as the result of mixing with ambient air (Fig. 8d).Based on isotope mass balance, we calculated the original δAr/N 2 values, which exhibited highly depleted values as low as −50 ‰ (Fig. 8e).The depleted δAr/N 2 in the lock-in zone provides a clue to the processes of gas loss in the firn (see later discussion).

Post-bubble close-off fractionation in firn:
process study Air bubbles in the polar firn or ice can be categorized into two types (Lipenkov, 2000): normal bubbles and microbubbles (< 50 µm).They can be distinguished as a bimodal distribution in ice cores (Lipenkov, 2000;Ueltzhöffer et al., 2010;Bendel et al., 2013).The air volume contribution of the microbubbles to the total air content is estimated to be 0.3 % in the Vostok ice core (Lipenkov, 2000), but the value is not known for Greenland ice cores.Importantly, the two types of bubbles have significantly different bubble pressure histories in the firn.The normal bubbles form at the bubble close-off depth.Most of the air in ice cores is captured as normal bub- bles, and the air-trapping processes are relatively well known (Schwander et al., 1997;Goujon et al., 2003;Mitchell et al., 2015).Normal bubble pressures build up according to increasing density (normal bubble process; Severinghaus and Battle, 2006).On the other hand, the microbubbles are believed to form near the surface (Lipenkov, 2000); so, they are highly pressurized and have a rounded shape by the time when the bubbles reach the bubble close-off depth (Lipenkov, 2000;Ueltzhöffer et al., 2010).As a result, the microbubbles are more sensitive to changes in the overloading pressure at the bubble close-off depth (pressure sensitive process).
Owing to the different bubble pressure histories in the firn, δAr/N 2 or δO 2 /N 2 in the microbubbles and normal bubbles are expected to be different due to the differential permeation of each molecule.In this study, we attempted to quantify two types of the gas loss processes -pressure sensitive process (microbubble) and normal bubble process -in the firn using a permeation model (Ikeda-Fukazawa et al., 2005) combined with the inputs from firn-densification, heat-diffusion models (Schwander et al., 1993;Spahni et al., 2003;Goujon et al., 2003).

Pressure sensitive process (microbubbles)
We first look into the pressure sensitive process as exemplified by the microbubbles.Microbubbles are believed to form in the shallow firn by sublimation-condensation processes (Lipenkov, 2000).These bubbles have smaller sizes, smoothed spherical surfaces, and can generally be found in the interior of the ice crystals (Lipenkov, 2000).The bubble pressure reaches a near-overloading pressure at the bubble close-off depth, so it is sensitive to changes in the overloading pressure.As the actual contribution of microbubbles and air content involved in the pressure sensitive processes is not known, we consider a 2 % contribution of air to the total air.As it will be discussed later, more air fractions than simply from microbubbles (0.3 % in Vostok) are likely involved in the pressure sensitive process.Therefore, we conducted additional calculations with 0.3, 1, and 3 % microbubble contributions, and assessed the impacts to the total δAr/N 2 .
To model the gas permeation process from the microbubbles, we assumed a steady state with given surface temperatures and accumulation rates and calculated the ages, firn densities, porosities, and overloading pressures at given depths, using a firn-densification, heat-diffusion model (Schwander et al., 1993;Spahni et al., 2003).Then, they are interpolated for annual layers in the firn for the following calculation.
Changes in the concentrations of species m were calculated according to Eq. ( 8), similar to Eq. (3).
where l is an annual layer from the surface to below the firn layer (e.g.time variable.At l = 1, the microbubbles in an annual layer are at the surface, although they are not active in terms of permeation at these depths (Fig. 9).With l increasing in a 1-year step, the microbubbles move deeper in the firn with l annual layers overlying.C(l) is a coefficient defining the gas concentration in annual layer l relative to the total air in ice.It is assumed that the pressure P i (l) in the microbubbles starts increasing with overloading pressure from the depth at which the normal bubbles generation initiates (firn density of around 0.7 g cm −3 ) (Fig. 9c) and that pressure changes were considered to be negligible above that depth (Lipenkov, 2000).Initial P i (0) was set at 0.065 MPa, similar to the atmospheric pressure at the Greenland summit (Schwander et al., 1993) with a 0.3 MPa lag from overloading pressure as in Fig. 9 (Lipenkov, 2000).We estimated the specific surface area (S/V (l)) in a layer l from the linear relationship between the specific surface areas (m −1 ) and densities ρ from the Greenland summit (Lomonaco et al., 2011) with an equation: S/V (l) (m −1 ) = −16 799 ρ(l) (g cm −3 ) + 14 957.The initial gas content in the microbubbles was set at 0.3- 3 % of the air content (6.53226 × 10 −5 mol × 0.01) per 1 mol of ice, and it is composed of nitrogen (78.084 %), oxygen (20.9476 %), and argon (0.934 %).The specific surface area S/V was multiplied by the open porosity ratio s o /s(l) (Spahni et al., 2003;Fig. 9a), as the gas loss occurs toward open pores.k m X m was calculated as for the post-coring fractionation, and we used the estimate Ar(II) for argon.
Figure 9 shows model results with a temperature of −31 • C, an accumulation rate of 0.25 m ice year −1 (similar to GISP2 condition), and 2 % microbubble contribution.It shows that the gas permeation from the microbubbles starts soon after the pressure was applied in the microbubbles (Fig. 9c, d).As oxygen has a larger permeability than of argon, δO 2 /N 2 depletion is larger than δAr/N 2 (Fig. 9b).At the temperature of −30 • C and accumulation rate of 0.25 m ice year −1 , the depletion reaches up to 133 ‰ for δAr/N 2 , and 243 ‰ for δO 2 /N 2 in the model, which corresponds to a 12 % gas loss from the original air content of the microbubbles (Fig. 9c).

Normal bubble process
Most of the air in ice cores is trapped as normal bubbles near the lock-in depth.As a result, bulk air pressure in the normal bubbles does not build up as high as the microbubbles in the lock-in zone (Lipenkov, 2000).We used Eq. ( 8) to model the permeation process for the normal bubbles.As for the microbubbles, we assumed a steady state with the given temperatures and accumulation rates.The general characters of the firn in various depths (ages, densities, porosities, loading pressures, bubble close-off depths) were calculated using the firn-densification, heat-diffusion model (Schwander et al., 1993;Spahni et al., 2003), and they were interpolated for annual layers as for the microbubbles.We first calculated how much bubble air is generated in each annual layer according to the increase in the closed porosity (s c ) with depth as the following equation: where V 0 (l) is newly trapped air in an annual layer l, ρ ice is the density of ice, and ρ(l) is the density at depth l. s c (l) is the closed porosity in an annual layer l, and a is a scaling coefficient.s c can be written as (Schwander, 1989;Spahni et al., 2003) where ρ co is the density at the depth in which the air is totally enclosed in bubbles.The sum of all the newly generated air ) is set to have the air content of 6.53 × 10 −5 mol per mole of ice.Then, V 0 (l) was scaled accordingly using the coefficient a and converted to the volume (m 3 ) with the atmospheric pressure (0.065 MPa) as in Fig. 10a.The normal bubbles start forming at approximately 40 m depths and the formation is maximum around the bubble close-off depth of 60-75 m at −31 • C and 0.24 m ice year −1 in the model (Fig. 10a).Then, the permeation from each annual layer was calculated according to Eq. ( 8).The difference from the microbubble permeation process is that the volume of the normal bubbles decreases according to increasing modeled density towards deeper depth, leading to a generally smaller pressure build-up and total permeation from the bubbles in the firn than that of the pressure sensitive process for the microbubbles (Fig. 10a).C(l) in Eq. ( 8) was calculated from V 0 (l) for each annual layer l by setting the sum of C(l) as 1 − 0.02 = 0.98 (if microbubble contribution is 2 %) (Fig. 10e).Other parameters in Eq. ( 8) were set to be the same as for the microbubbles.
Figure 10 shows the evolution of the normal bubble volumes, the nitrogen and argon concentrations, the δAr/N 2 in each annual layer, and the air content and bulk δAr/N 2 with depth at a temperature of −31 • C and an accumulation rate of 0.24 m ice year −1 as for the microbubbles for Fig. 9.A new generation of closed pore volumes in annual layers generally increases towards deeper depths (Fig. 10a).When open pore space disappears completely, we assume the gas permeation to the open pore stops.As argon (oxygen) permeation in ice is faster than nitrogen by 289 (479) % at −31 • C (Ar(II), Fig. 2), δAr/N 2 (δO 2 /N 2 ) within the bubbles decreases when the permeation proceeds.At the temperature of −31 • C and accumulation rate of 0.24 m ice year −1 , the δAr/N 2 depletion can reach about −5 ‰ for those bubbles formed at shallow depths (Fig. 10d).However, the amount of air contained in these bubbles is so small (Fig. 10a) that the influence on the total δAr/N 2 is limited (Fig. 10e vs. δAr/N 2 relationship of the total air from the normal bubbles (Fig. 10e) indicates that the total δAr/N 2 reaches the minimum of −0.39 ‰ at the middle of the bubble closeoff depth of 73.2 m.Then, the total δAr/N 2 increases to −0.29 ‰ as a large amount of ambient air with δAr/N 2 = 0 is trapped in these depths (Fig. 10a, d, e).

Total air in bubbles
The permeation models for the normal and microbubbles were run for various firn conditions with different surface temperatures, accumulation rates, and microbubble contributions to investigate their effects on the δAr/N 2 in the bubbles (Fig. 11, Table 3).The resultant air content (i.e., nitrogen, argon, and oxygen) for each annual layer from the microbubbles and normal bubbles was added to calculate the combined effects of the accumulation rates and temperatures on total δAr/N 2 (Fig. 11).Results show that the normal bubbles experience only limited δAr/N 2 depletion (> −0.5 ‰) by the different temperatures or accumulation rates we considered (Table 3).On the other hand, δAr/N 2 in the microbubbles varies with temperatures through thickening of the firn, leading to higher pressures in the bubbles and longer exposure to the gas loss in the firn (Table 3).A higher accumulation rate with the same temperatures induces more depletion as it is primarily controlled by the changes in loading pressure (Fig. 11c, Table 3).As a result, the total δAr/N 2 generally reflects the variation of δAr/N 2 in the microbubbles (r = 0.95, Table 3).Overall, the total δAr/N 2 have a higher correlation with temperatures (r = 0.97) than with accumulation rates (r = 0.57) in the model (Table 3).The modeled δAr/N 2 agrees with the observed δAr/N 2 corrected for the post-coring fractionation and is within their uncertainty ranges (Table 3).The extremely cold temperature in Dome Fuji with low accumulation rate induces a long (274 years) bubble exposure to the permeation in the firn, leading to a large depletion of δAr/N 2 in the microbubbles and in the total air (Table 3).The variations of δAr/N 2 in normal bubbles are limited and, clearly, microbubbles (or the pressure sensitive process) play a critical role for the variation of δAr/N 2 in ice cores.The δAr/N 2 minima in the firn ranges from −14 to −83 ‰ depending on the temperatures and accumulation rates.The most depleted δAr/N 2 with a temperature of −30 • C and accumulation rate of 0.2 m ice year −1 in Fig. 11c capture the highly depleted observation-based estimates of δAr/N 2 in the NGRIP ice core (Fig. 8e).As the normal bubble process alone does not produce such depleted values in the firn (Fig. 11a), the observed highly depleted δAr/N 2 (Fig. 8e) is an evidence for the involvement of the microbubble process (or pressure sensitive process).The total δAr/N 2 at the bubble close-off depth increase to less depleted values from the minimum owing to the rapid inclusion of the ambient air (Fig. 11c).
The calculated dependencies of the δAr/N 2 variations on the temperature (0.24 ‰ • C −1 for an accumu- lation rate of 0.25 m ice year −1 ) and accumulation rate (−0.05 ‰ (0.01 m ice year −1 ) −1 at −30 • C) with the 2 % microbubble contribution (Table 3) are lower than those of the observed ones in GISP2 ice cores (0.72 ± 0.1 ‰ • C −1 and −0.58 ± 0.09 ‰ (0.01 m ice year −1 ) −1 ), respectively.Considering a possibility of larger volume contributions on the pressure sensitive process, we calculated the permeation model with microbubble volume contributions of 0.3-3 % to the total air.The 3 % microbubble contribution induces more depletion in the total δAr/N 2 (Fig. 11d).Also, the dependencies of δAr/N 2 on temperatures and accumulation rates lin- early increase to 0.38 • C −1 with an accumulation rate of 0.25 m ice year −1 , and −0.11 ‰ (0.01 m ice year −1 ) −1 with a temperature at −30 • C, respectively.The fact that they are still lower than those of the observations indicates the involvement of larger air content as microbubbles and/or normal bubbles influenced by the pressure sensitive process.This is plausible considering the inhomogeneity of firn (Hörhold et al., 2012) and resultant differential pressurization of bubbles at the same depth.
Evidence of the larger air involvement in the pressure sensitive process is the significantly positive correlation between δAr/N 2 and air contents over the past 6000 years in GIPS2 (Fig. 1).This correlation indicates that the bubble air was squeezed out before close-off, resulting in smaller air contents when overloading pressure was higher, eventually inducing higher pressure in the bubbles and thus enhanced δAr/N 2 depletions.This observation is also consistent with recent findings that abrupt increases of accumulation rate at abrupt warming during the last glacial period induced reductions in air contents (Eicher et al., 2015).In addition, artificial sintering of snow with higher pressure has been shown to contain much smaller air content than ice cores, owing to the lack of time to develop spherical cavities by vapor transport (B. Stauffer, personal communication, 2015).These lines of evidence indicate that higher overloading pressures at the lock-in zone have impacts on normal bubbles and microbubbles.The inclusion of this process in the model is beyond the scope of the current paper, and we leave it for future studies.
We also investigated the observed lags of 68 and 38 years, respectively, for the δAr/N 2 variations in GISP2 from the changes in the surface temperatures and accumulation rates (Fig. 5).Presumably, the lags are introduced during the process of transferring surface temperature and accumulation rate signals into overloading pressure at the bubble close-off depths.Therefore, two transient simulations were conducted using a firn-densification and heat-diffusion model (Goujon et al., 2003).First, the model was run with a constant temperature (−30 • C) and accumulation rate (0.2 m ice year −1 ) over thousands of years to reach an equilibrium state.Then, surface temperature and accumulation rate anomalies of −35 • C and 0.26 m ice year −1 for 20 years were introduced, separately (Fig. 12a).The surface anomalies of the temperature and accumulation rate were set to induce similar δAr/N 2 changes of 3.5 ‰ from the relationship obtained by the multiple linear regressions on the δAr/N 2 of GISP2.
We found that the surface temperature anomaly takes 20 years to reach the minimum temperature at the bubble close-off depth (Fig. 12b).The cooling induces maximum firn thickening after 56 years.The accumulation rate anomaly also induces firn thickening with an 11-year lag (Fig. 12c).Overloading pressures at the bubble close-off depth reach similar maximum values with 85-and 21-year lags from the surface temperature and accumulation rate anomalies, respectively (Fig. 12d).Apparently, the surface temperature anomaly takes longer to reach the maximum increase in the overloading pressure than the accumulation rate anomaly, which is consistent with the observations (68 and 38 years, respectively).The accumulation rate anomaly is almost instantaneously and increasingly felt by the bubble close-off depth through overloading pressure, as compared to the temperature anomaly that takes decades to reach the bubble close-off depth.In addition, we note that similar magnitudes of the overloading pressure anomalies were induced by the temperature and accumulation rate anomalies (Fig. 12d).Therefore, we conclude that the overloading pressure is the carrier of the surface temperature and accumulation rate signals, linking the δAr/N 2 variations through the permeation.

Discussions
The processes responsible for the δAr/N 2 variations should also play similar roles on the variations of δO 2 /N 2 in ice cores but with larger magnitudes owing to the larger permeability of oxygen (Bender et al., 1995;Huber et al., 2006;Severinghaus and Battle, 2006;Battle et al., 2011).In ear-lier studies, causes of the δO 2 /N 2 variation were attributed to the metamorphisms of surface snow induced by local insolation changes (Bender, 2002;Kawamura et al., 2007).The altered snow properties remain until the snow reaches the bubble close-off depth and affects the preferential oxygen loss (Bender, 2002).Our work demonstrates that the permeation processes in the firn can be induced by changes in the surface temperature and the accumulation rate through the changes in overloading pressure, indicating a possibility that the δO 2 /N 2 variations in the orbital scale are also a result of the surface temperature and accumulation rate changes.We note that δAr/N 2 in GISP2 also shows a significant positive correlation (r = 0.38, p < 0.001 after linear detrending) with the air content (Kobashi et al., 2008b) over the past 6000 years, indicating a similar link between δO 2 /N 2 and air content in the orbital timescale (Raynaud et al., 2007;Lipenkov et al., 2011).As the timescale we considered in this study is different from the orbital scale variation, other mechanisms may play a role in controlling the δO 2 /N 2 variations in ice cores.However, the mechanisms discussed here must be considered in future studies.
Although the gas permeation from ice is generally believed to be a mass-independent process (no effects on isotopes), there is some evidence of isotopic fractionation (Bender et al., 1995;Severinghaus et al., 2003Severinghaus et al., , 2009;;Severinghaus and Battle, 2006;Kobashi et al., 2008b;Battle et al., 2011).In particular, poor quality ice cores often exhibit isotope fractionation (e.g, δ 18 O and δ 40 Ar) with highly depleted δO 2 /N 2 or δAr/N 2 (Bender et al., 1995;Severinghaus et al., 2009).This mass-dependent fractionation is explained by the existence of microcracks in poor quality ice samples that permit a relatively large air flow.On the other hand, slowly occurring gas permeations through ice crystals in good quality ice cores (e.g, NGRIP, GISP2, and Dome Fuji) appear to have small or non-existent effects on isotopes (Kobashi et al., 2008b;Suwa and Bender, 2008b).As small mass-dependent fractionation of δ 15 N and δ 40 Ar during the gas loss are similar to the gravitational fractionation (Kobashi et al., 2008b), the removal of the gravitational components also cancels the post-coring isotopic fractionation.As a result, the estimated temperature gradients in the firn are little affected by the gas loss (Kobashi et al., 2008b).
Another sign of isotopic fractionation during the gas loss is δ 40 Ar enrichment in ice cores, which produces calculated temperature gradients in the firn to be lower than expected from firn modeling (Kobashi et al., 2010(Kobashi et al., , 2011(Kobashi et al., , 2015)).The systematically higher δ 40 Ar is believed to be caused by processes during the bubble close-off; however, so far no clear evidence has been found in firn air studies (Huber et al., 2006;Severinghaus and Battle, 2006) except δ 18 O of O 2 (Battle et al., 2011).If the enrichment of δ 40 Ar occurs in the firn, it should be correlated with δAr/N 2 .Therefore, the corrections for the δ 40 Ar enrichment have been applied using δAr/N 2 (Kobashi et al., 2010(Kobashi et al., , 2011(Kobashi et al., , 2015)), δKr/Ar (Severinghaus et al., 2003), or a constant value (Orsi, 2013;Kobashi et al., 2015).All these methods of correction generate similar surface temperature histories (Kobashi et al., 2010(Kobashi et al., , 2015)).Other possible causes for the systematic offset are related to the standardizations to the atmosphere (in this case both nitrogen and argon isotopes can be affected) or methodological differences during the extraction from ice samples (Kobashi et al., 2008b).In these cases, a constant shift should be a better solution.
Some uncertainties remain regarding the bubble air pressures for the modeling of post-coring fractionation.First, Lipenkov (2000) reported that bubble air pressure increases toward deeper depths through the increase of ice loads, which should have induced a decrease in δAr/N 2 toward deeper depths.However, the δAr/N 2 data not exhibit any trends with depth (Fig. 7), indicating that some other processes (e.g, changes in bubble diameters, S/V , and relaxation of ice after coring especially at depths deeper than 300 m; Gow and Williamson, 1975) may have canceled the depth effect.At even deeper depths, where the exist as clathrate, the pressure between ice and clathrate boundaries can be estimated from the dissociation pressures of clathrates and should be independent of depth (Ikeda-Fukazawa et al., 2005).In future studies, it will be necessary to consider changes in each parameter in ice cores and investigate postcoring fractionation.Second, we identified that overloading pressure at the bubble close-off depth plays an important role in the post-bubble close-off fractionation in the firn.These pressure anomalies should also remain in ice cores and play some roles for the post-coring fractionation.For example, the relationship of δAr/N 2 with temperatures and accumulation rates in GISP2 may have been overestimated by the imprints of differential post-coring fractionations owing to the different bubble pressures induced by temperatures and accumulation rates at the time of the bubble close-off.Of course, the imprints of the post-coring fractionation increase if the duration of storage is longer at warmer temperatures, emphasizing the need for colder storage temperatures and the timing of measurements to recover the original signals.
For future studies on δAr/N 2 or δO 2 /N 2 in ice cores, the following suggestions should be taken into account.First, the solubility and diffusivity of argon, oxygen, and nitrogen in ice are not well constrained (Salamatin et al., 2001;Ikeda-Fukazawa et al., 2005;Bereiter et al., 2014).As precise δAr/N 2 or δO 2 /N 2 data from various ice cores are building up, the reanalyses from these cores could provide stronger constraints on the permeability.Second, although δAr/N 2 is less susceptible to the post-coring gas loss than δO 2 /N 2 , we have shown that ice core preservation is critical to retrieve the original δAr/N 2 signals.To preserve original signals, ice cores need to be stored at low temperatures (ideally < −50 • C) (Ikeda-Fukazawa et al., 2005;Bereiter et al., 2009;Landais et al., 2012) and/or to be analyzed soon after the coring.Third, we also found that the use of large ice samples (600-700 g) for each analysis reduced the noise in δO 2 /N 2 and δAr/N 2 substantially (Headly, 2008), compared to the data from smaller samples in GISP2 (Suwa and Bender, 2008b).This observation emphasizes the importance of the sample size.Fourth, observations of the bubbles in the firn and ice cores, especially of the microbubbles (e.g., numbers, volume contributions, pressure, and gas composition), are lacking, which are critical for further advances in understanding permeation in the firn and ice cores.Fifth, we have shown that δAr/N 2 could be estimated from local temperatures and accumulation rates.Therefore, combined with nitrogen and argon isotopes, it may be possible to retrieve the information of past temperatures and accumulation rates from δAr/N 2 in ice cores.Finally, the high-resolution analyses (10-20 years) provided key observations for the effects of the accumulation rates and temperatures on the permeation, which warrants further similar studies along with surface temperature reconstructions.

Conclusions
Gas fractionation after bubble close-off in the firn is complex and the associated processes are poorly understood, especially in ice cores.In this study, we investigated the gas permeation processes in the firn and ice cores using high-resolution δAr/N 2 data from the GISP2, NGRIP, and Dome Fuji ice cores for the past few millennia.We found that δAr/N 2 on the gas age in the GISP2 ice core is significantly negatively correlated with the accumulation rate and positively with air contents over the past 6000 years.Furthermore, the precise surface temperatures (Kobashi et al., 2011) and accumulation rates (Alley et al., 1997;Cuffey and Clow, 1997) over the past 4000 years from the GISP2 ice core have nearly equal control on the δAr/N 2 variations over the past 4000 years with the sensitivities of 0.72 (‰ • C −1 ) and −0.58 (‰ (0.01 m ice year −1 ) −1 ).To understand the processes of the δAr/N 2 fractionation, we applied a permeation model (Ikeda-Fukazawa et al., 2005) in which air in the bubbles leaks out by steric diffusion through ice crystals, driven by the pressure gradients between the bubbles and the atmosphere.The permeation model in the firn was applied considering two processes on the bubbles: pressure sensitive process (e.g., microbubbles) and normal bubble process.Microbubbles are believed to form near the surface.Therefore, by the time the microbubbles reach the bubble close-off depth, they develop pressures as high as overloading ice pressure that are strongly associated with changes in the accumulation rates at the surface.Several lines of evidence indicate that the pressure sensitive process occurs on a larger air fraction than that only from the microbubbles.On the other hand, the normal bubbles develop slightly higher pressures than that of the atmosphere at the bubble close-off depth such that the permeation in the firn is limited (> −0.5 ‰).The model also indicates that δAr/N 2 of the microbubbles is negatively correlated with changes in accumulation rates through increases in the overloading pres-sures, although it underestimates the magnitude observed in the GISP2 ice core.Colder temperatures are found to induce more depletions in δAr/N 2 through higher overloading pressure (thicker firn) and longer exposure time to the permeation, which explains a larger depletion in the Dome Fuji ice core.Further understanding of the gas permeation processes in the firn may lead to a new tool for estimating the past accumulation rates and/or surface temperatures.

Figure 3 .
Figure 3.Comparison of δAr/N 2 for shallow GISP2 cores (124-214 m) measured in different periods.Color data points (blue triangles, orange squares, and green diamonds) are individual data fromBender et al. (1995), and black data points with error bounds are fromKobashi et al. (2008b).We did not use shallower data ofBender et al. (1995) as they exhibit depletions similar to our shallow NGRIP data (Fig.8); an anomalous value (−16.91 ‰ at 145.4 m) in the Bender data set was also excluded.Squares, diamonds, and triangles represent the data measured after 1 week, 3 months, and 7 months of coring, respectively(Bender et al., 1995).The average difference between the Kobashi and Bender data sets is −1.51 ± 0.58 ‰, which we interpret as the post-coring fractionation for GISP2.

Figure 6 .
Figure 6.Observed and modeled δAr/N 2 over the Holocene, and decomposition of δAr/N 2 into the effects of accumulation rates and temperatures.Note that the observed δAr/N 2 are corrected values for the post-coring fractionation (1.5 ‰ added).(a) Observed and modeled δAr/N 2 .(b) Decomposition of δAr/N 2 into the effects of temperatures and accumulation rates using multiple linear regression (see text).

Figure 8 .
Figure 8. δ 15 N, δ 40 Ar/4, and δAr/N 2 in the NGRIP ice core from shallower depths (60-100 m).(a) δ 15 N, (b) δ 40 Ar/4, (c) δAr/N 2 , (d) estimated original air fractions, and (e) estimated original δAr/N 2 .The estimated original air fractions relative to the value at 75.6 m were calculated with a mass balance calculation, assuming that δ 15 N in the lock-in zone is constant with the value of 0.289 ‰ at 75.6 m and δ 15 N of the ambient air is 0.0 ‰.From the calculated original air fraction, the original δAr/N 2 were estimated again by the mass balance calculation, assuming that the ambient δAr/N 2 is 0.0 ‰.The green shaded area indicates the lock-in zone.Black dotted lines in δ 15 N, δ 40 Ar, and the estimated original air fraction are the values at 75.6 m (red dotted line).

Figure 9 .
Figure 9. Simulated δAr/N 2 vs. depth relationship in the microbubbles with a temperature of −31 • C, accumulation rate of 0.24 m ice year −1 , and microbubble contribution of 2 %.(a) Density and closed porosity (s c ).(b) δAr/N 2 and δO 2 /N 2 .(c) Air content and air pressure in the microbubbles.(d) Nitrogen and argon concentrations.

Figure 10 .
Figure 10.Traces of simulated δAr/N 2 changes for each annual layer for the normal bubbles.The model calculates bubble generation for each annual layer, gas permeation into open air, and finally trapping into ice (see text).The model is calculated assuming an equilibrium state with a temperature of −31 • C, accumulation rate of 0.24 m ice year −1 , and microbubble contribution of 2 % (same as for Fig. 9).(a) Changes in the volumes of the normal bubbles for each annual layer induced by density changes with depth.(b) Nitrogen concentrations as in (a).(c) Argon concentrations as in (a).(d) δAr/N 2 as in (a).(e) Air contents with depth, δAr/N 2 , and C(l) for the bulk of the normal bubbles (sum of the values in annual layers for each depth).Different colors (a-d) indicate values for each annual layer, showing how the bubbles generated in different annual layers evolve with time.
et al.: Post-bubble close-off fractionation of gases in polar firn and ice cores

Figure 11 .
Figure 11.The simulated δAr/N 2 fractionation with depth in the firn for the normal bubbles and microbubbles with different temperatures and accumulation rates.Microbubble contribution was set to 2 % except in (d); see also Table 3.(a) δAr/N 2 changes in the normal bubbles.(b) δAr/N 2 changes in the microbubbles.(c) Total δAr/N 2 changes (changes in the sum of the microbubbles and normal bubbles).(d) Total δAr/N 2 changes as in (c) but with different microbubble contributions (0.3 to 3 %) with a temperature of −30 • C and accumulation rate of 0.25 m ice year −1 .

Figure 12 .
Figure 12.Two model experiments for the effects of surface temperatures and accumulation rates on the overloading pressure at the bubble close-off depth.(a) Input data for the accumulation rate (0.2 m ice year −1 ) and surface temperature (−30 • C) with 20-year anomalies (+0.06 m year −1 and −5 • C) for the 1000-981 BP, respectively.When one input was used for an experiment, the other was set constant.Zero in (a) indicates the central year (model year 990 BP) of the anomalies.(b) Temperatures at the bubble close-off depth.(c) Firn thickness.(d) Overloading pressures at the bubble close-off depth.The orange line is the accumulation rate experiment, and the blue line is the temperature experiment.Numbers on peaks in (b-d) are the lag in years from the central of the initial anomalies in (a).

Table 2 .
(a) Estimated post-coring fractionation on δAr/N 2 .The values are averages over the past 2100 years for GISP2 and Dome Fuji.NGRIP shallow and deep are averages of the corresponding depths defined in the text.(b)kN 2 X N 2 and k Ar X Ar (m s −1 mol mol −1 ice MPa −1 ) at various temperatures.See also Fig.2.

Table 3 .
Modeled and observed δAr/N 2 in various conditions with microbubble contribution of 2 %.In the first column, T indicates temperature ( • C) and A indicates accumulation rate (m ice year −1 ).Duration is the time that bubbles experience from the depth of 20 % bubble closure to the depth of complete bubble close-off.Average pressure is the average overloading pressure between the depths of the 20 % bubble closure and complete bubble close-off.The average depth is the middle depth between the 20 % bubble closure and complete bubble close-off.Depth width is the depth range from 20 to 100 % bubble closed.Microbubbles, normal bubbles, and total δAr/N 2 are the values after all the bubbles are closed (i.e., in ice cores).Observed δAr/N 2 are the values corrected for the post-coring fractionation in Table2.