Atmospheric Chemistry and Physics Controls on the movement and composition of firn air at the West Antarctic Ice Sheet Divide

. We sampled interstitial air from the perennial snowpack (ﬁrn) at a site near the West Antarctic Ice Sheet Divide (WAIS-D) and analyzed the air samples for a wide variety of gas species and their isotopes. We ﬁnd limited convective inﬂuence (1.4–5.2 m, depending on detection method) in the shallow ﬁrn, gravitational enrichment of heavy species throughout the diffusive column in general agreement with theoretical expectations, a ∼ 10 m thick lock-in zone beginning at ∼ 67 m, and a total ﬁrn thickness consistent with pre-dictions of Kaspers et al. (2004). Our modeling work shows that the air has an age spread (spectral width) of 4.8 yr for CO 2 at the ﬁrn-ice transition. We also ﬁnd that advection of ﬁrn air due to the 22 cm yr − 1 ice-equivalent accumulation rate has a minor impact on ﬁrn air composition, caus-ing changes that are comparable to other modeling uncertainties and intrinsic sample variability. Furthermore, estimates of 1 age (the gas age/ice age difference) at WAIS-D appear to be largely unaffected by bubble closure above the lock-in zone. Within the lock-in zone, small gas species and their isotopes show evidence of size-dependent fractionation due to permeation through the ice lattice with a size threshold of 0 . 36 nm, as at other sites. We also see an unequivocal and unprecedented signal of oxygen isotope fractionation within the lock-in zone, which we interpret as the mass-dependent expression of a size-dependent fractionation process.


Introduction
Air extracted from the perennial snowpack (firn) is scientifically interesting for two reasons: it provides insight into the processes that determine the composition of air, both in the firn, and in the bubbles trapped in the underlying ice, and it provides a record of atmospheric composition in the recent past.
In December of 2005, we traveled to the US West Antarctic Ice Sheet Divide (WAIS-D) research station, where a deep ice coring program was commencing. Working from a pair of shallow holes, we incrementally extracted samples of firn air from the surface down to the firn-ice transition. The samples were returned to various laboratories in the United States and analyzed for a range of different species. The primary goal of this work was to quantify the nature of firn air movement at the site, thereby aiding interpretations of trapped gases extracted from the deep core. In addition, the measurements provide atmospheric histories of intermediate length and resolution for a wide range of species, allowing comparisons with existing records from other sites and emissions models.
The utility of the firn as an archive of past atmospheres was first realized by Schwander et al. (1993). Many subsequent studies have used firn air to establish or confirm recent atmospheric histories of O 2 , CO 2 , CH 4 and N 2 O, some of their isotopes, and numerous halocarbons and other trace species (e.g. Bender et al., 1994;Battle et al., 1996;Trudinger et al., 1997;Butler et al., 1999;Braunlich et al., 2001;Sturges et al., 2001;Montzka et al., 2004). The firn-based histories vary in length and resolution, depending on the local climate and the resulting glaciological properties of the firn from which the air was extracted (Kaspers et al., 2004). With occasional exceptions (Severinghaus et al., 2010), sites with low temperatures and high accumulation rates generally have the longest records. At each level in the firn, the air is a mixture of air of different ages, due to the gas diffusion process. Thus, the longest records also tend to have the widest age spreads (Trudinger et al., 1997), thereby limiting the time resolution with which past atmospheric histories may be reconstructed.
In the manuscript that follows, we focus on the use of gas measurements to illuminate the physical processes controlling the composition of air in the firn and in bubbles in the underlying ice. We begin with an overview of air movement in the firn and then describe the sampling site, collection methods, and modeling approach. We then use measurements of multiple atmospheric species and their isotopes to quantify the presence and relative importance of various transport and fractionation processes within the open pores of the firn. Finally, we consider fractionation of gases during the formation of bubbles, and briefly discuss the relation between the formation of bubbles and determination of the gas age/ice age difference ( age).

Air movement in firn
Air movement in the firn is driven by many processes on a variety of spatial and temporal scales. Concentration gradients within the firn drive gas fluxes via molecular diffusion (Craig et al., 1988). These concentration gradients can arise from changes in the composition of the overlying atmosphere, gravitational settling of heavy species, thermally driven migration of heavy species to cold strata within the firn (Severinghaus et al., 2001), and close-off related fluxes of relatively mobile species (Severinghaus and Battle, 2006).
Changes in barometric pressure and windflow over surface contours can also lead to bulk convective transport in the shallow firn (Colbeck, 1989;Sowers et al., 1992). Significant temperature gradients have the potential to cause convective instability and vertical overturning (Powers et al., 1985;Severinghaus et al., 2010), while the steady accumulation of snow results in downward motion of the firn and its air bubbles, inducing a downward motion of the air (Rommelaere et al., 1997). Finally, as air pockets are compressed and firn crystals sinter into ice, permeation may occur as gases escape through the crystal lattice in the newly forming ice bubbles (Severinghaus and Battle, 2006).
While advection and convection move all species through the firn indiscriminately, diffusion and permeation have the potential to fractionate gas mixtures. Different processes dominate in different parts of the firn (Sowers et al., 1992). As indicated schematically in Fig. 1, convection is confined to the uppermost region in the firn. Below this is an extensive zone where diffusion dominates. In the lower firn, many  showing generic features (on the right) and the transport mechanisms expected to be significant in each region (on the left). Also shown (in parentheses on the right) are the approximate depths of the features at WAIS-D. The figure is not drawn to scale, and not all processes will be present in every naturally occurring firn column.
sites (including WAIS-D) exhibit a lock-in zone, where vertical air movement is largely impeded by low-permeability seasonal layers, yet air samples can still be extracted from the remaining open pores in the intervening layers laid down in the complementary season (Martinerie et al., 1992). The firn column terminates at the firn-ice transition, where the last of the open pores fuse shut, forming ice bubbles. At the firn-ice transition, or immediately above it, the severely limited open porosity precludes the extraction of firn-air samples.

