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

Abstract

Gases in ice cores are invaluable archives of past environmental changes (e.g., the past atmosphere). However, gas fractionation processes after bubble closure in the firn are poorly understood, although increasing evidence indicates preferential leakages of smaller molecules (e.g., neon, oxygen, and argon) from the closed bubbles through the ice matrix. These fractionation processes are believed to be responsible for the observed millennial $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ variations in ice cores, linking ice core chronologies with orbital parameters. In this study, we investigated high-resolution $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ of the GISP2 (Greenland Ice Sheet Project 2), NGRIP (North Greenland Ice Core Project), and Dome Fuji ice cores for the past few thousand years. We find that $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ at multidecadal resolution on the “gas-age scale” in the GISP2 ice core has a significant negative correlation with accumulation rate and a positive correlation with air contents over the past 6000 years, indicating that changes in overloading pressure induced $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ fractionation in the firn. Furthermore, the GISP2 temperature and accumulation rate for the last 4000 years have nearly equal effects on $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ with sensitivities of 0.72 $\pm$ 0.1 ‰ ${}^{\circ }$C${}^{-\mathrm{1}}$ and$-$0.58 $\pm$ 0.09 ‰ (0.01 m ice year${}^{-\mathrm{1}}{)}^{-\mathrm{1}}$, respectively. To understand the fractionation processes, we applied a permeation model for two different processes of bubble pressure build-up in the firn, “pressure sensitive process” (e.g., microbubbles: 0.3–3 % of air contents) with a greater sensitivity to overloading pressures and “normal bubble process”. The model indicates that $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in the bubbles under the pressure sensitive process are negatively correlated with the accumulation rate due to changes in overloading pressure. On the other hand, the normal bubbles experience only limited depletion ($<$ 0.5 ‰) in the firn. Colder temperatures in the firn induce more depletion in $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ through thicker firn. The pressure sensitive bubbles are so depleted in $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ at the bubble close-off depth that they dominate the total $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ changes in spite of their smaller air contents. The model also indicates that $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ of ice cores should have experienced several per mil of depletion during the storage 14–18 years after coring. Further understanding of the $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ fractionation processes in the firn, combined with nitrogen and argon isotope data, may lead to a new proxy for the past temperature and accumulation rate.

Introduction

