Atmos. Chem. Phys., 12, 425492–77, 2012 www.atmos-chem-phys.net/12/4259/2012/ Atmospheric doi:10.5194/acp-12-4259-2012 Chemistry © Author(s) 2012. CC Attribution 3.0 License. and Physics Gas transport in firn: multiple-tracer characterisation and model intercomparison for NEEM, Northern Greenland C. Buizert1, P. Martinerie2, V. V. Petrenko3,*, J. P. Severinghau4s, C. M. Trudinger 5, E. Witrant 6, J. L. Rosen7, A. J. Orsi4, M. Rubino1,5, D. M. Etheridge1,5, L. P. Steele5, C. Hogan8, J. C. Laube8, W. T. Sturges8, V. A. Levchenko9, A. M. Smith9, I. Levin10, T. J. Conway11, E. J. Dlugokencky11, P. M. Lang11, K. Kawamura12, T. M. Jenk1, J. W. C. White3, T. Sowers13, J. Schwander14, and T. Blunier1 1Centre for Ice and Climate, Niels Bohr Institute, University of Copenhagen, Juliane Maries vej 30, 2100 Copenhagen Ø, Denmark 2Laboratoire de Glaciologie eteǴophysique de l’Environnement, CNRS, UniveérsJitoseph Fourier-Grenoble, BP 96, 38 402 Saint Martin d’H̀eres, France 3Institute of Arctic and Alpine Research, University of Colorado, Boulder, CO 80309, USA 4Scripps Institution of Oceanography, Univ. of California, San Diego, La Jolla, CA 92093, USA 5Centre for Australian Weather and Climate Research/ CSIRO Marine and Atmospheric Research, Aspendale, Victoria, Australia 6Grenoble Image Parole Signal Automatique, Unive´rJsoitseph Fourier/CNRS, Grenoble, France 7Department of Geosciences, Oregon State University, Corvallis, OR 97331-5506, USA 8School of Environmental Sciences, University of East Anglia, Norwich, NR4 7TJ, UK 9Australian Nuclear Science and Technology Organisation, Locked Bag 2001, Kirrawee DC NSW 2232, Australia 10Institut für Umweltphysik, University of Heidelberg, INF 229, 69120 Heidelberg, Germany 11NOAA Earth System Research Laboratory, Boulder, Colorado, USA 12National Institute of Polar Research, 10-3 Midorichou, Tachikawa, Tokyo 190-8518, Japan 13The Earth and Environmental Systems Institute, Penn State University, 317B EESB Building, University Park, PA 16802, USA 14Climate and Environmental Physics, Physics Institute, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland * current address: Earth and Environmental Sciences, University of Rochester, Rochester, NY 14627, USA Correspondence toC: . Buizert (christo@nbi.ku.dk) Received: 2 May 2011 – Published in Atmos. Chem. Phys. Discuss.: 26 May 2011 Revised: 10 April 2012 – Accepted: 23 April 2012 – Published: 14 May 2012 Abstract. Air was sampled from the porous firn layer at the single-tracer tuning. Six firn air transport models are tuned NEEM site in Northern Greenland. We use an ensemble ofto the NEEM site; all models successfully reproduce the data ten reference tracers of known atmospheric history to char-within a 1σ Gaussian distribution. A comparison between acterise the transport properties of the site. By analysing un-two replicate boreholes drilled 64 m apart shows differences certainties in both data and the reference gas atmospheriicn measured mixing ratio profiles that exceed the experimen- histories, we can objectively assign weights to each of thetal error. We find evidence that diffusivity does not vanish gases used for the depth-diffusivity reconstruction. We de-completely in the lock-in zone, as is commonly assumed. fine an objective root mean square criterion that is minimisedThe ice age- gas age differenc1ea(ge) at the firn-ice tran- in the model tuning procedure. Each tracer constrains the firnsition is calculated to be 18+239 yr. We further present the− profile differently through its unique atmospheric history and first intercomparison study of firn air models, where we intro- free air diffusivity, making our multiple-tracer characterisa- duce diagnostic scenarios designed to probe specific aspects tion method a clear improvement over the commonly usedof the model physics. Our results show that there are major Published by Copernicus Publications on behalf of the European Geosciences Union. 4260 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM differences in the way the models handle advective transportw. hole firn (Fabre et a,l.2000). Consequently, the diffusiv- Furthermore, diffusive fractionation of isotopes in the firn is ity profile with depth needs to be reconstructed for each firn poorly constrained by the models, which has consequenceasir site independently through an inverse methRodom( me- for attempts to reconstruct the isotopic composition of tracelaere et al,. 1997; Trudinger et al,. 2002; Sugawara et a,l. gases back in time using firn air and ice core records. 2003). The procedure consists of forcing a firn air transport model with the atmospheric history of a selected reference gas, often CO2, and subsequently optimising the fit to mea- 1 Introduction sured mixing ratios in the firn by adjusting (“tuning”) the effective diffusivity profile. The compacted snow (firn) found in the accumulation zone With few exceptions, the firn air modeling studies found of the major ice sheets acts as a unique archive of old air,in literature tune their effective diffusivity profile to a single preserving a continuous record of atmospheric compositiontracer. For the NEEM firn air site in Northern Greenland we from the present up to a century back in timBeat(tle et al,. have chosen to tune the effective diffusivity to an ensemble of 1996). Sampling of this archive has allowed for reconstruc- ten tracers, namely C2O, CH4, SF6, CFC-11, CFC-12, CFC- tion of the recent atmospheric history of many trace gas113, HFC-134a, CHCCl 14 153 3, CO2 andδ N2. The studies by species (e.gB. utler et al,. 1999; Sturrock et al,. 2002; Ay- Trudinger et al.(1997, 2002, 2004) also use a wide variety din et al., 2004; Montzka et al,.2004; Assonov et a,l.2007; of tracers, including halocarbons, greenhouse gases and ra- Martinerie et al,.2009) and their isotopologues (e.Fgr.ancey diocarbon 1( 14CO2), to characterise firn air transport. The et al., 1999; Ferretti et al,. 2005; Bernard et a,l.2006). Be- main difference with these works is that we show how multi- cause of its temporal range it naturally bridges the age gapple tracers can be combined in a methodical tuning process. between direct atmospheric observations and the ice coreBy systematically analysing all the uncertainties, both in the record (Etheridge et a,l.1998). Firn air analysis has some data and in the atmospheric histories of the references gases, significant advantages over the ice core technique. First, firnwe assign a unique weight to each data point in the tuning air can be sampled directly using a pumping linSech(wan- procedure. Using these uncertainty estimates we define an der et al,. 1993), making an ice extraction step unneces- objective root mean square deviation (RMSD) criterion that sary. Second, large sample sizes can be obtained making thiis minimised in the tuning. We introduce tracers that have method very suited for studying, e.g. recent changes in thenot previously been used in firn air studies, and provide at- isotopic composition of trace gaseSsu(gawara et a,l.2003; mospheric reconstructions for these species with realistic un- Röckmann et a,l.2003; Ishijima et al,. 2007). Third, because certainty estimates. bubble occlusion introduces additional smoothing, firn air Our approach has several advantages. A central difficulty records can achieve higher temporal resolution than ice coreisn reconstructing the diffusivity profile is that the problem is (Trudinger et al,.2004). under-determined, with (infinitely) many solutions optimis- Another important motivation for studying firn air is the ing the fit (Rommelaere et a,l1. 997). By adding more tracers, interpretation of ice core records. All air found in glacial ice each with a unique atmospheric history, the diffusivity profile has first traversed the firn, and has thus been affected by itiss constrained more strongly. The final reconstructed profile transport properties and the bubble close-off process. Somies a trade-off between the demands of the different tracers. of the most commonly described artifacts of firn air transportSecond, many trace gases of interest, such as halocarbons, are: gravitational fractionationC(raig et al,. 1988; Schwan- have a free air diffusivity that differs from C2Oby up to a fac- der, 1989), thermal fractionation in the presence of tempera- tor of 2. It is not a priori clear whether a profile tuned to 2CO ture gradientsS(everinghaus et a, l2. 001), diffusive smooth- alone provides a good solution for these gases. The tracers in ing of rapid atmospheric variationsSp(ahni et a,l. 2003), this study have a wide range of free air diffusivities, where molecular size fractionation during bubble close-oHffu(ber the fastest tracer (C4H) has a free air diffusivity that is three et al., 2006; Severinghaus and Bat,tle2006), diffusive iso- times that of the slowest one (CFC-113). Third, when using topic fractionation T( rudinger et al,.1997) and1age, the fi- ten tracers the available dataset is much larger, with several nite age difference between gas bubbles and their surroundd-ata points at each sampling depth. This makes the final re- ing ice (Schwander et a,l.1997). In certain cases the pecu- sult more robust, i.e. less susceptible to effects of outliers liarities of firn air transport actually give rise to new proxies, and analytical biases. Finally, our method is less sensitive to e.g. for temperatureS(everinghaus and Bro,o1k999; Landais errors in the reconstructed atmospheric history of individual et al., 2004; Dreyfus et al,. 2010) and local summer insola- reference gases. Our analysis shows that uncertainties in the tion (Kawamura et a,l.2007). atmospheric reconstruction are the largest source of potential Diffusion is the dominant mechanism by which variations error. in atmospheric composition are transferred into the firn. The Apart from presenting a new methodology for characteris- effective diffusivity decreases with depth as the pore spaceing firn air sites, this work will also serve as a reference for decreases; unfortunately, diffusivities measured on small firnother studies using NEEM firn air. The modeling in a number samples do not adequately describe the behaviour of the Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4261 of forthcoming publications in this ACP special issue will be a 0.9 based on the diffusivity reconstructions presented here. 0.8 We further present a model intercomparison between six state of the art one-dimensional firn air transport models. All 0.7 ×10−2 the models are tuned to the same dataset, using the same 0.6 s 8 physical firn parameters, such as porosity, free-air diffusion op coefficients, etc. To diagnose model performance we intro- 0.5 6 duce four synthetic scenarios that are designed to probe spe- 0.4 4 cific aspects of the model physics. s 2 For many sections in this work more information is avail- cl 0 able in the Supplement. Because of the vast amount of mate-b rial we have chosen to structure it the same as this work for 0.3 CZ easy referencing. Sections marked with an aste∗ri)sikn(this work have a corresponding section in the Supplement where 0.2 additional information can be found. The Supplement also includes all the firn air data and atmospheric reconstructions used in this study. 0.1 2 Methods 0 DZ LIZ 0 10 20 30 40 50 60 70 80 2.1 NEEM 2008 firn air campaign∗ Depth (m) The firn air used in this study was sampled during 14–30Fig. 1. (a)Firn densityρ (left axis), and the open and closed poros- July 2008 from two boreholes near the NEEM deep ice coreity sop andscl, respectively (right axis). Black dots show firn den- drilling site, Northern Greenland (77.◦4N5 51.06◦ W). The sity measurements in 0.55 m segments, the red line is the fit used in sampling site was 1.5 km outside of the main camp and cho-this study. Density data were kindly provided by S. J. Johnsen (per- sen to avoid contamination by going upwind of the prevailing sonal communication, 2009(b)) Gravitational enrichment as shown15 wind direction. The two boreholes, S2 and S3, were sepa-by δ N2, corrected for thermal effectsSe( veringhaus et a, l2.001). The blue line shows the barometric slo1pMe g(z−4.5)/RT , where rated by 64 m. Firn air was sampled by drilling through the we use1M = 1×10−3 kg mol−1. Convective zone (CZ), diffu- firn layers to the desired depth, and lowering the firn air sam-sive zone (DZ) and lock-in zone (LIZ) are indicated by changes in pling device into the borehole. The device consists of a purgeshading. and sample line, running through an inflatable bladder to seal off the firn air from the overlying atmosphere during sam- pling. S2 was sampled using the firn air system of the Uni-we use a combination of quadratic and exponential functions versity of Bern S( chwander et a,l.1993); S3 with the US firn to fit the data, keeping the derivative continuous over the tran- air system B( attle et al,. 1996). For this reason the boreholes sitions between the stages. The fit is generated o1nz a= are referred to as the “EU” and “US” holes, respectively. 0.2 m grid, which is the spatial resolution used in this study. The CO2 mixing ratio was monitored on the purge line with Following the temperature relationship givenSbcyhwander LICOR (US site) and MAIHAK (EU site) CO2 analysers to et al.(1997), a solid ice densityρice= 0.9206 g cm−3 is used. assess whether the firn air being pumped had been purged of We use an accumulation rate oAf= 0.216 m yr−1 ice modern air prior to sampling. equivalent, as derived from the 2007 NEEM S1 shallow core (D. Dahl-Jensen, personal communication, 2010). Since the 2.2 Physical characterisation of NEEM firn air site∗ model runs cover the period 1800–2008 CE, we used the av- erage rate over the same period as a best estimate. This long- The annual mean site temperature ◦−i2s8.9 C, as obtained term average is slightly lower than our best estimate for the from borehole thermometry on the EU hole; an isothermalcurrent day accumulation of 0.227 m−y1r. firn column is used for this modeling exercise. Atmospheric The total porositys = 1− ρ/ρice is divided into open and pressure is 745 hPa. We use the assumption commonly madcelosed porositys(op, scl) using the parameterisation Gofou- in firn air modeling that the firn density profile is in steady jon et al.(2003): state, and that the site has been climatically stable over the study period. For the firn density profi ( )ρle(z) we use an em- −7.6s pirical fit to the NEEM main ice core density measurementsscl = 0.37 s (1) sco averaged over 0.55 m segments (F1iga.). For each of the three stages of the densification process (grain boundary slidw- heresco is the mean close-off porosity. The deepest point ing, sintering and bubble compressioAnr;naud et al,. 2000) at which firn air could be sampled successfully wza=s www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 δ15N (‰) 2 ρ (g cm −3) porosity close−off depth lock−in depth 4262 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 77.75 m. The depth of total pore closurseop(= 0) was de- 2.4 Atmospheric histories of reference tracer∗s termined experimentally in the field by drilling the EU hole to a depth ofz= 83 m, and subsequently raising the blad- A best-estimate atmospheric history was reconstructed for der in 0.5 m steps. The deepest point at which air could each of the reference tracers used in the tuning. Additionally, be pulled wasz= 79 m, though the flow was insufficient we reconstructed δa13CO2 history which was used to con- for sampling. To be consistent with this observation we usevert the114CO2 history to a ppm scaleS(tuiver and Polac,h sco= 9.708×10−2 in Eq. (1), which leads to total pore clo- 1977). All the atmospheric histories start in the year 1800 sure (sop= 0) at depthz= 78.8 m, as shown in Fig1.a. and have monthly resolution. We use 2008.54 (mid July) as Traditionally, the firn column is divided into three zones, the decimal sampling date in the models. Four different types based on the gravitational enrichment with depSthow(ers of data were used in the reconstructions: et al., 1992). For NEEM theδ15N2 together with the zonal division is shown in Fig1. b. Using the barometric line fitting 1. Direct atmospheric measurements from sampling net- method K( awamura et a,l.2006), we obtain a convective zone works or archived air. Northern hemisphere high lati- (CZ) thickness of 4.5 m. The zone below is called the diffu- tude stations (i.e. Alert, Summit and Barrow) are used sive zone (DZ), as diffusion dominates the transport at these whenever available. When using data from other sta- depths. The lock-in depth is defined as the depth at which tions a correction should be made to account for the gravitational enrichment stops, which happens at63 m. latitudinal gradient; this is done whenever the gradi-z= Below we find the lock-in zone (LIZ), where advection with ent could be reliably determined. We used station data the ice matrix dominates the transport. This is also the region from the NOAA-ESRL network for CO2, CH4, SF6 and where the majority of the bubble occlusion occurs, as can be HFC-134a C( onway et al,. 1994; Dlugokencky et a,l. seen from the curve. 1994; Geller et al,. 1997; Montzka et al,. 1996); fromscl ALE/GAGE/AGAGE for CFC-11, CFC-12, CFC-113 2.3 Gas measurement∗s and CH3CCl3 (Prinn et al,. 2000, 2005; Cunnold et al,. 1997; Fraser et a,l. 1996); from CSIRO for δ13CO2 A total of 345 samples were taken from the boreholes, with (Francey et a,l.2003). All ALE/GAGE/AGAGE data atmospheric samples taken at three occasions during the were converted to the latest NOAA calibration scales sampling period. Different flask types were used (SilcoCan, that are also used in the firn air analysis; δa1l3lCO2 stainless steel flasks and glass flasks); for the tracers used in data are kept on the CSIRO calibration scale. Further- this study no effect of flask type could be found in data com- more, we used the SIO Mauna Loa record for 2CO parison. We use firn air data from six different laboratories; between 1958–1985 CEKe(eling et al,. 2009), and an overview is given in Tabl1e. the Cape Grim air archive foδr13CO2 between 1978– For CH4 and SF6, data from IUP were corrected for known 1991 (Francey et a,l. 1999). The 14CO2 history be- calibration offsets to place all data on the NOAA scales used tween 1963–1993 is based on Fruholmen, Norway, for in the atmospheric reconstructions; corrections were always its proximity to GreenlandN(ydal and L̈ovseth, 1996); smaller than 3 ppb and 0.06 ppt, respectively. After the cal- data from other stations were used to complete time ibration correction no systematic offsets were observed be- coverage M( anning and Melhuis,h1994; Levin and tween the different laboratories. As discussed in S2e.c6t,. Kromer, 2004; Levin et al., 2008). there are significant offsets between the two boreholes, and data from the holes are not combined but treated separately. 2. High resolution firn air/ice core measurements from the The114CO2 data and atmospheric reconstruction were both high accumulation Law Dome sites, Antarctica. The re- converted to a (mass conserving) ppm scale to a1ll4oCwO constructions of CO2 before 1958, and C4Hbefore 19832 to be modeled like a regular tracer. Wherever multiple data are based onEtheridge et al(.1996, 1998) andMacFar- 13 points were available for a certain depth, they were averaged ling Meure et al.(2006); δ CO2 before 1976 is based in the following order: first replicate measurements from onFrancey et al(.1999). one laboratory were averaged, then the resulting values from the different labs were averaged. In this way one composite 3. Dendrochronologically-dated tree-ring measurements dataset was created that contains 204 data points for the 10 of radiocarbon. The reconstruction of atmospheric14 tracers on the EU hole, and 77 data points for the 4 tracers CO2 before 1955 is based oRneimer et al.(2004). on the US hole. Theδ15N2 profile is the same for both holes, 4. Emission-based estimates using a 2-D atmospheric and a combination of EU and US depths are used in the final transport model that includes latitudinal source and sink dataset. This gives a final dataset of 260 data points. distribution (Martinerie et al. (2009) and references therein). This has been used to complete time coverage for SF6 and halocarbons before the onset of direct at- mospheric measurements. Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4263 Table 1.Overview of firn air data used in this study. Acronyms rep- resent: the School of Environmental Sciences at the University of a EU borehole East Anglia (UEA;Laube et al,.2010), NOAA Earth System Re- 370 US borehole search Laboratory, Boulder CO (NOAAC;onway et al,.1994; Dlu- 360 gokencky et a,l.1994), the Institut f̈ur Umweltphysik at the Uni- versity of Heidelberg (IUP;Levin et al., 2010, 2011), the Com- 350 monwealth Scientific and Industrial Research Organisation, Ma- rine and Atmospheric Research (CSIRFOra; ncey et a,l.2003), the 340 Australian Nuclear Science and Technology Organisation (ANSTO; Smith et al,. 1999), and Scripps Institution of Oceanography at the 330 University of California, San Diego in collaboration with the Na- tional Institute of Polar Research, Japan (SIO/NIPSRev;eringhaus et al., 2003). b 1800 Tracer EU borehole US borehole 1700 CO2 NOAA, CSIRO, IUP NOAA, IUP CH4 NOAA, CSIRO NOAA, IUP 1600 SF6 NOAA, IUP, UEA NOAA, IUP CFC-11 UEA – 1500 CFC-12 UEA – CFC-113 UEA – HFC-134a UEA – 1400 DZ LIZ CH3CCl3 UEA – 62 64 66 68 70 114CO2 ANSTO – 15 Depth (m)δ N2 SIO/NIPR SIO/NIPR δ86Kr ∗ SIO/NIPR SIO/NIPR Fig. 2.Comparison o(fa) CO2 and(b) CH4 data in the lock-in zone ∗ Used for gravitational correction only of the two boreholes. Data from three laboratories are averaged; the 1σ uncertainty bar is a combination of analytical and sampling uncertainties as specified in Se2c.t7.. Curves are modeled profiles 2.5 Gravitational correction using the CIC model, see Se3ct..3. The gravitational fractionation of gases in the firn column is well established both theoretically and experimentaCllrya(ig 2.6 Differences between the EU and US boreholes et al., 1988; Schwande,r1989; Sowers et a,l.1992). Before the modeling, all data were corrected for the effect of gravity:Though the boreholes are separated by a mere 64 m, within the lock-in zone we find differences in the mixing ra- [X]meas(z) tio profiles for CO2, CH4 (Fig. 2) and SF6 (not shown) [X]corr(z)= (2) 1M(δ (z)/1000 1) that exceed the estimated uncertainty of the combinedgrav + sampling-measurement procedure (indicated with errorbar, where[X]corr ([X]meas) is the gravity corrected (measured) see Sect2. .7for details). Note that we cannot make the same mixing ratio of gas species X1,M =MX −Mair is the mo- comparison for the other gas species since only data from the lar mass difference between gas X and air, δagnradv(z) is the EU hole are available. The estimated uncertainty on the depth gravitational fractionation per unit mass difference at depthmeasurements is 5 cm, whereas at certain depths an error of z (Supplement). Theδgrav(z) values are based on measure- 50–60 cm is required to explain the observed borehole off- ments ofδ86Kr (86Kr/82kr) with depth, and corrected for the set. So although depth measurement errors can certainly con- effect of thermal fractionationS(everinghaus et a, l2. 001). tribute to the offset, they cannot explain the full magnitude of The rationale for using Kr rather than2,Nis that its free-air it. We attribute the differences to lateral inhomogeneity in the diffusivity is closer to that of most of the tracers we use, sofirn stratigraphy, possibly originating from surface wind fea- it should represent the disequilibrium effects on gravitationaltures that have been preserved in the densification process. fractionation more accurately. The models are run with ei-Because of these differences we have chosen to model both ther the molecular weight of all gases set equal to that ofholes separately. air (M =Mair) or gravity set to zerog(= 0). This empiri- Furthermore, for S6Fwe observe a∼ 0.25 ppt offset in the cal gravitational correction ensures that the effect of gravity5–50 m depth range, contrary to the 2COand CH4 profiles is included correctly, and could potentially also be used tothat agree well between the holes at these depths. This ex- correct for mass-dependent sampling artifaSctesv(eringhaus cludes differences in age distribution of the air as the origin. and Battle, 2006). Alternative explanations we considered, such as incomplete www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 CH (ppb) CO (ppm) 4 2 4264 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM Table 2.Relative contributions to the total uncertainty, averaged over the firn column (EU borehole). Tracer Analytical Reconstruction Contamination Sampling Other CO 0.02 0.47 0.00 0.17 0.1a2 5, 0.19 b CH4 0.08 0.85 0.01 0.05 – SF6 0.07 0.54 0.04 0.08 0.2 c7 CFC–11 0.11 0.81 0.04 0.04 – CFC–12 0.35 0.48 0.06 0.11 – CFC–113 0.55 0.38 0.04 0.03 – HFC–134a 0.29 0.63 0.04 0.04 – CH3CCl3 0.10 0.73 0.00 0.17 – 114CO2 0.14 0.80 0.00 0.00 0.0 a6 δ15N2 0.91 – 0.05 0.04 – a In-situ artifactsb Undersampling of seasonal cyccleEU-US borehole S6Foffset flask flushing, sample contamination, procedural blanks and3 Modeling firn air transport at NEEM bladder outgassing, could all be excluded. Since we found no objective reason to reject data from either hole, we accoun3t .1 Tuning of the diffusivity profile ∗ for the discrepancy by assigning an additional errorbar to the SF data from both holes in this depth range (see below). How the diffusive transport changes with depth needs to be6 reconstructed through an inverse method for each firn air site 2.7 Full uncertainty estimation∗ independently. The procedure consists of forcing the trans- port models with the atmospheric history of one or several Having accurate uncertainty estimates is essential for ourselected reference tracer(s), and subsequently optimising the multiple-tracer methodology. A unique uncertainty estimatefit to the measured mixing ratios in the firn by adjusting the has been assigned to each of the 260 individual data pointesffective diffusivity values at each deptRho(mmelaere et a,l. used in this study based on the following seven potential1997; Trudinger et al,.2002; Sugawara et a,l.2003). sources of uncertainty: (1) analytical precision as specified Here we will use diffusion as a generic term for mass trans- by the laboratories; (2) uncertainty in atmospheric recon-fer processes that are driven by a concentration gradient, and structions; (3) contamination with modern air in the deep- well described by Fick’s law. The diffusivity profile we re- est firn samples; the estimates are based on HFC-1346a, SFconstruct is composed of two parts: and CFCs, which should be absent in the deepest samM- olecular diffusionis a microscopic process originating in ples; (4) sampling effects estimated from inter-laboratory andthe thermal motion of the molecules constituting the firn air. inter-borehole offsets; (5) possibility of in-situ C2Oartifacts The effective molecular diffusion coefficient of gases in the (e.g.Tschumi and Stauffe, r2000) and CO2 enrichment due open porosity decreases with depth as the diffusive path be- to close-off fractionationH(uber et al,. 2006; Severinghaus comes increasingly tortuous due to the densification process and Battle, 2006); (6) undersampling of seasonal cycle in (Schwander et a,l.1988). The diffusion coefficient of gaXs the monthly atmospheric reconstruction (2CoOnly); and (7) at depthz is given by large unexplained EU-US borehole difference in the diffusive zone (SF only). 0 D06 DX CO2 The errors are assumed to be independent of each otheDr.X(z)= sop = sopγX (3)τ(z) τ (z) The uncertainties in the atmospheric reconstructions are con- verted from a time scale to a depth scale by running themwhereD0X is the diffusion coefficient in free airτ,(z) is the through the CIC firn air model (Sec3t.2). Details are given tortuosity at depthz, andγX =DX/DCO2 is the ratio of dif- in the Supplement. The relative contribution of the sevenfusion coefficientsT(rudinger et al,.1997). It is clear from sources, averaged over the firn column, is given in Ta2b. le Eq. (3) that the molecular diffusivity profiles for the different For most tracers the atmospheric reconstruction is the largesgtas species scale linearly with each other. γTXheused in this contributor to the total uncertainty, followed by the analytical study are based on measurementsMbaytsunaga et a(l.1993, precision. 1998, 2002a,b), or a theoretical formulaC(hen and Othme, r 1962) for gas species where no experimental data are avail- able (Supplement). Eddy diffusionrefers to mass transfer caused by air flow pat- terns in the open porosity which the models cannot resolve directly, and are instead parameterised through the inclusion Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4265 Table 3.Overview of firn air transport models. See text for explanation. Model Coordinates Advection Convection Time stepping Tuning CIC static flux exponential Crank-Nic. automated CSIRO moving coordinates exponential Crank-Nic. automated INSTAAR static flux exponential explicit manual LGGE-GIPSA static flux tuned implicit automated OSU static boxes exponential Crank-Nic. manual SIO static boxes exponential explicit automated of a diffusion coefficientDeddy(z). Unlike molecular diffu- Research Organisation (CSIROT,rudinger et al,.1997), In- sion,Deddy is equal for all gases as the process is macro-stitute of Arctic and Alpine Research (INSTAAR), Labora- scopic in origin. The first contribution tDo eddy is convective toire de Glaciologie et Ǵeophysique de l’Environnement and mixing in the top few meters due to seasonal temperatureGrenoble Image Parole Signal Automatique (LGGE-GIPSA, gradients and wind pumpingC(olbeck, 1989; Severinghaus Witrant et al,. 2011), Oregon State University (OSU) and the et al., 2001; Kawamura et a,l.2006). Some of the models Scripps Institution of Oceanography (SIOSe, veringhaus and also include aDeddy term in the deep firn to represent disper- Battle, 2006). sive mass transfer. The classical example of dispersive mix- Table3 summarises a selection of relevant model charac- ing is Taylor dispersion, where (viscous) shear flow in a cir- teristics. The second column indicates the kind of coordinate cular tube enhances the effective diffusivity of a soluAtreis(, system used by each model. The CSIRO model is unique in 1956). Viscous flow can be induced in the deep firn by low using a coordinate system that moves downwards with de- frequency pressure fluctuations at the surfaScceh(wande,r scending ice layers (Lagrangian coordinates), giving a down- 1989), as well as air expulsion due to pore compactiRoonm( - ward ice velocitywice= 0 relative to the reference frame. melaere et a,l.1997). In soils, viscous flow induced by atmo- The other models use a static (Eulerian) coordinate system spheric pressure fluctuations has been shown to be an imw- here the surface stays zat= 0 and the ice layers move portant transport mechanismMa( ssmann and Farri,e1r992). down at a finite velocitywice= Aρ/ρice. The choice of coor- Pore compaction causes a non-uniform velocity distributiondinate system has consequences for the way advection of air between air expelled upwards, and air remaining behind inwith the ice matrix is handled (third column). In the CSIRO the pores. This spread in velocities results in additional (dis-model, the moving coordinate system automatically leads to persive) mixing. advective transport without the need to include it explicitly. The extensive uncertainty analysis we introduced inThe models expressed in Eulerian coordinates use either an Sect.2.7 assesses how reliable each data point is, and cona-dvective flux as described bRyommelaere et a(l.1997), or sequently how much weight should be assigned to it in theuse boxes of air that are shifted downwards at regular time tuning procedure. This allows us to combine multiple tracersintervals (Schwander et a,l.1993). Convection (fourth col- in the tuning, and the final reconstructed diffusivity profile umn) is handled in two different ways. The LGGE-GIPSA is a trade-off between the constraints placed by the differentmodel uses aDeddy term that is tuned by the RMSD min- tracers. All models in this study tuned their effective diffu- imisation algorithm, in combination with zero gravitational sivity profi(les to minimise th)e root mean square deviation fractionation forz < 4 m. The other models use the parame-∑ terisation ofKawamura et al.(2006) where aD1 eddy term is1 N 2(m d )2− included that falls off exponentially with depth. The fifth col- RMSD i i= (4) N 2 umn of Table3 indicates whether the time stepping used in i=1 σi solving the diffusion equation was either explicit (Euler for- where thed are the data for a given borehole (all tracers), ward), implicit (Euler backward) or mixed implicit-expliciti m are the linearly interpolated modeled values at the same(Crank-Nicolson method). Finally, some of the models werei depths, andσi give the total uncertainties. The indei xruns tuned through manual adjustments of the diffusivity profile, over all data points, wherNe = 204 for the EU hole, and while others used an automated control method (Reo.gm. - N = 77 for the US hole. melaere et a,l.1997; Trudinger et al,.2002). 3.2 Model description∗ 3.3 Fit of modeled profiles to the data∗ Six 1-D firn air transport models are tuned to the NEEM 2008 firn air data in the way outlined above. The modelsThe modeled profiles for all 10 tracers from the EU borehole (in alphabetical order) were developed at the Centre for Iceare shown in Fig.3, profiles from the US borehole can be and Climate (CIC), Commonwealth Scientific and Industrial found in the Supplement. www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 4266 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70 a 380 DZ LIZ f 250 DZ LIZ 200 360 CIC 150 CSIRO 340 INSTAAR 100 LGGE−GIPSA 320 OSU 50 SIO 300 0 b g 1800 500 400 1600 300 1400 200 100 1200 0 c h 6 80 60 4 40 2 20 0 0 d i 100 50 80 40 60 30 40 20 20 10 0 0 e 550 j 0.3 0.25 500 0.2 0.15 450 0.1 400 DZ LIZ 0.05 DZ LIZ 0 0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70 Depth (m) Depth (m) Fig. 3. (a–j)Modeled profiles for all 10 tracers from the EU borehole. With the exceptio(jn) odfata are gravity corrected and the models are run with gravity turned off. Errorbars correspond to fuσll 1uncertainty as defined in Sec2t..7. Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ 14CO (10−12 ppm) CH CCl (ppt) SF (ppt) CH (ppb) CO (ppm) 2 3 3 6 4 2 δ15N (‰) 2 HFC−134a (ppt) CFC−113 (ppt) CFC−12 (ppt) CFC−11 (ppt) C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4267 −4 −2 0 2 4 −4 −2 0 2 4 a 30 bCIC 30 CSIRO RMSD= 0.73 RMSD= 0.92 20 20 10 10 0 0 c 30 dINSTAAR 30 LGGE−GIPSA RMSD= 0.78 RMSD= 0.74 20 20 10 10 0 0 e 30 fOSU 30 SIO RMSD= 0.92 RMSD= 0.80 20 20 10 10 0 0 −4 −2 0 2 4 −4 −2 0 2 4 (Model−data) / σ (Model−data) / σ Fig. 4. (a–f)Histogram of(mi − di)/σi for the firn air transport models in this study using the EU borehole data. The black curve gives a Gaussian distribution of widtσh = 1, normalised to have equal surface to the histogram. The RMSD is calculated usin4g).Eq. ( For CO2 (Fig. 3a) we find a pronounced mismatch for these peaks that the divergence between the models is most z > 70 m which is reproduced consistently by all models easily visible. It must be noted, however, that equally large and in both boreholes. The atmospheric reconstruction ismodel differences are found for other tracers as well (e.g. based on the Law Dome recordEth( eridge et a,l.1996), CO2 and HFC-134a). Furthermore, models that fit the peak and consequently the largest source of uncertainty is the inh- eight exactly do not necessarily provide the best overall fit terhemispheric CO2 gradient. To explain the observed mis- to the data. match with an error in the atmospheric history would re- To quantitatively assess how well the modeled profiles fit quire a 6 ppm interhemispheric gradient in the 1950s, whichthe data, we make a histogram (omf i − di)/σi where the can be ruled outK(eeling et al,. 2010). Also contamination index i goes over all 204 data points of the EU borehole. with modern air can be ruled out for these flasks, basedTheσi are the unique uncertainties we assigned to each data on measurements of halocarbons. We cannot, however, exp-oint. This is shown in Fig4. together with a Gaussian dis- clude in-situ production of C2Ofrom, e.g. organic material tribution of width σ = 1 and a surface area equal to that or (bi)carbonates found in Greenland ice (eT.gsc. humi and of the histogram. The figure furthermore shows the RMSD Stauffer, 2000). There might also be a small close-off frac- from the data as given by Eq4.)(. For all models we find a tionation of CO2, of order 1 ‰ or less, based on its effective distribution that is more narrow than the Gaussian distribu- molecular diameterH(uber et al,. 2006; Severinghaus and tion, meaning that within the assigned uncertainties all data Battle, 2006). from the borehole can be modeled consistently. The good Both CH 143CCl3 and CO2 (Fig. 3d–e) have a peak within agreement gives confidence in the correctness of our recon- the LIZ. The peak height and position, which are sensitive tostructed atmospheric histories, which will be of use in future the shape and magnitude of the diffusivity profile, provide anfirn air studies. For the US hole the RMSD ranges between important constraint in the tuning procedure. It is exactly at0.60–1.06. www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 Counts Counts Counts Counts Counts Counts 4268 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 0 10 20 30 40 50 60 70 a DZ LIZ 600 no LIZ diffusion1.5 with LIZ diffusion 550 1 500 0.5 450 0 b 400 DZ LIZ −5 10 40 50 60 70 −6 10 Depth (m) −7 CIC Fig. 6. Modeled 14CO2 profiles both with, and without diffusivity10 CSIRO in the LIZ. Calculations done with the CIC firn air model. −8 10 INSTAAR LGGE−GIPSA −9 OSU 10 4 Model intercomparison and discussion SIO c 1 4.1 Diffusivity profiles∗ 0.8 All the models are tuned separately to optimise the fit to the data. Figure5a shows the reconstructed molecular diffusivity 0.6 profiles for CO2. The spread between the solutions is caused by two factors. First, it reflects differences in model physics, 0.4 which are compensated for by adjustments to the diffusivity profile. Second, it relates to the degree in which the inverse 0.2 problem of diffusivity reconstruction is under-determined. By using 10 different tracers the problem is more strongly 0 constrained than in the case of using only2C, Obut the solu- 0 10 20 30 40 50 60 70 tion is not unique. Depth (m) A second important observation is that all models require a non-vanishing diffusivity on the order of 210−9 m2s−1× Fig. 5. (a)CO2 molecular diffusivity profile with depthDCO2(z) for within the LIZ (see Fig.5b). This contradicts the notion that the EU borehole(.b) Semi-log plot of the total CO2 diffusivity pro- file . (c) Plot of . diffusivity vanishes completely in the LIZ due to the pres-Dtotal(z)=DCO2(z)+Deddy Deddy(z)/Dtotal(z) Eddy diffusion near the surface corresponds to convective mix-ence of impermeable layers, which is frequently expressed ing; in the LIZ some of the models include dispersive mixing. The in the firn air literature (e.gB. attle et al,. 1996; Trudinger CSIRO model has zero eddy diffusivity in the LIZ; the line is not et al., 2002; Assonov et a,l.2005). Our findings confirm a re- visible as the LGGE-GIPSA model is plotted on top of it. cent study at Megadunes, Antarctica, that also reports a finite (eddy) diffusivity in the LIZ (Severinghaus et a, l2. 010). We believe this result to be robust, since it is reproduced by all Two models stand out as having a larger RMSD of 0.92six firn air models. In particular, the well reproduced height (CSIRO and OSU). For the CSIRO model this is caused byof the14CO2 bomb spike, which falls entirely within the LIZ, the absence of a back flow due to pore compression in theprovides strong evidence for finite LIZ diffusion at NEEM; transport description (Sec4t..4.4). For the OSU model we letting diffusivity go to zero leads to1a4CO2 peak that is too attribute the larger RMSD to the tuning procedure for the narrow and overshoots the measured value∼s b7yσ (Fig.6). molecular diffusivity, which has fewer degrees of freedom How the LIZ diffusivity is parameterised, however, varies than the procedures used by the other models (Supplement)s.trongly among the models as shown in F5icg.. The models Differences in RMSD give information on the perfor- use either molecular diffusion, dispersive mixing (eddy dif- mance of models, or model configurations, relative to eachfusion) or a mixture of both. To fit thδe15N2 data, most mod- other; the RMSD as defined by Eq4.)(should not be inter- els require some dispersive mixing in the LIZ, since molec- preted in an absolute sense. ular diffusion alone would lead to continued gravitational Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ D / D 2 −5 eddy total D (m /s) D (10 m 2/s) total CO2 14CO (× 10−12 ppm) 2 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4269 a 14 a CIC z = 63 m −5 10 DZ LIZ 12 CSIRO INSTAAR 10 LGGE−GIPSA −610 OSU 8 SIO −7 6 10 4 −8 10 EU borehole 2 US borehole −9 0 10 1980 1985 1990 1995 2000 2005 2010 0 10 20 30 40 50 60 70 Year Depth (m) b b 12 4 z = 78 m 10 z = 63 m 3 8 6 2 4 z = 76 m 1 2 0 0 1880 1900 1920 1940 1960 1980 1900 1920 1940 1960 1980 2000 Year Year Fig. 7. EU borehole modeled C2Oage distribution densities for Fig. 8. Comparison between the EU and US boreho(lae)s.Semi- (a) depthz= 63 m (lock-in depth) and(b) depthz= 78 m (bot- log plot of the total CO2 diffusivity profile Dtotal(z)=DCO2(z)+ tom of LIZ). On the horizontal axis are calendar years C.E.; dec-Deddy(z) averaged over all six firn model(sb. ) Age distributions imal sampling year is 2008.54 (i.e. mid July). Age distributions at depthsz= 63 m (lock-in depth) andz= 76 m (bottom of LIZ) are calculated by applying a surface forcing which is unity for averaged over five firn models. The CSIRO model is excluded in the 0.2≤ t < 0.4 yr, and zero elsewhere. age distribution comparison as the US age distributiozn=at76 m is unreliable (Supplement). fractionation. The LGGE-GIPSA model circumvents this problem by using a combination of Fick’s and Darcy’s trans- 4.2 Gas age distributions∗ port to describe almost stagnant transport in the lower zones (Witrant et al., 2011), which allows for using purely molec- Firn air does not have a single age, but rather a distribution ular diffusion in the LIZ. The SIO model uses the ratio of of ages S( chwander et a,l.1993). The modeled gas age dis- eddy to molecular diffusion in the LIZ as a tuning parame- tribution, also referred to as Green’s function, transfer func- ter in the RMSD minimisation, and finds that the fit is op- tion, response function or age spectrum, provides a complete timised for 27 % eddy diffusion (Supplement). Models that description of the model transport propertieRsom( melaere reproduce the observations equally well can have completelyet al., 1997). Figure7 compares CO2 age distribution densi- different parameterisations, so our analysis does not tell usties for the models at the lock-in deptzh=( 63 m) and at the which scheme is more likely to be correct. These model dif-bottom of the LIZ (z= 78 m, close to the deepest EU hole ferences do have important consequences for the modelinsgample atz= 77.75 m). The relevant characteristics of the of isotopic ratios, as is discussed in Se4c.4t..1. age distributions are summarised in Ta4b.leAs a measure of the relative spread in model results, we list thσes2tandard deviation divided by the meanµ(). www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 Density ( ×10−2 yr−1) Density ( ×10−2 yr−1) 2 −1 Density ( ×10−2 yr−1) D (m s )total 4270 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM The LIZ lies between the two depths depicted here. FromTable 4.Mean age, Median age, Full Width at Half Maximum and the closed porosity parameterisation (E1q),. we find that Spectral Width1( , Eq. (1) inTrudinger et al,.2002) at the lock-in ∼ 95 % of the air is trapped in this zone. The gas age dis-depth (z= 63 m) and at the bottom of the LIZz (= 78 m). All values tribution found in NEEM ice below the full close-off depth given in years. Theσ2 standard deviation divided by the meµan (sop= 0) will be intermediate to these two end members, gives a measure of the spread in the model results. with a tail of older air trapped already in the DZ. Interest- ingly, at the lock-in depthz(= 63 m), estimates of both the Model Mean Median FWHM 1 age and the distribution width show a spread of around 25 % z= 63 m between models that reproduce the data equally well. Clearly CIC 9.2 7.4 7.9 4.5 these properties are not as well determined as one would CSIRO 10.3 8.8 9.3 4.2 expect based on the similarity between the modeled curves INSTAAR 8.7 7.0 7.4 4.2 of Fig. 3 alone. On traversing the LIZ tzo= 78 m the rela- LGGE-GIPSA 11.1 8.1 9.1 6.8 tive spread in modeled ages is reduced; the absolute spread OSU 7.9 6.5 7.0 3.7 grows slightly to about 5 yr, though. The spread in distribu- SIO 9.0 7.2 7.6 4.5 tion widths remains large. 2σ/µ 0.25 0.22 0.23 0.47 For the US borehole the spread in the calculated mean ages z= 78 m at z= 63 m is even larger (41 %). This large spread is due to the fact that we have only four tracers on the US borehole, CIC 70.9 69.4 30.1 9.8 and therefore the mean ages are less strongly constrained by CSIRO 67.6 66.0 24.8 8.3 the data. INSTAAR 68.3 66.1 37.1 12.4 We can use the gas age distributions to calculate the LGGE-GIPSA 72.8 70.4 36.8 12.5 OSU 67.5 65.7 33.3 10.8 NEEM modern day1age, i.e. the difference between gas age SIO 69.5 67.7 32.2 10.5 and ice age below the firn-ice transition. From annual layer 2σ/µ 0.06 0.06 0.28 0.30 counting and matching of reference horizons to the NGRIP GICC05 timescaleR(asmussen et a, l2.006), the ice at depth z= 63 (78) m is estimated to be 190.6 (252±.51) yr old (S. O. Rasmussen, personal communication, 2011). To calcu- late the true1age we would require the gas age in the closed porosity. With the model results of Tab4lewe can only cal- the US hole has undergone more diffusive mixing, leading culate1ageop, i.e. the difference between ice age and gas ageto a broader distribution width. Note that these differences in in the open pores. By averaging results from the six firn airtransport properties appear to have only a minor impact on models we find a1ageop of 181±2 and 183±2 yr atz= 63 the1age calculated in Sec4t.2. and 78 m, respectively. This shows that the LIZ, ice and air Polar firn is a layered medium that exhibits density vari- are aging at roughly the same rate with depth. An estimatedations with depth caused by seasonal variations in local cli- 4–5 % of the air is already trapped above 63 m, which biasesmatic conditions and deposited snow density (eH.ögr.hold our air age estimate young. This air fraction can introduce anet al., 2011). Snow drift and redeposition at the surface lead error of at most 8 yr. Based on these considerations, our bestot lateral inhomogeneities in the firn stratigraphy. The bore- estimate of the true1age is 182+3 yr (181+39 9 yr) for the EU hole comparison suggests that firn air transport in the LIZ is− − (US) borehole. very sensitive to these lateral variations, whereas transport in the DZ is not. This is not unexpected since the formation of 4.3 Borehole comparison the LIZ has been linked to the degree of layerinLgan(dais et al., 2006). If the amplitude of density variations is suf- Here we compare the reconstructed diffusivity profiles andficiently large, the high-density layers will reach the close- gas age distributions for the EU and US boreholes. Fig8aure off density before the low-density layers do, thus creating shows the total diffusivity with depth for C2Oon a semi- sealing layers that impede vertical diffusion and prevent fur- log scale, where we have averaged over the solutions fromther gravitational fractionation of the gas mixtuMrea(rtinerie the different firn air models. We observe that the largest di-et al., 1992; Schwander et a,l1. 993; Sturrock et al,.2002). In vergence between the boreholes occurs in the LIZ. Also theSect.4.1 we presented evidence for finite vertical mixing in gas age distributions show divergence mainly in the LIZ. Fig- the LIZ despite the presence of such impeding layers. The ure8b shows the age distributions at two depths, where againamount of vertical mixing is expected to depend on the hori- the curves represent an average over model output from thzeontal extent of the sealing layers, as well as the position of different firn air models. At the lock-in depthz=( 63 m) the (micro)cracks and passageways. This could lead to a strong age distributions of both boreholes are still very similar. Af- lateral variation in LIZ diffusion as observed. Note that such ter traversing the LIZz(= 76 m, chosen to be near the final lateral processes cannot be represented adequately in one- sampling depth on the US hole), we find that the firn air in dimensional transport models. Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4271 a 0 10 20 30 40 50 60 70 b 0 10 20 30 40 50 60 700 010 −0.02 −0.04 −110 −0.06 −0.08 −210 CIC −0.1 CSIRO INSTAAR −3 −0.12 LGGE−GIPSA 10 −0.14 OSU SIO −4 10 c 0.3 d 800 0.25 0.2 600 0.15 400 0.1 0.05 200 0 DZ LIZ DZ LIZ0 0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70 Depth (m) Depth (m) Fig. 9. Model comparison using the four diagnostic scenarios and diffusivity tuned to the EU borehol(ea)dSatcae.nario I: diffusive frac- tionation for a hypothetical monotonic C2Oincrease.(b) Scenario II: attenuation of a 15 yr period sinusoidal2CfOorcing with depth.(c) Scenario III: gravitational enrichment for gXaswith D0 = 0.025D0CO . Data show gravitational enrichment 1o5fN2 corrected for the effectX 2 of thermal diffusion. See Supplement for deta(ilds). Scenario IV: mean age of gaYs when using advective transport onlDy 0( = 0, i.e. Y diffusion is absent). With the exception of S-III all scenarios were run with the effect of gravity turned off. 4.4 Synthetic diagnostic scenario∗s isotopologue has a lower diffusivity than th12eCO2 isotopo- logue, it is transported less efficiently into the firn, leading to To diagnose model properties further, we developed foura depletion inδ13CO2 with depth. Ice core and firn air records synthetic scenarios that probe specific aspects of the modenleed to be corrected for the effect of diffusive fractionation physics. We present the individual scenarios below, and dis-(DF). For example, the (site specific) corrections for the re- cuss the model differences they reveal. cent anthropogenic increase a∼r0e.1 ‰ forδ13CO2 at DE08, Law Dome (Francey et a,l.1999), and 1.2 ‰ forδ13CH4 at 4.4.1 Scenario I: diffusive fractionation WAIS-Divide (Mischler et al,. 2009). Figure9a shows the modeled isotopic signal with depth for For gas species whose atmospheric mixing ratios changethe different models. Results in the DZ are consistent, how- over time, the disequilibrium in the firn leads to an isotopic ever they diverge strongly once the LIZ is reached.δF1o3Cr fractionation with depth, even in the absence of changes in ato- f CO2, the observed range between the different model so- mospheric isotopic compositionTr(udinger et al,.1997). As lutions in the deepest firn is 0.05 ‰, which is around 40 % a scenario we use a monotonic 2CiOncrease o∼f 1 ppm yr−1 of the total DF signal. For isotopologues of C4 Hthe DF without seasonal cycle, chosen to resemble the current ani-s larger, because of the larger relative mass difference be- thropogenic increaseK(eeling et al,. 2010). The surface iso- tween the isotopologues. Based on the observed DF model topic ratio is kept fixed aδt 13CO2 = 0. Because the13CO2 www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 δ (‰) 13 grav δ CO (‰)2 Age of air (yrs) Amplitude |c| 4272 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM disagreement we estimate an uncertainty of 0.5 ‰ for both4.4.2 Scenario II: attenuation of a sine wave with depth δ13C and δD of CH4. For δD of CH4 this uncertainty is smaller than the typical measurement uncertainty, as well asRapid atmospheric variations are attenuated in the firn, and the atmospheric signal observed in the fiBrnrä(unlich et al,. recorded in the ice with a reduced amplitudSepa( hni et a,l. 2001; Mischler et al,. 2009). However, for bothδ13CO2 and 2003). To compare the firn smoothing effect between the δ13CH4 the model disagreement exceeds the typical instru-models, we force them with a sinusoidal atmospheric vari- mental precision. At this moment we have no objective wayation in CO2 that has a 15 yr period. Due to diffusion, the of determining which of the models predicts the DF correctly. amplitude of the signal decreases with depth in the firn. For This additional uncertainty should therefore be recognisedmathematical convenience we impose two separate atmo- when interpreting ice core and firn gas isotopic records ofspheric scenarios: periods with rapid atmospheric variations. Sites with a thick LIZ, combined with a low accumulation rate, allow for re- [C1](t,z= 0)= 1+ sin(2πt/15) constructing atmospheric mixing ratios further back in time. [C2](t,z= 0)= 1+ cos(2πt/15) However, a thick LIZ also means a larger uncertainty in DF (and thereby isotopic reconstruction) due to the poorly con-wheret is the time in years. Note that both signals are at- strained LIZ transport. tenuated to the same degree, and thatπt/h2e phase angle Molecular diffusion leads to DF between isotopologues between them is preserved. To get the signal amplitude with because their diffusion coefficients are different. Advection,depth, we combine them in the following way: convection and dispersive mass transfer do not discriminate √ |[C](t,z)| = ([C ](t,z)−1)2+ ([C ](t,z)−1)2. (5) between isotopologues, and consequently do not fraction- 1 2 ate the gas mixture. The observed model disagreement origT-he outcome is shown in Fig9.b. We observe that signal at- inates in the way in which LIZ transport is parameterised tenuation in the LIZ is stronger than in the DZ. This is be- (Fig.5c). The CIC and OSU models both change from purelycause transport in the DZ is efficient; the mean2CaOge at molecular diffusion in the DZ to purely eddy diffusion in the lock-in depth is around 9 yr, which is less than the period the LIZ. Once the molecular diffusion approaches zero (andof the atmospheric oscillation. The mixing ratios in most of DeddyDX), the DF ceases to increase with depth giving athe DZ will be “in phase” with the atmospheric signal, giv- horizontal line. The LGGE-GIPSA model, on the other hand, ing little attenuation. By comparison, air in the LIZ is nearly uses molecular diffusion all the way down to the close-off stagnant. It takes around 60 yr to be transported over 15 m, depth, giving a continued fractionation with depth. One waycorresponding to 4 periods of the sinusoid. As air of differ- the DF uncertainty could be incorporated in future studies isent ages (i.e. different phases of the atmospheric cycle) gets by testing different parameterisations for LIZ diffusion, or by mixed, the signal amplitude is reduced. using transfer functions generated by different models. Within the DZ the models agree well. In the LIZ there are The Law Dome-based records δo1f3CO2 (Francey et a,l. two intervals (65< z < 68 m and 74< z < 78 m) where the 1999) andδ13CH4 (Ferretti et al,.2005) were both corrected CSIRO model has significantly less attenuation with depth. for the effects of DF. The correction should not be affectedThese intervals correspond to the depths where the recon- much by our findings here for several reasons. Most im-structed CSIRO diffusivity goes to zero (Fi5gb.). On the US portantly, the accumulation rate (and therefore the advectiveborehole this effect is even more pronounced (Supplement). flux) at the Law Dome DE08 and DE08-2 sites is very high. In Sect.4.1 we touched upon the question of whether or This strongly reduces DF in the LIZ. Furthermore, the LIZ is not there is continued diffusion in the LIZ; our diffusivity relatively thin (Trudinger, 2001). The study byFerretti et al. reconstructions indicate that there is at NEEM. Scenario II (2005) combines ice core samples from the aforementionedshows that this remnant LIZ diffusivity causes very strong high accumulation Law Dome sites with firn air measure- attenuation of the atmospheric signal. Therefore, the ques- ments from the DSSW20K site, which has a rather low accu-tion of LIZ diffusivity has important consequences for the mulation rate∼( 0.16 m yr−1 ice equivalent). After DF cor- temporal resolution at which we can hope to reconstruct at- rection the results agree well. Other published firn air and icemospheric variations back in time using ice cores. core isotopic records that cover the recent anthropogenic in- crease have also been corrected for the effects of DF (e.g4..4.3 Scenario III: balance of transport fluxes Bräunlich et al,. 2001; Sugawara et a,l.2003; Röckmann et al., 2003; Sowers et a,l.2005; Bernard et a,l.2006; Ishi- The total gas flux is the sum of the advective, diffusive and jima et al., 2007; Mischler et al,. 2009). Until the true nature convective fluxes. Within the DZ, molecular diffusion over- of LIZ mixing is better understood, there is no way of know- whelms the other transport mechanisms; in this scenario we ing whether the DF magnitude calculated in these studies isbring the different fluxes into balance by considering a hypo- correct. Our findings imply that the uncertainty in these at- thetical gasX with a very small diffusion coefficienDt 0X = mospheric reconstructions might have been underestimated0..025D0CO , and a masMs −3 −1 X =Mair+1×10 kg mol . The2 relative strength of the different mechanisms is shown by the Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4273 gravitational enrichment with depth. This is the only of the sion term is neglected in the CSIRO modTerlu(dinger et al,. four scenarios where gravity is turned on in the model runs. 1997), which leads to a gas age that equals the ice age (since Figure9c shows the model output together with the true wair = wice). Note also that the observed difference is not re- gravitational enrichment as measured1f5oNr 2. In thermody- lated to the choice of reference frame (static vs. moving). namic equilibrium the enrichment is given by the baromet- In the LIZ air transport is dominated by advection, so this ric equation:δgrav(z)= [exp(1Mgz/RT )−1] ≈1Mgz/RT is where we expect to see the effect of the back flow most (Sowers et a,l.1992). The advective and convective fluxes clearly. Indeed we see that in the LIZ the CSIRO model di- drive the air column out of equilibrium reducinδggrav(z) be- verges most strongly from the other models, e.g. in the total low the barometric case. For gaXsthe diffusive flux is not reconstructed diffusivity (Fig5.b) and the diffusive smooth- large enough to oppose the other two, leading to a δfignraavl ing of scenario II (Fig.9b). which is about half that o1f5N2. In this way the enrichment Neglecting the back flow term corresponds to a physical with depth represents the balance between diffusion on onemodel of the firn where all the air is transported downwards hand, and advection and convection on the other. with the ice layers. This would be possible if, e.g. the denser The differences in implementation of the convective mix- firn layers succeed in completely sealing off the air below ing are clearly visible. The LGGE-GIPSA model uses no (Martinerie et al,.1992; Schwander et a,l.1993). In the case gravitational enrichment in the top 4 m, followed by a tuned such impermeable layers are present in the LIZ, pressure will eddy diffusion. The convective mixing in the other mod- build up in the pore space below. At some of the Law Dome els, which all use the parameterisation Kbaywamura et al. sites (e.g. DSSW20K) outgassing was observed from freshly (2006), penetrates much deeper into the firn. In the DZ theredrilled boreholes, which could be caused by the venting of are also clear differences; here the gravitational slope mostlysuch pressurised pores to the atmosphere. However, at these represents the balance between diffusion and advection. Thseites the effects of wind pumping could not be excluded as CSIRO has a stronger advection term due to the absence othf e origin of the observed air flow. At NEEM, no outgassing pore compression back flow (see below); this is reflected in afrom pressurised pores was observed. More evidence for the more reduced slope. A higher advective flux can partially beabsence of sealing layers comes from the observation of dif- compensated in the diffusivity tuning by a lower diffusivity, fusive mixing in the LIZ (Sect4. .1). Finally, the total air con- and vice versa. This explains why the diffusivity profile re- tent implied by the assumption of fully pressurised pores in constructed through the CSIRO model is slightly lower thanthe LIZ is incompatible with observed air content in mature those of the other models (see F5iga.). ice (Martinerie et al,.1992). For these reasons we believe that at NEEM the back flow 4.4.4 Scenario IV: advective transport due to pore compression can not be neglected when describ- ing the firn air transport. This makes the assumptions in the Pore closure traps air in bubbles, which are subsequently adC- SIRO model invalid, explaining why it obtains a higher vected downwards with the ice matrix. A downward flux of RMSD than most other models (Fig4).. Preliminary tests air in the open pores compensates for this removal process isnhow that by including the advective back flux in the CSIRO the firn-ice transition. In the last scenario we study hypotheti-model the RMSD reduces to 0.74. cal gasY , which is not subject to diffusion or convection, and is transported into the firn by the advective flux alone. Both eddy and molecular diffusion are set to zero. Its atmospheric5 Summary and conclusions history is a linear decrease by 1 ppm−1y,r until it reaches zero at the final model date. This leads to a mixing ratio ofWe presented a new multi-tracer method for characterising gasY in the firn that equals the age of the air at each depth.the firn air transport properties of a site. We applied the The results are shown in Fi9gd. . In the absence of diffusion, method to NEEM, Northern Greenland, by tuning six firn most models find an age of around 900 yr at the bottom ofair models to an ensemble of ten tracers. The firn air models the firn column. The CSIRO model stands out with a muchused in this study were able to fit the data within σa 1n-or- lower gas age. Differences in total air flux, and thereby themal distribution, meaning that the site can be modeled con- total air content in mature ice, affect the advective transport.sistently within the estimated uncertainties. Each of the refer- Also, artificial numerical diffusion in the models reduces the ence tracers constrains the firn profile differently through its calculated gas ages. unique atmospheric history and free air diffusivity, making There are significant differences in the way advectiveour multiple-tracer characterisation method an improvement transport is treated by the models. The key difference isover the commonly used single-tracer tuning. whether or not the model includes an upward air flow rel- Six different firn air transport models were tuned to the ative to descending ice layers to account for compression oNf EEM site in this fashion, allowing a direct model compari- open poresR(ommelaere et a,l1. 997; Sugawara et a,l2. 003). son. We find that all the models require a non-vanishing dif- Note that this term is smaller tha|wn ice|, and that the net air fusivity in the lock-in zone to reproduce the measured pro- flow in the open pores is always downwards. The compres-files of the reference gases. This is contrary to the commonly www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 4274 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM held notion that diffusive mixing is absent in the lock-in (FIST), France (IPEV, CNRS/INSU, CEA and ANR), Germany zone. A comparison between replicate boreholes located(AWI), Iceland (RannIs), Japan (NIPR), Korea (KOPRI), The 64 m apart suggests that this lock-in zone diffusion is sen-Netherlands (NWO/ALW), Sweden (VR), Switzerland (SNF), sitive to lateral inhomogeneities in firn stratigraphy. Further- United Kingdom (NERC) and the USA (US NSF, Office of Polar more, we find that despite the fact that the models reproducePrograms). C. Buizert would like to thank Bruce Vaughn (IN- the measured profiles equally well, the gas age distributionsSTAAR, University of Colorado) for his hospitality and support; V. V. Petrenko has been supported by the NOAA Postdoctoral they calculate can differ substantially. Often-used quantities,Fellowship in Climate and Global Change and NSF grants 0632222 such as the mean age and distribution width, differed amongand 0806387 (White); A. J. Orsi and J. P. Severinghaus were the models by up to 25 %. The modern d1aayge was calcu- supported by NSF 0806377; UEA acknowledges NERC for awards lated to be 18+23 −9 yr. NE/F021194/1, NE/F015585/1 and a PhD Studentship (Hogan); Finally, we introduced four diagnostic scenarios that arethe CSIRO contributions to this work were partly supported by the designed to probe specific aspects of the model physics, fromAustralian Climate Change Science Program, funded jointly by the which we can draw two main conclusions. First, we studiedDepartment of Climate Change and Energy Efficiency, the Bureau the consequences of the different ways in which the mod-of Meteorology and CSIRO; the French contribution to this study els implement advection, in particular whether they includewas further supported by CNRS through INSIS/PEPS-automatique a back flow term to account for pore compression. Includ-and INSU/LEFE programs. ing the back flow term significantly improves the fit to the Edited by: R. van de Wal data, suggesting this term can not be neglected, as is done in the CSIRO model. Second, we find that the effect of isotopic diffusive fractionation is poorly constrained by the models. Near the close-off depth, the calculated diffusive fractiona-References tion differs by 40 % between the models. These differences are related to the way in which lock-in zone diffusion is pa- Aris, R.: On the dispersion of a solute in a fluid flowing through a rameterised. Molecular diffusion leads to continued fraction- tube, Proc. R. Soc. Lond. A, 235, 67–77, 1956. ation, whereas dispersive transport (eddy diffusion) does notA. rnaud, L., Barnola, J. M., and Duval, P.: Physical modeling of the Further work is needed to elucidate the true nature of lock-in densification of snow/firn and ice in the upper part of polar ice zone diffusion. Meanwhile, when interpreting firn air and ice sheets, in: Physics of Ice Core Records, edited by Hondoh, T., 285–305, Hokkaido University Press, 2000. core isotopic records, an additional uncertainty needs to beAssonov, S. S., Brenninkmeijer, C. A. M., and Jockel, P.: The O-18 included to account for the model discrepancy found here. isotope exchange rate between firn air2CaOnd the firn matrix at three Antarctic sites, J. Geophys. Res.-Atmos., 110, D18310, doi:10.1029/2005jd00576, 92005. Supplementary material related to this article is Assonov, S. S., Brenninkmeijer, C. A. M.ö,cJkel, P., Mulvaney, R., available online at:http://www.atmos-chem-phys.net/12/ Bernard, S., and Chappellaz, J.: Evidence for a CO increase in 4259/2012/acp-12-4259-2012-supplement..zip the SH during the 20th century based on firn air samples from Berkner Island, Antarctica, Atmos. Chem. Phys., 7, 295–308, doi:10.5194/acp-7-295-200, 27007. Aydin, M., Saltzman, E. S., De Bruyn, W. J., Montzka, S. A., Butler, J. H., and Battle, M.: Atmospheric variability of methyl chloride AcknowledgementsW. e would like to thank Stephen Montzka, during the last 300 years from an Antarctic ice core and firn air, Bradley Hall and Geoff Dutton from NOAA-ESRL for providing Geophys. Res. Lett., 31, L02109d,oi:10.1029/2003gl01875, 0 halocarbon and S6Fatmospheric data and for their help on gas 2004. calibration scales; Ray Weiss and Ray Wang for their kind help onBattle, M., Bender, M., Sowers, T., Tans, P. P., Butler, J. H., Elkins, AGAGE atmospheric data and calibration scales; Paul Krummel J. W., Ellis, J. T., Conway, T., Zhang, N., Lang, P., and Clar- for providing NOAA/AGAGE(SIO) concentration ratios for ket, A. D.: Atmospheric gas concentrations over the past century halocarbon species based on flask and in-situ measurements; measured in air from firn at the South Pole, Nature, 383, 231– Bradley Hall for providing IHALACE results prior to their 235,doi:10.1038/383231a,01996. publication. Comments and suggestions by three anonymousBernard, S., R̈ockmann, T., Kaiser, J., Barnola, J.-M., Fischer, referees improved the manuscript considerably. This work ben- H., Blunier, T., and Chappellaz, J.: Constraints o2nONbud- efited greatly from numerous data available in public databases, get changes since pre-industrial time from new firn air and ice in particular: AGAGE h( ttp://cdiac.ornl.gov/ndps/alegage.h)t,ml core isotope measurements, Atmos. Chem. Phys., 6, 493–503, NOAA ESRL (http://www.esrl.noaa.gov/gm)d and EDGAR doi:10.5194/acp-6-493-200, 26006. (http://www.rivm.nl/edga)r. We thank the people involved in the Bräunlich, M., Aballanin, O., Marik, T., Jockel, P., Brenninkmei- acquiring, analysis, and accessibility of these data. NEEM is jer, C. A. M., Chappellaz, J., Barnola, J. M., Mulvaney, R., and directed and organised by the Center for Ice and Climate at the Sturges, W. T.: Changes in the global atmospheric methane bud- Niels Bohr Institute and US NSF, Office of Polar Programs. It get over the last decades inferred from C-13 and D isotopic anal- is supported by funding agencies and institutions in Belgium ysis of Antarctic firn air, J. Geophys. Res.-Atmos., 106, 20465– (FNRS-CFB and FWO), Canada (GSC), China (CAS), Denmark 20481,doi:10.1029/2001JD90019, 20001. Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4275 Butler, J. H., Battle, M., Bender, M. L., Montzka, S. A., Clarke, Program Australia 1999–2000, edited by: Tindale, N. W., Derek, A. D., Saltzman, E. S., Sucher, C. M., Severinghaus, J. P., and N., and Fraser, P. J., 42–53, Bureau of Meteorology, Melbourne, Elkins, J. W.: A record of atmospheric halocarbons during the Victoria, Australia, 2003. twentieth century from polar firn air, Nature, 399, 749–755, Fraser, P., Cunnold, D., Alyea, F., Weiss, R., Prinn, R., Simmonds, 1999. P., Miller, B., and Langenfelds, R.: Lifetime and emission es- Chen, N. H. and Othmer, D. F.: New Generalized Equation for Gas timates of 1,1,2-trichlorotrifluorethane (CFC-113) from daily Diffusion Coefficient, J. Chem. Eng. Data, 7, 37–41, 1962. global background observations June 1982 June 1994, J. Geo- Colbeck, S. C.: Air movement in snow due to windpumping, J. phys. Res.-Atmos., 101, 12585–1259d9o,i:10.1029/96JD0057,4 Glaciol., 35, 209–213, 1989. 1996. Conway, T. J., Tans, P. P., Waterman, S., L., and Thoning,Geller, L., Elkins, J., Lobert, J., Clarke, A., Hurst, D., But- K. W.: Evidence for interannual variability of the carbon-cycle ler, J., and Myers, R.: Tropospheric 6S:FObserved latitudi- from the National-Oceanic-and-Atmospheric-Administration nal distribution and trends, derived emissions and interhemi- Climate-Monitoring-and-Diagnostics-Laboratory Global-air- spheric exchange time, Geophys. Res. Lett., 24, 675–678, sampling-network, J. Geophys. Res.-Atmos., 99, 22831–22855, doi:10.1029/97GL0052,31997. 1994. Goujon, C., Barnola, J. M., and Ritz, C.: Modeling the densifica- Craig, H., Horibe, Y., and Sowers, T.: Gravitational Separation of tion of polar firn including heat diffusion: Application to close- Gases and Isotopes in Polar Ice Caps, Science, 242, 1675–1678, off characteristics and gas isotopic fractionation for Antarc- 1988. tica and Greenland sites, J. Geophys. Res.-Atmos., 108, 4792, Cunnold, D., Weiss, R., Prinn, R., Hartley, D., Simmonds, P., Fraser, doi:10.1029/2002jd00331, 92003. P., Miller, B., Alyea, F., and Porter, L.: GAGE/AGAGE mea- Hörhold, M. W., Kipfstuhl, S., Wilhelms, F., Freitag, J., and Fren- surements indicating reductions in global emissions of3CFCl zel, A.: The densification of layered polar firn, J. Geophys. Res.- and CC2l F2 in 1992–1994, J. Geophys. Res.-Atmos., 102, 1259– Earth., 116, F01001d,oi:10.1029/2009JF0016,320011. 1269,doi:10.1029/96JD0297,31997. Huber, C., Beyerle, U., Leuenberger, M., Schwander, J., Kipfer, R., Dlugokencky, E. J., Steele, L. P., Lang, P. M., and Masarie, K. A.: Spahni, R., Severinghaus, J. P., and Weiler, K.: Evidence for The growth-rate and distribution of atmospheric methane, J. Geo- molecular size dependent gas fractionation in firn air derived phys. Res.-Atmos., 99, 17021–17043, 1994. from noble gases, oxygen, and nitrogen measurements, Earth Dreyfus, G. B., Jouzel, J., Bender, M. L., Landais, A., Masson- Planet. Sci. Lett., 243, 61–73d,oi:10.1016/j.epsl.2005.12.0,36 Delmotte, V., and Leuenberger, M.: Firn processes δa1n5dN: 2006. potential for a gas-phase climate proxy, Quaternary Sci. Rev., 29I,shijima, K., Sugawara, S., Kawamura, K., Hashida, G., Morimoto, 28–42,doi:10.1016/j.quascirev.2009.10.0,122010. S., Murayama, S., Aoki, S., and Nakazawa, T.: Temporal vari- Etheridge, D. M., Steele, L. P., Langenfelds, R. L., Francey, R. J., ations of the atmospheric nitrous oxide concentration and its Barnola, J. M., and Morgan, V. I.: Natural and anthropogenic delta N-15 and delta O-18 for the latter half of the 20th century changes in atmospheric C2Oover the last 1000 years from air reconstructed from firn air analyses, J. Geophys. Res.-Atmos., in Antarctic ice and firn, J. Geophys. Res., 101, 4115–4128, 112(12), D03305d, oi:10.1029/2006jd00720, 82007. doi:10.1029/95JD0341,01996. Kawamura, K., Severinghaus, J. P., Ishidoya, S., Sugawara, S., Etheridge, D. M., Steele, L., Francey, R., and Langenfelds, R.: At- Hashida, G., Motoyama, H., Fujii, Y., Aoki, S., and Nakazawa, mospheric methane between 1000 AD and present: Evidence T.: Convective mixing of air in firn at four polar sites, Earth of anthropogenic emissions and climatic variability, J. Geophys. Planet. Sci. Lett., 244, 672–68d2o, i:10.1016/j.epsl.2006.02.0,17 Res.-Atmos., 103, 15979–1599d3o,i:10.1029/98JD0092,31998. 2006. Fabre, A., Barnola, J. M., Arnaud, L., and Chappellaz, J.: Determi-Kawamura, K., Parrenin, F., Lisiecki, L., Uemura, R., Vimeux, nation of gas diffusivity in polar firn: Comparison between ex- F., Severinghaus, J. P., Hutterli, M. A., Nakazawa, T., Aoki, perimental measurements and inverse modeling, Geophys. Res. S., Jouzel, J., Raymo, M. E., Matsumoto, K., Nakata, H., Mo- Lett., 27, 557–560d,oi:10.1029/1999GL01078, 02000. toyama, H., Fujita, S., Goto-Azuma, K., Fujii, Y., and Watan- Ferretti, D. F., Miller, J. B., White, J. W. C., Etheridge, D. M., abe, O.: Northern Hemisphere forcing of climatic cycles in Lassey, K. R., Lowe, D. C., Meure, C. M. M., Dreier, Antarctica over the past 360,000 years, Nature, 448, 912–916, M. F., Trudinger, C. M., van Ommen, T. D., and Lan- doi:10.1038/nature0601, 52007. genfelds, R. L.: Unexpected changes to the global methaneKeeling, R. F., Piper, S. C., Bollenbacher, A. F., and Walker, J. S.: budget over the past 2000 years, Science, 309, 1714–1717, Atmospheric CO2 records from sites in the SIO air sampling net- doi:10.1126/science.11151,923005. work, in: Trends: A Compendium of Data on Global Change, Francey, R., Allison, C., Etheridge, D., Trudinger, C., Enting, I., Carbon Dioxide Information Analysis Center, Oak Ridge Na- Leuenberger, M., Langenfelds, R., Michel, E., and Steele, L.: tional Laboratory, U.S. Department of Energy, Oak Ridge, Tenn., A 1000-year high precision record of delta C-13 in atmospheric USA, 2009. CO2, Tellus B, 51, 170–193d,oi:10.1034/j.1600-0889.1999.t01- Keeling, C. D., Piper, S. C., Whorf, T. P., and Keeling, R. F.: 1-00005.x, 1999. Evolution of natural and anthropogenic fluxes of atmospheric Francey, R. J., Steele, L. P., Spencer, D. A., Langenfelds, R. L., CO2 from 1957 to 2003, Tellus B, 63, 1–2d2o, i:10.1111/j.1600- Law, R. M., Krummel, P. B., Fraser, P. J., Etheridge, D. M., 0889.2010.00507,.x2010. Derek, N., Coram, S. A., Coopen, L. N., Allison, C. E., Porter, Landais, A., Caillon, N., Goujon, C., Grachev, A. M., Barnola, L., and Baly, S.: The CSIRO (Australia) measurement of green- J. M., Chappellaz, J., Jouzel, J., Masson-Delmotte, V., and house gases in the global atmosphere, in: Baseline Atmospheric Leuenberger, M.: Quantification of rapid temperature change www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012 4276 C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM during DO event 12 and phasing with methane inferred from air 791, 1992. isotopic measurements, Earth Planet. Sci. Lett., 225, 221–232M, atsunaga, N., Hori, M., and Nagashima, A.: Mutual diffusion doi:10.1016/j.epsl.2004.06.0,029004. coefficients of halogenated-hydrocarbon refrigerant-air systems, Landais, A., Barnola, J. M., Kawamura, K., Caillon, N., Delmotte, High Temp.-High Press., 25, 185–192, 1993. M., Van Ommen, T., Dreyfus, G., Jouzel, J., Masson-Delmotte,Matsunaga, N., Hori, M., and Nagashima, A.: Diffusion coefficients V., Minster, B., Freitag, J., Leuenberger, M., Schwander, J., Hu- of global warming gases into air and its component gases, High ber, C., Etheridge, D., and Morgan, V.: Firn air delta N-15 in Temp.-High Press., 30, 77–83, 1998. modern polar sites and glacial-interglacial ice: a model-data mis-Matsunaga, N., Hori, M., and Nagashima, A.: Measurements of the match during glacial periods in Antarctica?, Quaternary Sci. mutual diffusion coefficients of gases by the Taylor method (7th Rev., 25, 49–62, 2006. Report, measurements on the6S-aFir, SF6-N2, SF6-O2, CFC12- Laube, J. C., Martinerie, P., Witrant, E., Blunier, T., Schwan- N2, CFC12-O2, HCFC22-N2 and HCFC22-O2 systems), Trans. der, J., Brenninkmeijer, C. A. M., Schuck, T. J., Bolder, Jpn. Soc. Mech. Eng. B., 68, 550–555, 2002a. M., Röckmann, T., van der Veen, C.,öBnisch, H., Engel, Matsunaga, N., Hori, M., and Nagashima, A.: Measurements of A., Mills, G. P., Newland, M. J., Oram, D. E., Reeves, C. the mutual diffusion coefficients of gases by the Taylor method E., and Sturges, W. T.: Accelerating growth of HFC-227ea (8th Report, measurements on the HFC32-air, HCFC124-air, (1,1,1,2,3,3,3-heptafluoropropane) in the atmosphere, Atmos. HCFC125-air, HCFC143a-air, and HFC43-10mee-air systems), Chem. Phys., 10, 5903–5910d,oi:10.5194/acp-10-5903-20,10 Trans. Jpn. Soc. Mech. Eng. B., 68, 550–555, 2002b. 2010. Mischler, J. A., Sowers, T. A., Alley, R. B., Battle, M., Mc- Levin, I. and Kromer, B.: The tropospheric (C2)O-C-14 level in Connell, J. R., Mitchell, L., Popp, T., Sofen, E., and Spencer, mid-latitudes of the Northern Hemisphere (1959–2003), Radio- M. K.: Carbon and hydrogen isotopic composition of methane carbon, 46, 1261–1272, 2004. over the last 1000 years, Global Biogeochem. Cy., 23, GB4024, Levin, I., Hammer, S., Kromer, B., and Meinhardt, F.: Ra- doi:10.1029/2009GB00346, 02009. diocarbon observations in atmospheric 2C: ODetermining Montzka, S., Myers, R., Butler, J., Elkins, J., Lock, L., fossil fuel CO2 over Europe using Jungfraujoch observa- Clarke, A., and Goldstein, A.: Observations of HFC-134a in tions as background, Sci. Total Environ., 391, 211–216, the remote troposphere, Geophys. Res. Lett., 23, 169–172, doi:10.1016/j.scitotenv.2007.10.0,129008. doi:10.1029/95GL0359,01996. Levin, I., Naegler, T., Heinz, R., Osusko, D., Cuevas, E., Engel, A.,Montzka, S. A., Aydin, M., Battle, M., Butler, J. H., Saltzman, E. S., Ilmberger, J., Langenfelds, R. L., Neininger, B., Rohden, C. v., Hall, B. D., Clarke, A. D., Mondeel, D., and Elkins, J. W.: A Steele, L. P., Weller, R., Worthy, D. E., and Zimov, S. A.: The 350-year atmospheric history for carbonyl sulfide inferred from global SF6 source inferred from long-term high precision atmo- Antarctic firn air and air trapped in ice, J. Geophys. Res.-Atmos., spheric measurements and its comparison with emission invento- 109(11), D22302d, oi:10.1029/2004JD00468, 26004. ries, Atmos. Chem. Phys., 10, 2655–266d2o,i:10.5194/acp-10- Nydal, R. and L̈ovseth, K.: Carbon-14 Measurements in Atmo- 2655-2010, 2010. spheric CO2 from Northern and Southern Hemisphere Sites, Levin, I., Hammer, S., Eichelmann, E., and Vogel, F.: Verification 1962–1993, Carbon Dioxide Information Analysis Center, Oak of greenhouse gas emission reductions: The prospect of atmo- Ridge National Laboratory, Oak Ridge, Tennessee, 1996. spheric monitoring in polluted areas, Philos. T. Roy. Soc. A, 369,Prinn, R., Weiss, R., Fraser, P., Simmonds, P., Cunnold, D., Alyea, 1906–1924, 2011. F., O’Doherty, S., Salameh, P., Miller, B., Huang, J., Wang, R., MacFarling Meure, C., Etheridge, D., Trudinger, C., Steele, P., Hartley, D., Harth, C., Steele, L., Sturrock, G., Midgley, P., and Langenfelds, R., van Ommen, T., Smith, A., and Elkins, McCulloch, A.: A history of chemically and radiatively impor- J.: Law Dome CO2, CH4 and N2O ice core records ex- tant gases in air deduced from ALE/GAGE/AGAGE, J. Geophys. tended to 2000 years BP, Geophys. Res. Lett., 33, L14810 Res.-Atmos., 105, 17751–1779d2o, i:10.1029/2000JD90014, 1 doi:10.1029/2006gl02615, 2006. 2000. Manning, M. and Melhuish, W. H.: Atmospheric Delta 14C record Prinn, R., Huang, J., Weiss, R., Cunnold, D., Fraser, P., Sim- from Wellington, in: Trends: A Compendium of Data on Global monds, P., McCulloch, A., Harth, C., Reimann, S., Salameh, Change, Carbon Dioxide Information Analysis Center, Oak P., O’Doherty, S., Wang, R., Porter, L., Miller, B., and Krum- Ridge National Laboratory, U.S. Department of Energy, Oak mel, P.: Evidence for variability of atmospheric hydroxyl radicals Ridge, Tenn., USA, 1994. over the past quarter century, Geophys. Res. Lett., 32, L07809, Martinerie, P., Raynaud, D., Etheridge, D. M., Barnola, J. M., and doi:10.1029/2004GL02222, 82005. Mazaudier, D.: Physical and climatic parameters which influenceRasmussen, S., Andersen, K., Svensson, A., Steffensen, J., Vinther, the air content in polar ice, Earth Planet. Sci. Lett., 112, 1–13, B., Clausen, H., Siggaard-Andersen, M., Johnsen, S., Larsen, doi:10.1016/0012-821X(92)90002,-D1992. L., Dahl-Jensen, D., Bigler, M., Rothlisberger, R., Fischer, H., Martinerie, P., Nourtier-Mazauric, E., Barnola, J.-M., Sturges, W. Goto-Azuma, K., Hansson, M., and Ruth, U.: A new Greenland T., Worton, D. R., Atlas, E., Gohar, L. K., Shine, K. P., ice core chronology for the last glacial termination, J. Geophys. and Brasseur, G. P.: Long-lived halocarbon trends and budgets Res.-Atmos., 111, D06102d,oi:10.1029/2005JD00607, 29006. from atmospheric chemistry modelling constrained with mea- Reimer, P., Baillie, M., Bard, E., Bayliss, A., Beck, J., Bertrand, surements in polar firn, Atmos. Chem. Phys., 9, 3911–3934, C., Blackwell, P., Buck, C., Burr, G., Cutler, K., Damon, P., Ed- doi:10.5194/acp-9-3911-20,029009. wards, R., Fairbanks, R., Friedrich, M., Guilderson, T., Hogg, A., Massmann, J. and Farrier, D. F.: Effects of atmospheric pressures Hughen, K., Kromer, B., McCormac, G., Manning, S., Ramsey, on gas transport in the vadose zone, Water Resour. Res., 28, 777– C., Reimer, R., Remmele, S., Southon, J., Stuiver, M., Talamo, Atmos. Chem. Phys., 12, 42594–277, 2012 www.atmos-chem-phys.net/12/4259/2012/ C. Buizert et al.: Gas transport in firn: multiple-tracer characterisation for NEEM 4277 S., Taylor, F., van der Plicht, J., and Weyhenmeyer, C.: IntCal04Sowers, T., Bender, M., Raynaud, D., and Korotkevich, Y. S.: Delta- terrestrial radiocarbon age calibration, 0-26 cal kyr BP, Radio- N-15 of N2 in air trapped in polar ice – a tracer of gas-transport in carbon, 46, 1029–1058, 2004. the firn and a possible constraint on ice age-gas age-differences, Röckmann, T., Kaiser, J., and Brenninkmeijer, C. A. M.: The iso- J. Geophys. Res.-Atmos., 97, 15683–15697, 1992. topic fingerprint of the pre-industrial and the anthropogen2icON Sowers, T., Bernard, S., Aballain, O., Chappellaz, J., Barnola, source, Atmos. Chem. Phys., 3, 315–32d3o,i:10.5194/acp-3- J., and Marik, T.: Records of the delta C-13 of atmo- 315-2003, 2003. spheric CH4 over the last 2 centuries as recorded in Antarc- Rommelaere, V., Arnaud, L., and Barnola, J. M.: Reconstructing tic snow and ice, Global Biogeochem. Cy., 19, GB2002, recent atmospheric trace gas concentrations from polar firn and doi:10.1029/2004GB00240, 82005. bubbly ice data by inverse methods, J. Geophys. Res.-Atmos.S, pahni, R., Schwander, J.,üFclkiger, J., Stauffer, B., Chappellaz, 102, 30069–30083d,oi:10.1029/97JD0265,31997. J., and Raynaud, D.: The attenuation of fast atmospheri4c CH Schwander, J.: The transformation of snow to ice and the occlusion variations recorded in polar ice cores, Geophys. Res. Lett., 30, of gases, in: The Environmental record in glaciers and ice sheets, 1571,doi:10.1029/2003GL01709, 32003. edited by: Oescher, H. and Langway, C., 53–67, John Wiley, NewStuiver, M. and Polach, H. A.: Reporting of C-14 data – discussion, York, USA, 1989. Radiocarbon, 19, 355–363, 1977. Schwander, J., Stauffer, B., and Sigg, A.: Air mixing in firn and the Sturrock, G. A., Etheridge, D. M., Trudinger, C. M., Fraser, P. J., age of the air at pore close-off, Ann. Glaciol., 10, 141–145, 1988. and Smith, A. M.: Atmospheric histories of halocarbons from Schwander, J., Barnola, J. M., Andrie, C., Leuenberger, M., Ludin, analysis of Antarctic firn air: Major Montreal Protocol species, J. A., Raynaud, D., and Stauffer, B.: The age of the air in the firn Geophys. Res.-Atmos., 107, 476d5o,i:10.1029/2002JD00254, 8 and the ice at Summit, Greenland, J. Geophys. Res.-Atmos., 98, 2002. 2831–2838, 1993. Sugawara, S., Kawamura, K., Aoki, S., Nakazawa, T., and Schwander, J., Sowers, T., Barnola, J. M., Blunier, T., Fuchs, A., Hashida, G.: Reconstruction of past variations of delta C-13 and Malaize, B.: Age scale of the air in the summit ice: Impli- in atmospheric CO2 from its vertical distribution observed in cation for glacial-interglacial temperature change, J. Geophys. the firn at Dome Fuji, Antarctica, Tellus B, 55, 159–169, Res.-Atmos., 102, 19483–1949d3o,i:10.1029/97JD0130,91997. doi:10.1034/j.1600-0889.2003.0002,32.x003. Severinghaus, J. P. and Battle, M. O.: Fractionation of gases in polaTr rudinger, C.: The carbon cycle over the last 1000 years inferred lee during bubble close-off: New constraints from firn air Ne, from inversion of ice core data, Ph.D. thesis, 2001. Kr and Xe observations, Earth Planet. Sci. Lett., 244, 474–500,Trudinger, C. M., Enting, I. G., Etheridge, D. M., Francey, R. J., doi:10.1016/j.epsl.2006.01.0,32006. Levchenko, V. A., Steele, L. P., Raynaud, D., and Arnaud, L.: Severinghaus, J. P. and Brook, E. J.: Abrupt climate change at the Modeling air movement and bubble trapping in firn, J. Geophys. end of the last glacial period inferred from trapped air in polar Res.-Atmos., 102, 6747–676d3o, i:10.1029/96JD0338,21997. ice, Science, 286, 930–93d4o, i:10.1126/science.286.5441.9,30 Trudinger, C. M., Etheridge, D. M., Rayner, P. J., Enting, I. G., 1999. Sturrock, G. A., and Langenfelds, R. L.: Reconstructing atmo- Severinghaus, J. P., Grachev, A., and Battle, M.: Ther- spheric histories from measurements of air composition in firn, J. mal fractionation of air in polar firn by seasonal tem- Geophys. Res.-Atmos., 107, 478d0o,i:10.1029/2002JD00254, 5 perature gradients, Geochem. Geophy. Geosy., 2, 1048, 2002. doi:10.1029/2000GC00014, 26001. Trudinger, C. M., Etheridge, D. M., Sturrock, G. A., Fraser, P. J., Severinghaus, J. P., Grachev, A., Luz, B., and Caillon, N.: A method Krummel, P. B., and McCulloch, A.: Atmospheric histories of for precise measurement of argon 40/36 and krypton/argon ra- halocarbons from analysis of Antarctic firn air: Methyl bromide, tios in trapped air in polar ice with applications to past firn methyl chloride, chloroform, and dichloromethane, J. Geophys. thickness and abrupt climate change in Greenland and at Siple Res.-Atmos., 109, D22310d,oi:10.1029/2004JD00493, 2004. Dome, Antarctica, Geochim. Cosmochim. Ac., 67, 325–343, Tschumi, J. and Stauffer, B.: Reconstructing past atmospher2ic CO doi:10.1016/S0016-7037(02)00965, 2-1003. concentration based on ice-core analyses: open questions due to Severinghaus, J. P., Albert, M. R., Courville, Z. R., Fahnestock, in-situ production of CO2 in the ice, J. Glaciol., 46, 45–53, 2000. M. A., Kawamura, K., Montzka, S. A., Muhle, J., Scambos, Witrant, E., Martinerie, P., Hogan, C., Laube, J. C., Kawamura, T. A., Shields, E., Shuman, C. A., Suwa, M., Tans, P., and Weiss, K., Capron, E., Montzka, S. A., Dlugokencky, E. J., Etheridge, R. F.: Deep air convection in the firn at a zero-accumulation D., Blunier, T., and Sturges, W. T.: A new multi-gas constrained site, central Antarctica, Earth Planet. Sci. Lett., 293, 359–367, model of trace gas non-homogeneous transport in firn: evaluation doi:10.1016/j.epsl.2010.03.0,023010. and behavior at eleven polar sites, Atmos. Chem. Phys. Discuss., Smith, A. M., Levchenko, V. A., Etheridge, D. M., Lowe, D. C., 11, 23029–23080d,oi:10.5194/acpd-11-23029-20,121011. Hua, Q., Trudinger, C. M., Zoppi, U., and Elcheikh, A.: In search of in-situ radiocarbon in Law Dome ice and firn, in: 8th Interna- tional Conference on Accelerator Mass Spectrometry, 610–622, Elsevier Science Bv, Vienna, Austria, 1999. www.atmos-chem-phys.net/12/4259/2012/ Atmos. Chem. Phys., 12, 4245297–7, 2012