Site characteristics and significance
We collected air samples from two holes (WDC05-B and WDC05-C), each 10.8 cm in diameter, separated by roughly 20 m, drilled at 79 • 27.78 S 112 • 07.39 W, at an altitude of 1766 m (Conway and Rasmussen, 2009). This site was 2 km upwind of the main camp and deep core location, in a previously-undisturbed region. Mean annual temperature at the site, measured in the borehole, is −31 • C with an accumulation rate of 22 cm yr −1 ice-equivalent (Banta et al., 2008). These characterisitics make the site an intermediate one in many respects. It is neither very cold, nor is it warm enough to have melting episodes. The accumulation rate is comparable to Summit Greenland (∼ 25 cm yr −1 ) and Law Dome (∼ 16 cm yr −1 at DSSW20K) (Smith et al., 2000), higher than South Pole and Vostok (< 10 cm yr −1 ), but significantly lower than other Law Dome sites (65 cm yr −1 at DSS and 1.19 m yr −1 at DE08) (Etheridge et al., 1996).

Sample collection and analysis
Drilling and air-sampling occurred between 21 and 29 December 2005, generally following standard protocols described by Butler et al. (1999). At every sampling depth we collected pairs of PTFE-sealed glass flasks for analysis by the Carbon Cycle Greenhouse Gas group of NOAA and the Scripps Institution of Oceanography, pairs of PTFE-sealed glass flasks for Penn. State University, single pressurized PTFE-sealed glass flasks and single pressurized all-stainless steel flasks for the NOAA HATS group, and paired PTFEsealed glass flasks for the University of California, Irvine group. A full list of the species analyzed and the responsible institutions are given in Table 1. All samples were collected through Dekoron ™ tubing using a Metal Bellows MB-158 pump (Senior Aerospace, Sharon MA), with a KNF-NO5 pump with Viton valves and diaphragm (KNF Neuberger, Trenton NJ) added in series for the UC Irvine and NOAA HATS flasks. As described by Severinghaus et al. (2010), we used a system in which samples were collected through a continuous piece of tubing, which was free of couplings and joints from the intake at the bottom of the borehole to the low-pressure side of our sampling pump.
Samples were collected at complementary depths in the two holes, yielding a dataset with collections every 5 m in the diffusive column and roughly every 1 m in the lock-in zone. In addition, we collected two sets of surface samples (through the firn sampling equipment) and another set from each hole at a common depth (nominally 70 m). As in nearly all (Buizert et al., 2011) previous campaigns, compositions of the samples from the two holes agreed within uncertainties.

Modeling of firn-air movement and composition
To interpret our measurements of firn air, we employ a 1dimensional model built around a finite-difference approximation of Fick's 2nd law. We have three versions of the model. All versions operate with a 4.8 h time step and include molecular diffusion and (if desired) gravitational fractionation. The first, described previously by Battle et al. (1996), uses a fixed coordinate system and divides the firn column into 1 m thick boxes. The second follows the approach of Trudinger et al. (1997), and uses a moving coordinate system with boxes of equal ice mass, ranging in thickness from 1.03 m at the surface to 0.53 m at the firn-ice transition. The movingcoordinate version allows for convection in the shallow firn using an exponentially damped eddy diffusivity term  and implicitly accounts for downward advection of firn air due to the annual accumulation of snow at the surface. In every other way, it is identical to the first.
The third version of the model, similar to that used by Severinghaus et al. (2010), incorporates thermal fractionation and heat transport. This one includes the same formulation of convection as the second model, and also includes downward advection, but as an explicit process, set within a fixed coordinate system in the open firn. Within the lockin zone, advection is treated by shifting concentrations down a temporally-uniform grid to minimize numerical diffusion (Severinghaus et al., 2010).
None of these models include the bubble overpressure that drives permeation. Except where noted, all model results quoted below are derived from the second version.

Large-scale regimes in the firn
Species with temporally and spatially constant atmospheric mixing ratios can be powerful tools for investigating air movement in the firn, if they can be measured with sufficient precision. We begin by using δ 15 N of N 2 , shown in Fig. 2, to identify the various regimes depicted schematically in Fig. 1. While exact depth ranges for each regime are difficult to determine, approximate values (justified in Sects. 6.1 through 8) are indicated in Fig. 1. These data show many of the qualitative and quantitative features we expect.

The convective zone
The thickness of the convective zone at the top of the firn column can be challenging to identify. One approach is to extrapolate a linear fit to the data clearly in the diffusive zone and below the range of thermal influence (see Sects. 6.2 and 8) back to the depth at which δ 15 N = 0. Using only NOAA flasks between 34 and 56 m depth, this method indicates a convective zone of 5.2±1.6 m. A second approach, known as the barometric slope method (Bender et al., 1994;Severinghaus et al., 2010), fixes the slope of the line at the theoretical prediction (4.88 per meg m −1 ) and adjusts the intercept until the line and the deeper observations optimally agree (Bender et al., 1994). Again using only the NOAA flasks between 34 and 56 m, this method yields a convective zone 3.1 ± 0.4 m deep (assuming no uncertainty in the theoretically based slope). These two approaches are shown in Fig. 2. A third approach involves adjusting an eddy diffusivity term in the model until it optimally matches the observations in the diffusive column. The convective zone thickness is then defined by the depth at which eddy diffusivity equals molecular diffusivity . Using the second (i.e. moving coordinate) version of the model with the formulation of eddy diffusivity of Kawamura et al., we find a convective zone thickness of 1.4 m. A fourth approach (Severinghaus et al., 2010) again adjusts eddy diffusivity for optimal data-model agreement, but this time uses the version of the model that includes thermal fractionation. It also defines the convective zone slightly differently. Specifically, Values are reported relative to air collected in La Jolla, CA. Measurements made on sets of flasks from NOAA and Penn. State University (see text) showed no systematic differences. The inset simply enlarges the scale on those collections made in, and immediately above, the lock-in zone. Model results are from the third model, incorporating thermal gradients and advection, as well as gravitation and eddy diffusion (see Sect. 5). Also shown are linear fits to the NOAA flasks collected between 34 and 56 m, with the slope fixed to the theoretical value, and with the slope as a free parameter. the thickness of the convective zone is determined from the difference between the depths of a discrete set of model δ 15 N values located clearly within the diffusive zone, modeled with and without convection. This fourth approach yields 2.1 m as the convective zone thickness.
This range of results probably reflects the fact that each is defining the thickness of the convective zone in slightly different ways and only the last method includes thermal fractionation. Regardless of the method chosen, it seems clear that the convective zone is shallow, and bulk air motion originating at the surface has little influence on the composition of air deeper in the firn.