Atmospheric gases trapped in the firn layer (unconsolidated snow layer; $\sim$ 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 ($\mathrm{\Delta }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 $\sim$ 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 bubble-close-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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ or $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ variations in the time domain (i.e., ice cores) are still poorly understood.

$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ vs. accumulation rates or air contents in GISP2 over the past 6000 years. Note that $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ data. Two $\mathit{\sigma }$ error bounds are shown, which were estimated by 1000 iterations of Monte Carlo simulations. Accumulation rates (m ice year${}^{-\mathrm{1}})$ (black line) were filtered by 21-year RMs. Note that the $y$ axis for the accumulation rate is reversed. $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ vs. accumulation rates are significantly negatively correlated over the past 6000 years ($r=-0.35$, $p=0.03$). $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variation can explain 7–14 % of the total variance of the air contents.

[Figure omitted. See PDF]

On a longer timescale (i.e, orbital), variations in $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ fractionation (insolation hypothesis) (Bender, 2002; Fujita et al., 2009). In addition, air content in ice cores is also found to co-vary with $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ on the orbital timescale, indicating common causes (Raynaud et al., 2007; Lipenkov et al., 2011). According to this hypothesis, the orbital signals in $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variability in three ice cores (GISP2, NGRIP, and Dome Fuji) as well as gas loss processes during storage. $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations are generally highly correlated with $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ in ice cores, suggesting similar processes driving these variations (Bender et al., 1995). As $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ that varied in the atmosphere by $\sim$ 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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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.

Data description

$\mathit{\delta }$Ar$/$N${}_{\mathrm{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, 2015), and Dome Fuji data are new. Importantly, storage histories of these cores (i.e., temperatures) are known and methods for measuring $\mathit{\delta }$Ar$/$N${}_{\mathrm{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, 2010) and around the 8.2 ka event (8100 $\pm$ 500 years before present (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 ${}^{\circ }$C (Table 1). However, the accumulation rate of NGRIP ($\sim$ 0.19 m ice year${}^{-\mathrm{1}}$) is 20 % less than that of GISP2 (0.24 m ice year${}^{-\mathrm{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 ${}^{\circ }$C and a mean accumulation rate of $\sim$ 0.03 m ice year${}^{-\mathrm{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 $\sim$ 10 % of the estimated gas age–ice age difference (Goujon et al., 2003). To investigate the $\mathit{\delta }$Ar$/$N${}_{\mathrm{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., $\mathit{\delta }$Ar$/$N${}_{\mathrm{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.

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
				(${}^{\circ }$C)	(m ice year${}^{-\mathrm{1}}$)	(m ice year${}^{-\mathrm{1}}$)
GISP2	72.59${}^{\circ }$ N	38.46${}^{\circ }$ W	3203	$-$31.0	0.24	0.013
NGRIP	75.1${}^{\circ }$ N	42.32${}^{\circ }$ W	3230	$-$31.5	0.19	0.008
Dome Fuji	77.32${}^{\circ }$ S	39.67${}^{\circ }$ E	3810	$-$54.3	0.03	NA

GISP2 and NGRIP ice cores were analyzed for $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ $\sim$ 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 ${}^{\circ }$C in a commercial freezer until they were moved to a freezer ($-$36 ${}^{\circ }$C) at the National Ice Core Laboratory (NICL) in February 1993 (G. Hargreaves, personal communication, 2015). The ice samples were then cut and moved to the Scripps Institution of Oceanography, where $\mathit{\delta }$Ar$/$N${}_{\mathrm{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.2 m) were drilled in summer 1999 (Dahl-Jensen et al., 2002). Shallower parts (64.6–254.4 m) were stored in a freezer at the University of Copenhagen at around $-$24 ${}^{\circ }$C (J. P. Steffensen, personal communication, 2015), and deeper parts (255.5–445.2 m) were in a freezer of a commercial facility rented by the Alfred Wegener Institute (AWI) at $-$30 ${}^{\circ }$C (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 ${}^{\circ }$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 ${}^{\circ }$C, with a short period (2.5 months) at $<$ $-$25 ${}^{\circ }$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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ as follows: $$\mathit{\delta }\text{Ar}/{\text{N}}_{\mathrm{2}}=\left[\right(\text{Ar}/{\text{N}}_{\mathrm{2}}{)}_{\text{sample}}/(\text{Ar}/{\text{N}}_{\mathrm{2}}{)}_{\text{standard}}-\mathrm{1}]{10}^{\mathrm{3}}\left(\mathrm{‰}\right),$$ 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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$. All $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ data presented in this study were corrected for gravitational and thermal fractionations in the firn using the conventional method based on ${\mathit{\delta }}^{15}$N (Bender et al., 1995; Severinghaus and Battle, 2006; Severinghaus et al., 2009) as follows: $$\mathit{\delta }\text{Ar}/{\text{N}}_{\mathrm{2}\text{corr}}=\mathit{\delta }\text{Ar}/{\text{N}}_{\mathrm{2}}-11{\mathit{\delta }}^{15}\text{N}.$$ The coefficient 11 arises because the mass difference of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (${}^{40}$Ar and ${}^{29}$N${}_{\mathrm{2}})$ is 11 times larger than that of the nitrogen isotopes (${}^{29}$N${}_{\mathrm{2}}$ and ${}^{28}$N${}_{\mathrm{2}})$. This coefficient is replaced with 12 for the calculation of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}\text{corr}}$ for NGRIP and Dome Fuji because the mass difference between ${}^{40}$Ar and ${}^{28}$N${}_{\mathrm{2}}$ is 12. As the temperature sensitivities of ${\mathit{\delta }}^{15}$N and $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}\text{corr}}$ variations in the ice cores can be attributed only to gas loss. $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}\text{corr}}$ of the GISP2 data using the mass 28 or 29 leads to negligible differences (an average difference is 0.4 $\times$ 10${}^{-\mathrm{3}}$ ‰ and the standard deviation is 0.94 $\times$ 10${}^{-\mathrm{3}}$ ‰) which are much smaller than the measurement uncertainty of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$(1$\mathit{\sigma }$ $<$ 0.7 ‰). We also note that the standard deviation (0.07 ‰) of ${\mathit{\delta }}^{15}$N $\times$ 11 in GISP2 is much smaller than standard deviation of raw $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (1.33 ‰) over the past 6000 years, indicating that the variations of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}\text{corr}}$ mostly originate from the raw $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ and not from ${\mathit{\delta }}^{15}$N. For the sake of simplicity, we denote all the $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}\text{corr}}$ as $\mathit{\delta }$Ar$/$N${}_{\mathrm{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$\mathit{\sigma }$.

Post-coring fractionation

Before evaluating $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in ice cores for the changes that have occurred in the firn, it is necessary to consider the post-coring 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 $\sim$ 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${}_{\text{ice}}^{-\mathrm{1}}$) 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): $$U_m=U_m^{\mathrm{0}}-k_{\text{m}}{X}_m\left(P^i,Z_m^i,-,P^a,Z_m^s\right)S/Vt,$$ where $U_m^{\mathrm{0}}$ (mol mol${}_{\mathrm{ice}}^{-\mathrm{1}}$) is the original concentration of species $m$. $k_{\text{m}}$ (m s${}^{-\mathrm{1}})$ is the mass transfer coefficient, and it is equal to $D_m/\mathrm{\Delta }l$, where $D_m$ (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}})$ is the diffusion coefficient of the species $m$ and $\mathrm{\Delta }l$ (m) is the thickness of the surface layer of ice (Ikeda-Fukazawa et al., 2005). $X_m$ (mol mol${}_{\mathrm{ice}}^{-\mathrm{1}}$ M Pa${}^{-\mathrm{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_m^i$ and $Z_m^s$ are molar fractions of species $m$ in the bubbles and in the atmosphere, respectively. $S$ (m${}^{\mathrm{2}})$ and $V$ (m${}^{\mathrm{3}})$ represent the surface area and the volume of an ice sample such that $S/V$ can be understood as specific surface area (m${}^{-\mathrm{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 $\times$ 10${}^{-\mathrm{5}}$ mol in 1 mol of ice (a typical air content in ice cores). $U_m^{\mathrm{0}}$ for each gas is calculated from the total air content multiplied by the atmospheric molar ratio of each gas. In this case, $Z_m^i$ and $Z_m^s$ 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${}^{-\mathrm{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 $\sim$ 5 mm before the analyses (Kobashi et al., 2008b, 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_{\text{Ar}})$ and solubility ($X_{\text{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_{\text{Ar}}{X}_{\text{Ar}}$ (i.e., $D_{\text{Ar}}/\mathrm{\Delta }l\cdot X_{\text{Ar}})$ in relation to those for nitrogen and oxygen (Fig. 2). $K_{{\text{N}}_{\mathrm{2}}}{X}_{{\text{N}}_{\mathrm{2}}}$ and $k_{{\text{O}}_{\mathrm{2}}}{X}_{{\text{O}}_{\mathrm{2}}}$ in different temperatures can be estimated using Eqs. (4) and (2) with $\mathrm{\Delta }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).

$k_m\cdot X_m$ for oxygen, argon, and nitrogen for different temperatures. Ar(I) and Ar(II) were calculated from Eqs. (4) and (5) and Eqs. (6), respectively (see text).

[Figure omitted. See PDF]

The first estimate of Ar(I) uses a diffusion coefficient ($D_{\text{Ar}}$; 4.0 $\times$ 10${}^{-11}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}})$ of argon at 270 K calculated from molecular dynamic simulations with those of nitrogen ($D_{{\text{N}}_{\mathrm{2}}}$; 2.1 $\times$ 10${}^{-11}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}})$ and oxygen ($D_{{\text{O}}_{\mathrm{2}}}$; 4.7 $\times$ 10${}^{-11}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}})$ (Ikeda-Fukazawa et al., 2004). Owing to the molecular–size-dependent 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_{\text{Ar}}$ and $X_{\text{Ar}}$ were estimated assuming that the geometrical relationship between $D_{{\text{N}}_{\mathrm{2}}},D_{\text{Ar}},$ and $D_{{\text{O}}_{\mathrm{2}}}$ at 270 K from the molecular dynamics simulations holds for $k_{\text{Ar}}$ and $X_{\text{Ar}}$ at different temperatures as follows: $$\begin{array}{cc}& & k_{\text{Ar}}=k_{{\text{O}}_{\mathrm{2}}}-(D_{{\text{O}}_{\mathrm{2}}\text{at}\mathrm{\hspace{0.17em}270}\text{K}}-D_{\text{Ar}\text{at}\mathrm{\hspace{0.17em}270}\text{K}})/& (D_{{\text{O}}_{\mathrm{2}}\text{at}\mathrm{\hspace{0.17em}270}\text{K}}-D_{{\text{N}}_{\mathrm{2}}\text{at}\mathrm{\hspace{0.17em}270}\text{K}})(k_{{\text{O}}_{\mathrm{2}}}-k_{{\text{N}}_{\mathrm{2}}}),\\ & & X_{\text{Ar}}=X_{{\text{O}}_{\mathrm{2}}}-(D_{{\text{O}}_{\mathrm{2}}\text{at}\mathrm{\hspace{0.17em}270}\text{K}}-D_{\text{Ar}\text{at}\mathrm{\hspace{0.17em}270}\text{K}})/& (D_{{\text{O}}_{\mathrm{2}}\text{at}\mathrm{\hspace{0.17em}270}\text{K}}-D_{{\text{N}}_{\mathrm{2}}\text{at}\mathrm{\hspace{0.17em}270}\text{K}})(X_{{\text{O}}_{\mathrm{2}}}-X_{{\text{N}}_{\mathrm{2}}}).\end{array}$$ Second, we estimated Ar(II) from an observation that $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in ice is often depleted about half of $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ in ice cores (Bender et al., 1995). To satisfy this condition, $k_{\text{Ar}}{X}_{\text{Ar}}$ can be written as $$k_{\text{Ar}}{X}_{\text{Ar}}=(k_{{\text{N}}_{\mathrm{2}}}{X}_{{\text{N}}_{\mathrm{2}}}+k_{{\text{O}}_{\mathrm{2}}}{X}_{{\text{O}}_{\mathrm{2}}})/2.$$ Estimated $k_{\text{Ar}}{X}_{\text{Ar}}$ for Ar(I) and Ar(II) are higher than $k_{{\text{N}}_{\mathrm{2}}}{X}_{\text{N}\mathrm{2}}$ and increase with temperatures, resulting in a general depletion of $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ than that of Ar(II) owing to faster permeation of argon. With the two estimates of $k_{\text{Ar}}{X}_{\text{Ar}}$, we explore the range of uncertainties associated with argon permeation.

(a) Estimated post-coring fractionation on $\mathit{\delta }$Ar$/$N${}_{\mathrm{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) $k_{{\mathrm{N}}_{\mathrm{2}}}{X}_{{\mathrm{N}}_{\mathrm{2}}}$ and $k_{\text{Ar}}{X}_{\text{Ar}}$ (m s${}^{-\mathrm{1}}$ mol mol${}_{\text{ice}}^{-\mathrm{1}}$ MPa${}^{-\mathrm{1}})$ at various temperatures. See also Fig. 2.

(a)	Period 1			Period 2			$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (‰)
	Duration	Temp.		Duration	Temp.		Est.	Observation	Est. average
	(years)	(${}^{\circ }$C)		(years)	(${}^{\circ }$C)		post-coring	in	values
							depletion	ice cores	before coring
GISP2	2	$-$29		12	$-$36		1.5 $\pm$ 0.6	$-$3.9 $\pm$ 0.2	$-$2.4 $\pm$ 0.6
NGRIP shallow	12	$-$24		2	$-$30		3.0 $\pm$ 1.2	$-$6.3 $\pm$ 0.1	$-$3.3 $\pm$ 1.2
NGRIP deep	14	$-$30		–	–		2.5 $\pm$ 1.0	$-$5.9 $\pm$ 0.2	$-$3.4 $\pm$ 1.0
Dome Fuji	0.2	$-$25		18	$-$50		1.5 $\pm$ 0.7	$-$7.8 $\pm$ 0.3	$-$6.3 $\pm$ 0.8
(b)	$-$24 ${}^{\circ }$C	$-$25 ${}^{\circ }$C	$-$29 ${}^{\circ }$C	$-$30 ${}^{\circ }$C	$-$36 ${}^{\circ }$C	$-$50 ${}^{\circ }$C
$k_{{\mathrm{N}}_{\mathrm{2}}}{X}_{{\mathrm{N}}_{\mathrm{2}}}$	7.54 $\times$ 10${}^{-18}$	7.33 $\times$ 10${}^{-18}$	6.54 $\times$ 10${}^{-18}$	6.36 $\times$ 10${}^{-18}$	5.31 $\times$ 10${}^{-18}$	3.37 $\times$ 10${}^{-18}$
$k_{Ar}{X}_{Ar}$ (Ar I)	2.78 $\times$ 10${}^{-17}$	2.69 $\times$ 10${}^{-17}$	2.34 $\times$ 10${}^{-17}$	2.26 $\times$ 10${}^{-17}$	1.82 $\times$ 10${}^{-17}$	1.05 $\times$ 10$-$${}^{17}$
$k_{\text{Ar}}{X}_{\text{Ar}}$ (Ar II)	2.27 $\times$ 10${}^{-17}$	2.19 $\times$ 10${}^{-17}$	1.91 $\times$ 10${}^{-17}$	1.84 $\times$ 10${}^{-17}$	1.49 $\times$ 10${}^{-17}$	8.71 $\times$ 10${}^{-18}$

In a pioneering study by Bender et al. (1995), $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in a shallow core of GISP2 was analyzed after 1 week, 3 months, and 7 months of drilling in 1989 to study the time-dependent 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${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in GISP2 to be $-$1.5 $\pm$ 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 $\pm$ 0.2 MPa and 0.8 $\pm$ 0.3 MP for two different estimates of $k_{\text{Ar}}{X}_{\text{Ar}}$ of Ar(I) and Ar(II), respectively, which agree with the normal bubble pressure at 150–200 m depths in Vostok (Lipenkov, 2000). Using the estimated bubble air pressure and aforementioned parameters, the amounts of depletion in $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ after coring are estimated as $-$3.0 $\pm$ 1.2, $-$2.5 $\pm$ 1.0, and 1.5 $\pm$ 0.7 ‰ for NGRIP shallow, NGRIP deep, and Dome Fuji, respectively (Table 2). As a result, it is possible to derive the original $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ values before coring for GISP2, NGRIP shallow, and NGRIP deep, and Dome Fuji as $-$2.4 $\pm$ 0.6, $-$3.3 $\pm$ 1.2, $-$3.4 $\pm$ 1.0, and $-$6.3 $\pm$ 0.8 ‰, respectively (Table 2).

Comparison of $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 from Bender et al. (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 $\pm$ 0.58 ‰, which we interpret as the post-coring fractionation for GISP2.

[Figure omitted. See PDF]

${\mathit{\delta }}^{15}$N, ${\mathit{\delta }}^{40}$Ar/4, and $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ from the GISP2 ice core over the Holocene (Kobashi et al., 2008b). The grey arrow indicates the brittle zone (Gow et al., 1997).

[Figure omitted. See PDF]

Post-bubble close-off fractionation in firn: empirical evidence

GISP2 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variation over the Holocene

The $\mathit{\delta }$Ar$/$N${}_{\mathrm{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, $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in these depths are highly variable (Fig. 4). It is noted that ${\mathit{\delta }}^{15}$N and ${\mathit{\delta }}^{40}$Ar exhibit little influences from the anomalous $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ fractionation, indicating that the processes are mass independent 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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations. We found a significant negative correlation ($r=-0.35$, $p=0.03$) between $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ on the gas-age scale and the accumulation rate for the past 6000 years, a time interval in which the abnormal $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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.

Observed and modeled $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$, surface temperatures, and accumulation rates from the GISP2 ice core over the past 4000 years. Note that the observed $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ is corrected for the post-coring fractionation (1.5 ‰ added). (a) $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ and surface temperatures (Kobashi et al., 2011). Ages of the temperatures were adjusted for the lag (68 years). (b) $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ from multiple linear regression (see text). (d) Observed and modeled $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ of multiple linear regression using ${\mathit{\delta }}^{18}$O${}_{\text{ice}}$ as a temperature proxy (see text).

[Figure omitted. See PDF]

The subset of the GISP2 data covering the past 4000 years provides a unique opportunity to investigate $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$. It is found that the surface temperature is positively correlated with $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ with a multidecadal lag. On the other hand, the accumulation rate is negatively correlated with $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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., $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ over the past 4000 years. As expected, the model of the multiple linear regression captures the $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ indicates that changes in the surface temperature and accumulation rate play important roles in controlling the $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations. From the multiple linear regression, $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ on the gas ages in GISP2 can be expressed by temperature (${}^{\circ }$C) and accumulation rate (m ice year${}^{-\mathrm{1}}$) as a function of time after adjusting for the lags: $$\begin{array}{cc}& \mathit{\delta }\text{Ar}/{\text{N}}_{\mathrm{2}}\left(t\right)& =A\times \text{temperature}(t+t_{\text{temp}})& +B\times \text{accumulation}(t+t_{\text{accm}})+C,\end{array}$$ where $A=0.72\pm 0.06$ ‰ ${}^{\circ }$C${}^{-\mathrm{1}}$, $B=-58.8\pm 4.3$ ‰ (m year${}^{-\mathrm{1}}{)}^{-\mathrm{1}}$, $C=32.7\pm 1.8$ ‰, and $t,t_{\text{temp}}$, and $t_{\text{accm}}$ are time, lag for temperature, and lag for accumulation rate, respectively (all in years).

Next, we attempted to use oxygen isotopes of ice (${\mathit{\delta }}^{18}$O${}_{\text{ice}})$ as a temperature proxy for the same regression analyses of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ since we do not have the N${}_{\mathrm{2}}$/Ar isotope-based temperature information before 4000 years BP. Temperature records derived from ${\mathit{\delta }}^{18}$O${}_{\text{ice}}$ can be quite noisy but stacking several ${\mathit{\delta }}^{18}$O${}_{\text{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ performs not as well as the one with the temperature and accumulation rate records for the last 4000 years based on Ar/N${}_{\mathrm{2}}$ isotope values (Fig. 5c) but does slightly better than the correlations with the temperature or accumulation rate individually (Fig. 5a, b).

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

[Figure omitted. See PDF]

The $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ regression model with the ${\mathit{\delta }}^{18}$O${}_{\text{ice}}$-based temperatures and accumulation rates can span the entire Holocene, including the periods when the observed $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ exhibit rather constant values of $-$1 to $-$3 ‰ during the Holocene (Fig. 6). Interestingly, the model indicates that the constant $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations remained higher or noisier from the early Holocene to $\sim$ 6000 BP than for the later period, which probably make it difficult to decipher the original multidecadal to centennial signals in $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (Fig. 6).

Surface temperatures, accumulation rates, ${\mathit{\delta }}^{15}$N, $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ for GISP2, NGRIP, and Dome Fuji over the past 2100 years. (a) Surface temperatures for GISP2 (black) and NGRIP (blue) (Kobashi et al., 2015). (b) Accumulation rates in 21-year RMs for GISP2 (black: Alley et al., 1997; Cuffey and Clow, 1997) and NGRIP (blue: Gkinis et al., 2014). (c) Raw ${\mathit{\delta }}^{15}$N and spline for NGRIP and GISP2 (Kobashi et al., 2010, 2015). (d–f) $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ and the values corrected for the post-coring fractionation for GISP2, NGRIP, and Dome Fuji. Blue and black lines are the raw and corrected values for the post-coring fractionations, respectively. A red point with error bounds (2$\mathit{\sigma })$ indicates estimated $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ for Dome Fuji using Eq. (7).

[Figure omitted. See PDF]

NGRIP and Dome Fuji $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variation over the past 2100 years

The $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ data of NGRIP ice cores allows for a good comparison with the GISP2 data (Fig. 7). The averages of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ for the past 2100 years are $-$3.36 and $-$2.40 ‰ for NGRIP and GISP2, respectively (Fig. 7). The $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variability in NGRIP (1$\mathit{\sigma }=0.91$ ‰) over the past 2100 years is 24 % smaller than that of GISP2 (1$\mathit{\sigma }=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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ than the one used for GISP2 (Kobashi et al., 2008b, 2015). $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 errors $=$ 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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ and the temperature or accumulation rate of GISP2 in Eq. (7), we can calculate expected $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ from Eq. (7) were calculated as 0.3 $\pm$ 1.3 and $-$6.4 $\pm$ 1.2 ‰, respectively. The value for NGRIP is significantly higher than the observed value of $-$3.3 $\pm$ 1.2 ‰ corrected for the post-coring fractionation (Table 2). For Dome Fuji, the value is similar to the observed $-$6.3 $\pm$ 0.8 ‰ corrected for the post-coring fractionation (Fig. 7, Table 2). This may indicate that the relationship of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ with the temperature and accumulation rate becomes non-linear when the firn thickness becomes thinner than that of GISP2 as $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ is not expected to be positive without the existence of clathrate (see later discussion).

$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ ice core data of NGRIP from the depth range 64.6–80 m exhibit some interesting features (Fig. 8). The depth from $\sim$ 60 to 78 m corresponds to the lock-in zone in NGRIP, where vertical mixing of gas is limited so that ${\mathit{\delta }}^{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 ${\mathit{\delta }}^{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 ${\mathit{\delta }}^{15}$N and ${\mathit{\delta }}^{40}$Ar toward shallower depths in the lock-in-zone as the result of mixing with ambient air (Fig. 8d). Based on isotope mass balance, we calculated the original $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ values, which exhibited highly depleted values as low as $-$50 ‰ (Fig. 8e). The depleted $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in the lock-in zone provides a clue to the processes of gas loss in the firn (see later discussion).

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

[Figure omitted. See PDF]

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 $\mathrm{\mu }$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 bubbles, 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, $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ or $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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). $$\begin{array}{cc}& & U_m(l+\mathrm{1})=& U_m\left(l\right)-k_m{X}_m\left(P^i,\left(l\right),Z_m^i,\left(l\right),-,P^a,Z_m^s\right)\\ & \cdot \left({\displaystyle \frac{S}{V}},\left(,l,\right)\right)s_{\text{o}}/s\left(l\right)tC\left(l\right)\end{array}$$

$$Z_m^i\left(l\right)={\displaystyle \frac{U_m\left(l\right)}{U_{\text{Ar}}\left(l\right)+U_{{\mathrm{O}}_{\mathrm{2}}}\left(l\right)+U_{{\mathrm{N}}_{\mathrm{2}}}}}\left(l\right),$$ where $l$ is an annual layer from the surface to below the firn layer (e.g., $l=\mathrm{1}$ to 2000), and $s_{\text{o}}/s\left(l\right)$ in a layer $l$ is the open porosity ratio. $s,$ $s_{\text{c}}$, and $s_{\text{o}}$ are the total, closed, and open porosities ($s_{\text{o}}=s-s_{\text{c}})$, respectively (Spahni et al., 2003; see also Sect. 5.2). In a steady state, $l$ can be considered as a time variable. At $l=\mathrm{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\left(l\right)$ 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\left(l\right)$ 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${}^{-\mathrm{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\left(l\right))$ in a layer $l$ from the linear relationship between the specific surface areas (m${}^{-\mathrm{1}})$ and densities $\mathit{\rho }$ from the Greenland summit (Lomonaco et al., 2011) with an equation: $S/V\left(l\right)$ (m${}^{-\mathrm{1}}$) $=$ $-$16 799 $\mathit{\rho }\left(l\right)$ (g cm${}^{-\mathrm{3}}$) $+$ 14 957. The initial gas content in the microbubbles was set at 0.3–3 % of the air content (6.53226 $\times$ 10${}^{-\mathrm{5}}$ mol $\times$ 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_{\text{o}}/s\left(l\right)$ (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.

Simulated $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ vs. depth relationship in the microbubbles with a temperature of $-$31 ${}^{\circ }$C, accumulation rate of 0.24 m ice year${}^{-\mathrm{1}}$, and microbubble contribution of 2 %. (a) Density and closed porosity ($s_{\text{c}})$. (b) $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ and $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$. (c) Air content and air pressure in the microbubbles. (d) Nitrogen and argon concentrations.

[Figure omitted. See PDF]

Figure 9 shows model results with a temperature of $-$31 ${}^{\circ }$C, an accumulation rate of 0.25 m ice year${}^{-\mathrm{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, $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ depletion is larger than $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (Fig. 9b). At the temperature of $-$30 ${}^{\circ }$C and accumulation rate of 0.25 m ice year${}^{-\mathrm{1}}$, the depletion reaches up to 133 ‰ for $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$, and 243 ‰ for $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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_{\text{c}})$ with depth as the following equation: $$\begin{array}{cc}& V_{\mathrm{0}}(l+\mathrm{1})& =a\left\{s_{\text{c}},(,l,+,\mathrm{1},),-,s_{\text{c}},\left(,l,\right),(,{\mathit{\rho }}_{\text{ice}},-,\mathit{\rho },(,l,+,\mathrm{1},\left),\right),/\right.& \left..(,{\mathit{\rho }}_{\text{ice}},-,\mathit{\rho },(,l,\left),\right)\right\}/\mathit{\rho }(l+\mathrm{1}),\end{array}$$ where $V_{\mathrm{0}}\left(l\right)$ is newly trapped air in an annual layer $l$, ${\mathit{\rho }}_{\text{ice}}$ is the density of ice, and $\mathit{\rho }\left(l\right)$ is the density at depth $l$. $s_{\text{c}}\left(l\right)$ is the closed porosity in an annual layer $l$, and $a$ is a scaling coefficient. $s_{\text{c}}$ can be written as (Schwander, 1989; Spahni et al., 2003) $$s_{\text{c}}=\left\{\begin{array}{ll}s\cdot \exp \left[75,\cdot ,\left({\displaystyle \frac{\mathit{\rho }\left(l\right)}{{\mathit{\rho }}_{\text{co}}}},-,\mathrm{1}\right)\right],& \mathrm{0}<\mathit{\rho }\left(l\right)<{\mathit{\rho }}_{\text{co}},\\ s,& \mathit{\rho }\left(l\right)>{\mathit{\rho }}_{\text{co}},\end{array}\right.$$ where ${\mathit{\rho }}_{\text{co}}$ is the density at the depth in which the air is totally enclosed in bubbles. The sum of all the newly generated air ($\sum _{l=\mathrm{1}}^{2000}{V}_{\mathrm{0}}\left(l\right))$ is set to have the air content of $6.53\times {10}^{-\mathrm{5}}$ mol per mole of ice. Then, $V_{\mathrm{0}}\left(l\right)$ was scaled accordingly using the coefficient $a$ and converted to the volume (m${}^{\mathrm{3}})$ with the atmospheric pressure (0.065 MPa) as in Fig. 10a.

Traces of simulated $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 ${}^{\circ }$C, accumulation rate of 0.24 m ice year${}^{-\mathrm{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) $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ as in (a). (e) Air contents with depth, $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$, and $C\left(l\right)$ 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.

[Figure omitted. See PDF]

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 ${}^{\circ }$C and 0.24 m ice year${}^{-\mathrm{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\left(l\right)$ in Eq. (8) was calculated from $V_{\mathrm{0}}\left(l\right)$ for each annual layer $l$ by setting the sum of $C\left(l\right)$ as $\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in each annual layer, and the air content and bulk $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ with depth at a temperature of $-$31 ${}^{\circ }$C and an accumulation rate of 0.24 m ice year${}^{-\mathrm{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 ${}^{\circ }$C (Ar(II), Fig. 2), $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ ($\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}})$ within the bubbles decreases when the permeation proceeds. At the temperature of $-$31 ${}^{\circ }$C and accumulation rate of 0.24 m ice year${}^{-\mathrm{1}}$, the $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ is limited (Fig. 10e). The depth vs. $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ relationship of the total air from the normal bubbles (Fig. 10e) indicates that the total $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ reaches the minimum of $-$0.39 ‰ at the middle of the bubble close-off depth of 73.2 m. Then, the total $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ increases to $-$0.29 ‰ as a large amount of ambient air with $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ $=$ 0 is trapped in these depths (Fig. 10a, d, e).

Modeled and observed $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in various conditions with microbubble contribution of 2 %. In the first column, $T$ indicates temperature (${}^{\circ }$C) and $A$ indicates accumulation rate (m ice year${}^{-\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ are the values after all the bubbles are closed (i.e., in ice cores). Observed $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ are the values corrected for the post-coring fractionation in Table 2.

	Surface condition			Model							Observed
	Temp.	Accumulation		Duration	Avg.	Avg.	Depth	Micro-	Normal	Total	$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$
	(${}^{\circ }$C)	(m ice year${}^{-\mathrm{1}}$)		(year)	pressure	depth	width	bubble	bubble	$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$	(‰)
					(MPa)	(m)	(m)	$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$	$\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$	(‰)
								(‰)	(‰)
GISP2	$-$31	0.24		31	0.54	72.8	8.7	$-$133.2	$-$0.31	$-$2.74	$-$2.4 $\pm$ 0.6
NGRIP	$-$31.5	0.19		36	0.50	68.0	8.3	$-$120.8	$-$0.34	$-$2.56	$-$3.3 $\pm$ 1.2
Dome Fuji	$-$54.3	0.03		274	0.62	89.6	10.3	$-$599.8	$-$0.43	$-$6.43	$-$6.3 $\pm$ 0.8
$T$: $-$25, $A$: 0.2	$-$25	0.2		27	0.42	54.4	6.5	$-$41.8	$-$0.33	$-$1.13	–
$T$: $-$25, $A$: 0.25	$-$25	0.25		23	0.45	59.1	6.9	$-$60.6	$-$0.30	$-$1.45	–
$T$: $-$25, $A$: 0.3	$-$25	0.3		21	0.48	63.4	7.5	$-$74.4	$-$0.28	$-$1.68	–
$T$: $-$30, $A$: 0.2	$-$30	0.2		33	0.49	65.4	8.0	$-$103.3	$-$0.33	$-$2.25	–
$T$: $-$30, $A$: 0.25	$-$30	0.25		29	0.53	71.3	8.8	$-$121.7	$-$0.30	$-$2.54	–
$T$: $-$30, $A$: 0.3	$-$30	0.3		27	0.56	76.6	9.8	$-$135.4	$-$0.27	$-$2.75	–
$T$: $-$35, $A$: 0.2	$-$35	0.2		40	0.58	79.5	9.8	$-$184.2	$-$0.33	$-$3.58	–
$T$: $-$35, $A$:0.25	$-$35	0.25		36	0.63	86.9	10.9	$-$203.0	$-$0.30	$-$3.83	–
$T$: $-$35, $A$: 0.3	$-$35	0.3		34	0.67	93.5	12.4	$-$216.6	$-$0.27	$-$4.01	–

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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (Fig. 11). Results show that the normal bubbles experience only limited $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ depletion ($>$ $-$0.5 ‰) by the different temperatures or accumulation rates we considered (Table 3). On the other hand, $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ generally reflects the variation of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in the microbubbles ($r=0.95$, Table 3). Overall, the total $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ have a higher correlation with temperatures ($r=0.97$) than with accumulation rates ($r=0.57$) in the model (Table 3).

The simulated $\mathit{\delta }$Ar$/$N${}_{\mathrm{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) $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ changes in the normal bubbles. (b) $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ changes in the microbubbles. (c) Total $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ changes (changes in the sum of the microbubbles and normal bubbles). (d) Total $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ changes as in (c) but with different microbubble contributions (0.3 to 3 %) with a temperature of $-$30 ${}^{\circ }$C and accumulation rate of 0.25 m ice year${}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

The modeled $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ agrees with the observed $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in the microbubbles and in the total air (Table 3). The variations of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in normal bubbles are limited and, clearly, microbubbles (or the pressure sensitive process) play a critical role for the variation of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ in ice cores. The $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ minima in the firn ranges from $-$14 to $-$83 ‰ depending on the temperatures and accumulation rates. The most depleted $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ with a temperature of $-$30 ${}^{\circ }$C and accumulation rate of 0.2 m ice year${}^{-\mathrm{1}}$ in Fig. 11c capture the highly depleted observation-based estimates of $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (Fig. 8e) is an evidence for the involvement of the microbubble process (or pressure sensitive process). The total $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations on the temperature (0.24 ‰ ${}^{\circ }$C${}^{-\mathrm{1}}$ for an accumulation rate of 0.25 m ice year${}^{-\mathrm{1}})$ and accumulation rate ($-$0.05 ‰ (0.01 m ice year${}^{-\mathrm{1}}{)}^{-\mathrm{1}}$ at $-$30 ${}^{\circ }$C) with the 2 % microbubble contribution (Table 3) are lower than those of the observed ones in GISP2 ice cores (0.72 $\pm$ 0.1 ‰ ${}^{\circ }$C${}^{-\mathrm{1}}$ and $-$0.58 $\pm$ 0.09 ‰ (0.01 m ice year${}^{-\mathrm{1}}{)}^{-\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (Fig. 11d). Also, the dependencies of $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ on temperatures and accumulation rates linearly increase to 0.38 ‰ ${}^{\circ }$C${}^{-\mathrm{1}}$ with an accumulation rate of 0.25 m ice year${}^{-\mathrm{1}}$, and $-$0.11 ‰ (0.01 m ice year${}^{-\mathrm{1}}{)}^{-\mathrm{1}}$ with a temperature at $-$30 ${}^{\circ }$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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 ${}^{\circ }$C) and accumulation rate (0.2 m ice year${}^{-\mathrm{1}})$ over thousands of years to reach an equilibrium state. Then, surface temperature and accumulation rate anomalies of $-$35 ${}^{\circ }$C and 0.26 m ice year${}^{-\mathrm{1}}$ for 20 years were introduced, separately (Fig. 12a). The surface anomalies of the temperature and accumulation rate were set to induce similar $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ changes of 3.5 ‰ from the relationship obtained by the multiple linear regressions on the $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ of GISP2.

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${}^{-\mathrm{1}})$ and surface temperature ($-$30 ${}^{\circ }$C) with 20-year anomalies ($+$0.06 m year${}^{-\mathrm{1}}$ and $-$5 ${}^{\circ }$C) for the years 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 year of the initial anomalies in (a).

[Figure omitted. See PDF]

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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations through the permeation.

Discussions

The processes responsible for the $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations should also play similar roles on the variations of $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 earlier studies, causes of the $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ variations in the orbital scale are also a result of the surface temperature and accumulation rate changes. We note that $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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., 2003, 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, ${\mathit{\delta }}^{18}$O and ${\mathit{\delta }}^{40}$Ar) with highly depleted $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ or $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 ${\mathit{\delta }}^{15}$N and ${\mathit{\delta }}^{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 ${\mathit{\delta }}^{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, 2011, 2015). The systematically higher ${\mathit{\delta }}^{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 ${\mathit{\delta }}^{18}$O of O${}_{\mathrm{2}}$ (Battle et al., 2011). If the enrichment of ${\mathit{\delta }}^{40}$Ar occurs in the firn, it should be correlated with $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$. Therefore, the corrections for the ${\mathit{\delta }}^{40}$Ar enrichment have been applied using $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ (Kobashi et al., 2010, 2011, 2015), $\mathit{\delta }$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, 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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ toward deeper depths. However, the $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ data do 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 bubbles 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 post-coring 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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ or $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ or $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ data from various ice cores are building up, the reanalyses from these cores could provide stronger constraints on the permeability. Second, although $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ is less susceptible to the post-coring gas loss than $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$, we have shown that ice core preservation is critical to retrieve the original $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ signals. To preserve original signals, ice cores need to be stored at low temperatures (ideally $<$ $-$50 ${}^{\circ }$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 $\mathit{\delta }$O${}_{\mathrm{2}}$$/$N${}_{\mathrm{2}}$ and $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ data from the GISP2, NGRIP, and Dome Fuji ice cores for the past few millennia. We found that $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ variations over the past 4000 years with the sensitivities of 0.72 (‰ ${}^{\circ }$C${}^{-\mathrm{1}})$ and $-$0.58 (‰ (0.01 m ice year${}^{-\mathrm{1}}{)}^{-\mathrm{1}})$. To understand the processes of the $\mathit{\delta }$Ar$/$N${}_{\mathrm{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 $\mathit{\delta }$Ar$/$N${}_{\mathrm{2}}$ of the microbubbles is negatively correlated with changes in accumulation rates through increases in the overloading pressures, although it underestimates the magnitude observed in the GISP2 ice core. Colder temperatures are found to induce more depletions in $\mathit{\delta }$Ar$/$N${}_{\mathrm{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.

Acknowledgements

We are grateful to G. Hargreaves at the US National Ice Core Laboratory, J. P. Steffensen at the Center for Ice and Climate, S. Kipfstuhl at the Alfred Wegener Institute, and S. Fujita at the National Institute of Polar Research for supplying ice core information. We thank T. Uchida, and B. Vinther for discussions, and R. Spahni for help on firn modeling. This project is supported by KAKENHI 23710020, 25740007, 22221002, 21221002, and 21671001, and a EU Marie Curie fellowship for T. Kobashi. The NGRIP and Dome Fuji ice cores were analyzed at the National Institute of Polar Research by T. Kobashi supported by the Senshin project and NIPR ice core center. The GISP2 ice core was analyzed by T. Kobashi at Scripps Institution of Oceanography, supported by J. Severinghaus. The early version of the paper was written when T. Kobashi was visiting the Centre for Ice and Climate (CIC), University of Copenhagen, in spring 2015 hosted by T. Blunier and B. Vinther. Finally, we thank two anonymous reviewers, whose comments substantially improved this manuscript. Edited by: W. T. Sturges