The diffusive zone
Considering only the NOAA CCGG flasks, samples collected at 5 unique depths between 34 and 56 m exhibit a progressive enrichment in the heavy isotope, increasing linearly at 5.2 ± 0.2 per meg m −1 , (1σ errors) where 1 per meg = 0.001 ‰ = 0.0001 %. Including the samples collected at 30 m and/or those between 56 and 65 m does not siginificantly change this result. Similary, including the PSU flasks has no significant impact on the measured enrichment, despite the fact that much of the gas in these flasks was removed for other analyses prior to sampling for δ 15 N. The slope of δ 18 O O 2 data over this same depth interval is about twice as large: 10.2 ± 0.4 per meg m −1 . This is consistent with the δ 15 N enrichment on a per-mass-unit basis. These values are not completely consistent with the theoretical prediction of 4.88 meg mass unit −1 m −1 for pure gravitational/diffusive equilibrium given by the barometric equation (Craig et al., 1988) using a temperature of −31 • C and a local gravitational acceleration of 9.825 m s −2 (Lide and Weast, 1985). We believe that the discrepancy is probably due to the influence of seasonal thermal signals extending more than 30 m into the firn (see Sect. 8).

The lock-in zone
As shown in the inset of Fig. 2, δ 15 N increases slightly between 66.84 and 68.26 m, but at far less than the rate expected from gravitational enrichment. Below 68.26 m, δ 15 N is either constant or decreasing. The lack of fully realized gravitational enrichment indicates the bottom of the diffusive column and beginning of the lock-in zone (also known as the "sealing depth", Etheridge et al., 1996). The exact onset of the lock-in is difficult to pinpoint without more precise and abundant δ 15 N measurements. Data in this region may also reflect collection fractionation artifacts (see Sect. 10).

Firn-ice transition
The lock-in zone extends down to 76.54 m (the "total closeoff depth") , where the firn-ice transition precludes collection of any more samples. This depth is consistent with the modeling work of Kaspers et al. (2004). Using the densification model of Herron and Langway (1980), Kaspers et al. predicted a total close-off depth of ∼ 75 m for a site with the temperature and accumulation rate of WAIS-D (where close-off density is 803 kg m −3 ). The agreement provides a valuable check of the Kaspers et al. model for a previously untested combination of temperature and accumulation rate.

Advection due to accumulation
In addition to gravitation and convection, firn air will be advected downward. Trudinger et al. (1997) showed this is an important effect at high accumulation sites, while Battle et al. (1996) found it negligible at a low-accumulation site. In simplest terms, we seek to answer the question "Can a firn model reproduce the data at WAIS-D, even if it ignores advection?" Because we don't really know the diffusivity-depth relationship in the firn, this question can be restated "Does the freedom to choose a diffusivity profile enable us to compensate for the absence of advection in our model?" To answer this question, we use the first and second versions of our model (see Sect. 5). We infer a diffusivity profile through a tuning process that begins with density measurements of the firn core made at the time of collection. We then used these density data and the density-diffusivity relationship of Schwander et al. (1988) (Eq. 2) to create a preliminary diffusivity profile. This profile was then iteratively adjusted (piecewise-continuously) to improve the data-model agreement for CO 2 . Finally, the tuning was checked with CH 4 and F11. In each case, the model uses a diffusivity-depth profile that was tuned using the model itself (rather than a single profile for both models).
Comparisons of paired model runs for two species with very different atmospheric histories and molecular diffusivities (CO 2 and HFC-134a) are shown in Fig. 3. Also shown are measured values of these species. At some depths, there are clear differences between the model results with and without advection. These differences are generally comparable to (or larger than) the uncertainty in our mixing ratio measurements at any given depth. On the other hand, at some depths, single-depth measurement precision and advection are both small compared to other influences and uncertainties. For example, around the top of the lock-in zone, the exact location and sharpness of the diffusive/lock-in transition (64-68 m at this site) clearly have a greater impact on the data-model agreement than the presence or absence of advection. Similarly, in the lower part of the diffusive column (45-65 m here), the data-model differences are comparable to the model-model differences. These data-model discrepancies may be due to problems with our tuning of the diffusivitydepth profile (common to both models), or perhaps with the atmospheric forcing history we are using. The same comment applies for CO 2 in the lock-in zone, where additional "roughness" in the CO 2 -depth profile points to real variability in the firn air that is large compared to advective effects.
Uncertainties associated with the tuning of the diffusivitydepth profile can be reduced by more sophisticated multispecies techniques (Buizert et al., 2011), but real variability in the composition of firn air extracted from the lock-in zone is a more difficult problem.
In short, this exercise shows that the impact of advection at an intermediate-accumulation site like WAIS-D is likely smaller than, or comparable to, other apects of our models which are poorly constrained. Furthermore, as models of firn processes improve, as more precise datasets are assembled, or when working at sites with substantially higher accumulation rates, it might become necessary to take advection into account. Each model was run with its own diffusivity-depth profile. Two species with very different atmospheric histories (Montzka et al., 1996;Etheridge et al., 1996;Conway et al., 2010) were chosen for illustrative purposes. Also shown are our measurements of the firn air, with error bars. The error bars are based on the pair-wise agreement of replicate flasks collected at each depth, and are usually smaller than the markers on the plot.
here, with our noble gas data showing responses to thermal gradients that are qualitatively similar to other sites (Fig. 2). A quantitative analysis of thermal signals in the WAIS-D firn air will be presented elsewhere (Orsi et al., 2011).

Movement through the lattice and size-dependent fractionation
As discussed by Severinghaus and Battle (2006) and Huber et al. (2006), some species of gases appear to escape from recently-closed bubbles into the open pores immediately adjacent. It is hypothesized that this permeation leakage occurs through the ice lattice itself (Ikeda-Fukazawa et al., 2005). A molecular diameter of 0.36 nm appears to be a threshold value, with larger species contained (or all escaping at roughly the same slow rate), and smaller species escaping at rates that depend on their molecular size. The high mobility of O 2 , relative to N 2 , as predicted by Ikeda-Fukazawa et al. (2004), is corroborated by the observed decline over time of O 2 /N 2 in ice core samples stored in freezers at −25 • C and by the observed increase in O 2 /N 2 in the lock-in zones of virtually all firn air study sites (Severinghaus and Battle, 2006).
At WAIS-D, as at other sites, this size-dependent effect is evident in the substantial enrichment of O 2 /N 2 and Ne/N 2 above and within the lock-in zone, as well as the nonenrichment of Xe/N 2 and Kr/N 2 . We follow Severinghaus and Battle (2006) and plot various species versus O 2 /N 2 within the lock-in zone in Fig. 4. As in this earlier work, we make O 2 /N 2 behave as if it were temporally constant at the surface of the firn by subtracting a modeled component of O 2 /N 2 variation with depth (an O 2 /N 2 "correction"). This correction was generated by forcing the model with an atmospheric history of O 2 /N 2 derived from direct measurements at the South Pole since 1992 (Manning and Keeling, 2006), and in earlier times, fossil fuel combustion records (Marland et al., 2007) and estimates of biospheric carbon release (Bruno and Joos, 1997). While the uncertainty in this correction is hard to fully quantify, model output from differing scenarios gives some sense of the magnitude of the potential error. δO 2 /N 2 values in the deep firn for model runs assuming a balanced biosphere are ∼ 19 per meg higher than those using our best estimate of biospheric behavior. Even if this extreme range of input assumptions really does indicate the uncertainty in our correction, it is small compared to the total correction applied (roughly 1 ‰) and smaller still compared to the range of δO 2 /N 2 observed and modeled within the lock-in zone.
As at other sites (Severinghaus and Battle, 2006;Huber et al., 2006), Kr/N 2 and Xe/N 2 show essentially no covariation with O 2 /N 2 . Slopes for these species are 0.04±0.03 and 0.16±0.07 ‰/‰ respectively (Fig. 4 inset). We note that the exact values we measure for these slopes depend upon the details of the correction for gravitational and pressure-gradient fractionation we apply to the raw data (see Sect. 10 and Severinghaus and Battle, 2006). This dependence is particularly acute for Kr/N 2 and Xe/N 2 , on account of their large mass differences.
In contrast, Ar/N 2 shows a small but unequivocal covariation (0.27 ± 0.01), while Ne/N 2 shows a strong signal (24.7 ± 0.4). Interestingly, these slopes in the WAIS-D lock-in zone differ from the slopes measured at South Pole (0.15 ± 0.07 and 34.0 ± 1.5 for Ar/N 2 and Ne/N 2 , respectively) (Severinghaus and Battle, 2006). All uncertainties given here for slopes vs. O 2 /N 2 are 95 % confidence intervals.
Why should Ar/N 2 and Ne/N 2 have different slopes at WAIS-D and the South pole? The explanation may lie in the temperature-dependence of the mechanisms by which Ne and O 2 move through the lattice. Ikeda-Fukazawa et al. (2002 suggested that the movement of Ne in ice occurs by way of "hops" between interstitial sites. In contrast, O 2 moves primarily (though not exclusively) through a mechanism in which hydrogen bonds in the lattice are broken as the mobile impurity (in this case O 2 ) migrates (Ikeda-Fukazawa et al., 2004). Ar, being closer to the 0.36 nm threshold than O 2 , is still more dependent on bond-breaking. activitation energy, R is the gas constant, and T is the temperature. However, O 2 and Ar are expected to have a larger activitation energy Q than Ne, due to the strong temperature dependence of the breaking of hydrogen bonds. Thus, the warmer temperature at WAIS-D relative to South Pole (−31 vs. −51 • C) probably explains the lower observed slope of Ne/N 2 vs. O 2 /N 2 and the slightly higher slope of Ar/N 2 vs. O 2 /N 2 . This is also consistent with the relatively large slope of Ar/N 2 vs. O 2 /N 2 at Siple Dome: Slope = 0.33 ± 0.02, temperature = −25 • C (Severinghaus and Battle, 2006).

Anomalous fractionation in the lock-in zone:
The mobility of Ne, O 2 , Ar, Xe and Kr (relative to N 2 ) can be understood in terms of a size-dependent mechanism. However, at WAIS-D, we also observe very clear evidence of isotopic fractionation in the firn-air record of δ 18 O O 2 . In particular, there is a pronounced depletion of 18 O 16 O with depth in the lock-in zone (Fig. 5). As discussed in Sect. 10, δ 15 N is believed to vary in the lock-in zone solely due to pressuregradient fractionation at the time of collection. The possible departure from linearity in Fig. 5  Clearly, there must be some other mechanism at work. One possible explanation is that isotopic fractionation occurs during the permeation of O 2 through the ice lattice. The large molecules N 2 , Kr and Xe appear to be unfractionated by permeation despite large mass differences (Severinghaus and Battle, 2006). However, as discussed above, O 2 is smaller and therefore may move through the ice lattice by both bond-breaking (emphasized in Sect. 9.1 above) and the interstitial hop mechanism. The modeled 3-fold greater permeation rate of O 2 relative to N 2 reflects this contribution from hopping (Ikeda-Fukazawa et al., 2004;Severinghaus and Battle, 2006). The hops should be velocity-dependent, so a fractionation related to √ mass might be expected, and with it, a trend in δ 18 O O 2 : the mass-dependent expression of a size-dependent process.
Here, we explore this possibility by calculating the relative permeation coefficients implied by the observed trends of δ 18 O O 2 and O 2 /N 2 in the lock-in zone. Unfortunately, it is not straightforward to infer these physical parameters from the raw data because N 2 is also permeating at the same time as O 2 , effectively adding to the denominator as well as the numerator of the O 2 /N 2 ratio. We account for this with the following assumptions: 5. the composition of the bubbles does not evolve significantly during the course of the evolution of the gas in the open pores, so that the bubbles may be treated as having a constant composition.
Using these assumptions, our measurements imply the ratio of permeabilities of 18 O 16 O and 16 O 16 O is 0.9942 ± 0.0005, corresponding to a fractionation factor of −5.8 ± 0.5 ‰ (see Appendix B). This is much less than expected for a simple velocity-dependent process, for which we would expect a fractionation of √ mass = √ 32/34 = −30 ‰. We conclude that the processes occurring in the lock-in zone are more complex than a simple velocity-dependent transport. Our result also reinforces the view that effusion cannot be the dominant process in close-off fractionation because effusion also predicts a √ mass dependence (Ikeda-Fukazawa et al., 2005).
To our knowledge this is the first reported instance of natural isotopic fractionation by the bubble close-off process. This observation was made possible in part by significant advances in analytic precision (see Appendix A). If corroborated, it has implications for isotopic reconstructions of atmospheric gases from both firn air and bubble air in ice cores. Some prior studies have already attempted to make empirical corrections to measured δ 18 O O 2 of O 2 in ice core records using co-measured δO 2 /N 2 and Ar/N 2 (Severinghaus et al., 2009;Landais et al., 2010). While most of this fractionation is an artifact that occurred due to gas loss during core retrieval and handling (e.g. Bender et al., 1995), some small component of this fractionation may have been due to natural close-off fractionation. Our result supports the use of this empirical correction.
Given that oxygen isotopes are fractionated as the (indirect) result of size-dependent O 2 /N 2 fractionation, one might expect that argon isotopes might also be fractionated as the result of size-dependent Ar/N 2 fractionation. However, this does not appear to be the case at WAIS-D, although this conclusion is more tenuous than for δ 15 N.
Ar/N 2 is modestly enriched in the firn air (Fig. 4), suggesting that Ar is indeed expelled preferentially during bubble close-off like O 2 /N 2 , but more weakly. The fact that the atomic diameter of Ar is less than the threshold value of 0.36 nm, but more than that of O 2 , is consistent with the observed Ar/N 2 and O 2 /N 2 data. Nonetheless, there is no clear sign of isotopic fractionation in Ar: compare Fig. 7 and Fig. 6.
It is conceivable that the physics of the O 2 isotopic fractionation are not operative for the isotopes of a monatomic noble gas; in other words, mass-dependence might only exist for diatomic molecules that have the zero-point energy changes associated with covalent bonds. This speculation should be explored further. We also note that while δ 40 Ar/ 36 Ar is known to be fractionated by artifactual gas loss from ice core samples (Severinghaus et al., 2003), it is not necessarily the case that artifactual loss involves the same physical processes as natural close-off fractionation.
An obvious next step to test our hypothesis of isotopic close-off fractionation would be to measure neon isotopes in firn air, as Ne/N 2 is known to be heavily impacted by closeoff fractionation (Fig. 4 and Severinghaus and Battle, 2006). In addition, Ne has no significant atmospheric change, simplifying the interpretation. Because the small atoms neon and helium are thought to permeate through the ice lattice by the interstitial hopping mechanism, we expect that there will be substantial isotopic fractionation. We likewise anticipate that hydrogen isotopes of H 2 will be affected by the close-off fractionation, requiring a correction to atmospheric reconstructions, because the collision diameter of H 2 is similar to that of neon (Reid et al., 1987).

Collection fractionation at depth
In contrast to O 2 just discussed, there is no evidence that istopes of N 2 are fractionated during bubble close-off (Severinghaus and Battle, 2006). This is consistent with the hypothesis of Ikeda-Fukazawa et al. (2005) that hydrogenbond breaking limits the rate of diffusion through the ice lattice of N 2 and other gases with molecular diameters greater than 0.36 nm (Severinghaus and Battle, 2006;Huber et al., 2006). In addition, this hydrogen-bond breaking mechanism is insensitive to the mass of the molecule, implying that Kr and Xe and N 2 should not be fractionated with respect to each other despite large mass differences. This is consistent with observations at other firn sites (Severinghaus and Battle, 2006). As discussed by Severinghaus et al. (2009), there are also several other lines of evidence that N 2 is not fractionated during bubble close-off. Despite the lack of evidence for isotopic fractionation of N 2 during close-off, there is an unequivocal depthdependence of δ 15 N within the lock-in zone at WAIS-D (Fig. 2), as well as at other sites (Severinghaus and Battle, 2006). This δ 15 N signal has been attributed to the act of extracting air from the open porosity during sampling. When open channels are very restricted (as they are in the lock-in zone) the suction required to extract air may lead to substantial pressure gradients within the firn, accompanied by mass-dependent pressure-gradient fractionation (Chapman and Cowling, 1970). Such a collection artifact was invoked by Severinghaus and Battle (2006) to explain the drops in δ 15 N, δKr/N 2 and δXe/N 2 observed within the lock-in zone at both Siple Dome and South Pole.
While the absolute decrease in δ 15 N in the WAIS-D dataset is difficult to quantify, if pressure-gradient fractionation is present, we expect noble gases to show heavy isotope depletion in the lock-in zone, covarying with δ 15 N on a 1:1 line when normalized by mass difference. The same should be true of elemental ratios (Severinghaus and Battle, 2006).
Unfortunately, our measurements are too few and too imprecise to make a definitive statement about the presence (or absence) of pressure-gradient fractionation at WAIS-D. In particular, slopes fitted to our data are unduly influenced by the values measured at the single deepest collection point.
It is worth noting that pressure-gradient fractionation may not be present in all datasets (e.g. Huber et al., 2006) the development of pressure-gradient fractionation will depend on the specifics of flow rates and the geometry of the sampling equipment where air is being extracted from the firn. That said, regardless of whether it is weakly or strongly realized, pressure-gradient fractionation should manifest itself as a 1:1 covariation in species-vs.-δ 15 N plots.

Age distributions of firn air
In addition to the species presented here, the WAIS-D firn air samples were analyzed for a wide range of hydrocarbons, halocarbons and other geochemically significant species. Some records have already been published (methane isotopes (Mischler et al., 2009), HFC-23 , CFC-12  and ethane (Aydin et al., 2011)), while publication is forthcoming for others such as methane (Mitchell et al., 2011). Common to these analyses, and important for linking firn air records to underlying ice-core records, are the age distributions of the firn air at various depths. As discussed in detail by Schwander et al. (1993) and Trudinger et al. (2002), the composition of the air extracted from the firn at any given depth reflects the convolution of the atmospheric history with diffusive mixing subject to the depth-dependent diffusivity of the firn column itelf. This statement is true only to the extent that fractionating processes discussed above are either negligibly small, or have been corrected for. Assuming this is the case, we can follow earlier authors (Schwander et al., 1993;Rommelaere et al., 1997;Trudinger et al., 2003) and use one-year, normalized pulse-response functions from our 1-dimensional moving-coordinate firn model (described above) to characterize the age of the firn air as a function of depth. Since each type of molecule diffuses at a different speed, distributions are species-specific. Results of this exercise for CO 2 are shown in Fig. 8 and are consistent with other sites with similar glaciology (e.g. Trudinger et al., 2003). Air from shallow depths is young with a sharp peak, reflecting the brief time over which diffusion has been acting. As depth increases, the air becomes older, with a broader age spread. Once within the lock-in zone, the modeled air ages at the rate that would be expected from the accumulation rate. In reality, the firn air within the lock-in may age slightly less rapidly than the surrounding ice due to expulsion of air driven by compaction of the firn. This expulsion is neglected in our models. The coarse nature of the model curves within the lock-in zone in Fig. 8 has no physical significance and simply reflects the discretization chosen for the model treatment of the lock-in.
Quantifying the age of the air is challenging since various metrics may be chosen to characterize these non-Gaussian age distributions (Trudinger et al., 2002). For comparison with other sites, we calculate the 2nd moment, or "spectral width" of the age distribution at close-off, defined by the equation where z is depth, is the mean of the distribution, and G is the age distribution (Hall and Plumb, 1994;Hall and Waugh, 1997;Andrews et al., 1999). Using this metric for CO 2 , we find = 4.8 yr (with = 45.5 yr at 75.2 m; close to the bottom of the lockin zone) at WAIS-D. This compares to = 5 and = 17 yr for Law Dome (DSSW20K) and South Pole, respectively (Trudinger et al., 2002). Given the similarities in glaciology of DSSW20K and WAIS-D, similar spectral widths are to be expected.

Age of air in ice bubbles
Thus far, we have focussed on the composition of air extracted from open pores in the firn. Of equal or greater interest is the composition of air trapped in closed pores, since this bubble air is the basis for so much of our current understanding of paleoclimate. One quantity of importance is the difference in age between bubble air and the surrounding ice ( age). Mischler et al. (2009) calculated age = 205 yr for WAIS-D, assuming that essentially all bubbles close within the lockin zone. However, the exact location of bubble close-off remains somewhat uncertain (Schwander, 1989 Fig. 8. Normalized age distributions for CO 2 in firn air at selected depths (given in the legend) from the moving-coordinate firn model (see text). The step-like structure within the lockin zone is a result of model discretization and does not reflect any underlying physics. The inset shows the mean age as a function of depth, with points indicating the means for the particular distributions selected for the main plot. The dashed line is included only to highlight the departure from linearity. and Battle, 2006;Aydin et al., 2010). If a modest fraction of bubbles actually close within the diffusive column, the age of bubble air derived using these assumptions will be biased young. In addition, the temporal resolution of the record from bubble air will be reduced. We can attempt to understand the amount of bubble formation within the diffusive column by exploiting the steady gravitational enrichment of 15 N 14 N in this part of the firn. In particular, if substantial bubble formation is occurring within the diffusive column, bubble air from below the firn-ice transition should have a somewhat lower δ 15 N than the firn air at the top of the lock-in zone (deeper in the lock-in zone, δ 15 N may be reduced by the pressure-gradient fractionation associated with collection and described above).
The two firn air samples from 68.26 m depth have δ 15 N = 0.3086 ‰ and 0.3057 ‰ (or 0.307±0.002). The 130 bubbleair measurements made between 129 and 256 m yield δ 15 N = 0.3019 ± 0.0006 ‰ (1 standard error) on the same analytic scale. The agreement seems to imply, at most, a modest amount of bubble closure above the lock-in zone.
To make this more quantitative, we consider the parameterizations of bubble closure presented by Schwander (1989, Eq. 3) and Severinghaus and Battle (2006, Eq. A9). These relations, along with densities measured on the WAIS-D firn cores in the field, are used as input for a closed-pore model that we include in the firn diffusion model described above. This bubble trapping in each box of the model begins with the assumption that air in the open pores is at atmospheric pressure. For the purposes of this analysis, we also assume that no fractionation occurs during trapping. As each box moves downwards, some of the open porosity closes (the amount is an input to the model), with these new bubbles containing air of the same composition as the immediate (still open) surroundings. Consequently, the mixing ratio of all closed bubbles at a given depth will be the weighted average of the new bubble air mixing ratio and the preexisting bubble air mixing ratio (where the weights are the respective molar abundances of the newly-and previously-trapped air).
The Schwander and Severinghaus parameterizations of bubble closure imply only 8 % and 5 % of pores have closed above the lock-in zone, respectively. Using these as inputs to the model, we find a maximum difference in δ 15 N between open and closed pores of 0.9 and 0.5 per meg, respectively; too small to detect, given the scatter in our firn air (openpore) measurements.
As a simple sensitivity study, we can calculate how much closed porosity would have to form above the lock-in zone in order to give a measureable difference between firn air and bubble air in δ 15 N. We start with a parameterization of bubble formation and then simply scale the amount of porosity that closes above the lock-in zone by a multiplicative factor. We then adjust this factor until the model predicts a significant difference between δ 15 N of air in open and closed pores. We define a significant difference as non-zero at the 99 % confidence level, or about 5.4 per meg, where this value comes from the uncertainties in our measurements of both the bubble air and firn air. We find that unless at least 40 % of the pores close above the lock-in zone, the difference between the firn and bubble air δ 15 N values will be essentially undetectable.
Naturally, this result depends strongly on the location of bubble formation above the lock-in zone. Both the Severinghaus and Schwander parameterizations place the majority of pore closure quite close to the lock-in zone. If instead, the closure is occurring farther up the column, a smaller fraction of the total pores closing will yield a discernable signal in δ 15 N.
Despite these caveats, the calculations above make two things clear. First, δ 15 N provides a rather weak constraint on the amount of bubble formation occurring above the lockin zone. In other words, only if a substantial amount of pore closure is occurring above the lock-in zone, will we be able to see signs of it in δ 15 N.
Secondly, assuming that all bubble formation occurs in the lock-in zone probably does not lead to a significant error in age. This is because the age of the firn air increases with depth in a non-linear fashion in the diffusive column (see Fig. 8, inset). In particular, the strongest dependence occurs closest to the lock-in zone, where diffusivity is decreasing most rapidly. This is in contrast to δ 15 N, which increases with depth at a nearly constant rate, set by the barometric equation (Craig et al., 1988). In addition, the firn air-bubble air δ 15 N difference is relatively insensitive to up-hole bubble formation. Taken together, these factors mean that a value of age determined simply from the difference between the age of the firn air from the lock-in zone and the age of the surrounding ice is likely to be essentially unbiased. This will be true even in the presence of bubble formation well above the depths predicted by Schwander and Severinghaus.

Conclusions
In many ways, the firn air sampling and analysis presented here is a confirmation of firn air work conducted at other sites. We have used established sampling methods and models that are similar to those already developed, at a site that has neither unusually high, nor unusually low accumulation rates and temperatures. The absence of significant convective influence, the thermal signal in the shallow firn, the expected gravitational enrichment in the diffusive column, the depth of the firn-ice transition, the presence and extent of a lock-in zone and the limited impact of accumulation-driven advection are all consistent with results at other sites. Taken together, these speak to the robustness of our understanding of processes operating in the firn, and the integrity of the dataset presented here. In addition, using δ 15 N data from firn air as well as ice bubbles, along with a simple model of bubble formation, we find that estimates of age (the gas age/ice age difference) are unlikely to be biased by the formation of bubbles above the lock-in zone.
Within the lock-in zone, we also find signs of processes that have been observed elsewhere. In particular, molecules that are smaller than 0.36 nm in diameter (such as Ne, O 2 and Ar) are enriched in the open pore space, with greater enrichment corresponding to smaller size. Like past authors, we attribute this phenomenon to permeation through the ice lattice. Quantitatively, the enrichments we see differ from those observed at South Pole, probably due to the relatively warmer temperatures at our current study site.
One unprecedented signal in our dataset is the clear variation of δ 18 O O 2 within the lock-in zone. This cannot be explained purely in terms of direct size-dependent permeation, but is at least conceptually consistent with a velocitydependence (and by extension, mass-dependence) in the interstitial hopping that seems to be occurring during the sizedependent permeation of O 2 through the ice lattice. If this explanation is correct, we might expect to see similar variation in the lock-in zone of Ne and Ar isotopes. Unfortunately, we do not yet have Ne isotopic analyses for these samples (in which we would expect a very large signal), and our Ar isotopic data are equivocally inconsistent with the simplest hypothesized mechanism. We look forward to future datasets that will shed light on the processes that give rise to the clear δ 18 O O 2 signal we see here.

Appendix A Correction for mass spectrometer "Chemical Slope"
It has long been recognized that the mass spectrometric measurement of an isotopic ratio in a mixed gas such as air can be subject to an artifactual sensitivity of the measurement to the major gas abundance ratios (Sowers et al., 1989;Severinghaus et al., 2003). The cause of the artifact is not well understood but probably involves ion-molecule reactions in the mass spectrometer source, which lead to small differences in the ionization efficiencies of the isotopologues. For δ 18 O of O 2 measurements, the N 2 abundance in the sample is the controlling variable for this effect.
A correction for this effect is routinely applied to isotopic measurements made in mixed gases, known as the "chemical slope" correction. The basis for the correction is that the mass spectrometer can be demonstrated to respond in a stable and linear fashion to variations in the major gas ratio that are induced experimentally by progressive additions of a pure gas (e.g., commercially obtained N 2 ) to aliquots of an air standard. Because the air standard does not change, any observed isotopic variation must be artifactual, and the slope of the observed variation is recorded and used to make the correction to the subsequent set of unknown samples.
This correction is crucial to the problem of detecting isotopic fractionation of O 2 in the lock-in zone, because large changes in O 2 /N 2 are associated with the small isotopic signal that is sought. Prior observations of δ 18 O O 2 of O 2 trends in lock-in zones at various firn air study sites (e.g. South Pole, Siple Dome) were discounted because of the very large chemical slope correction and its associated uncertainty.
For the present study, we made a special effort to tune the Scripps mass spectrometer to minimize its chemical slope. This was accomplished by varying source parameters such as the ion extraction voltage. As a consequence the uncertainty in the chemical slope correction is small relative to the signal observed, and therefore we have confidence in the result. The measured slope is 8 × 10 −5 ‰ for each ‰ increase in the amount of N 2 added (Fig. A1). This results in a maximum chemical slope correction of 6 × 10 −4 ‰ This is much smaller than our analytical precision of 0.004 ‰.

Calculation of the physical fractionation factor
We start by gravitationally correcting the measured δ 18 O of O 2 (thermal fractionation is ignored) using the canonical formula δ 18 O O 2 gravcorr = δ 18 O O 2 measured − 2δ 15 N measured (Sowers et al., 1989). If there is any sampling artifact, this correction will suffice to remove it (Severinghaus and Battle, 2006). Note that sufficiency of this correction depends on the validity of assumption 2 in Sect. 9.2.  A1. Results of the "chemical slope" experiment used for correcting oxygen isotope data. Pure commercially-obtained N 2 was added to aliquots of an air standard. Note that N 2 is in the numerator in this plot. The resulting corrections applied to the data were 0.0006 ‰ or less.
As described in the main text, δ 18 O O 2 gravcorr decreases with increasing δO 2 /N 2 (Fig. 7). A geometric-mean regression gives an empirical slope m of −0.0090 ± 0.0007 ‰ per ‰ (95 % confidence interval). From the definition of δ, and considering the addition of a small number of moles of gas d 18 O to the initial pore space inventory 18  Here, R STD is the ratio in the respective standard used to report the delta value, in our case the free atmosphere in La Jolla. Because the initial gas composition in the pores is very close to atmospheric, the approximation used in the secondfrom-right term is appropriate.
Following Severinghaus and Battle (2006), the amount of permeating gas that is added can be written in terms of the permeability P , the mole fraction in the bubble x, the bubble overpressure p, the total area of wall A available for permeation, and the effective mean wall thickness z L through which permeation occurs. The permeation is assumed to occur for a small increment of time dt. As p, A, z L and dt are the same for all gases, they may be lumped together in the parameter γ . γ has units of Pa m 2 m −1 s,  Removal of gas from the open pore spaces occurs continuously within the lock-in zone via closure of additional bubbles and compaction with expulsion of air to higher layers in the firn. We ignore this effect in this calculation, because the gas removed will have a composition very similar to that in the open pores. Thus, the removal has virtually no effect on the relative abundances of the gas species (although it certainly has an effect on the absolute inventory). This follows because gas-phase diffusion in the interconnected open pores is rapid and will homogenize the open pore inventory, quickly eliminating any compositional differences within a given isolated but internally-connected layer.
The number of moles of air molecules in the open pores can be represented by C. The ratio γ /C is then a rough measure of how much the air in the open pores has evolved due to permeation. For example, to explain an O 2 /N 2 enrichment of 7 ‰ in the lock-in zone would require a value of γ /C = 10 18 Pa m −1 m 2 s mol −1 .
The initial number of moles of each gas is the initial mole fraction times C: The physical fractionation factor that we seek is the ratio of permeabilities P 18 /P 16 . Substituting Eqs. (B2-B7) into (B1) and rearranging, and taking advantage of the fact that R STD equals the ratio of the mole fractions in both the bubbles and in the initial open pore air (e.g. x 18 /x 16 ), we get: P 18 P 16 = P 16 γ + C P N γ + C m 1 − P N P 16 + 1 = [δO 2 /N 2 + 1]m 1 − P N P 16 + 1 where we have rewritten the left-most term in square brackets as [δO 2 /N 2 + 1]. This term has a value of 1.007 for an O 2 /N 2 enrichment of 7 ‰, so it may be neglected with an error much smaller than our slope fit error. This fact also indicates that our use of a linear fit, rather than a curve, is an acceptable approximation. This leaves the final equation for the permeation fractionation factor α p : α p ≡ P 18 P 16 m 1 − P N P 16 + 1 (B9) From Ikeda-Fukazawa et al. (2005) we have P N /P 16 = 0.356, giving the fractionation factor α p as 0.9942 ± 0.0005, or a 5.8±0.5 ‰ fractionation. These uncertainties do not include possible errors in the permeation coefficients. We have no information on the error on these, but as an example, a 10% error in P N (with no error in P 16 ) would translate into an error in α p of 0.0003. Fortunately, because we employ ratios of permeation coefficients, systematic error in both values will tend to cancel.