Glacial–interglacial dynamics of Antarctic ﬁrn columns: comparison between simulations and ice core air- δ 15 N measurements

. Correct estimation of the ﬁrn lock-in depth is es-sential for correctly linking gas and ice chronologies in ice core studies. Here, two approaches to constrain the ﬁrn depth evolution in Antarctica are presented over the last deglaciation: outputs of a ﬁrn densiﬁcation model, and measurements of δ 15 N of N 2 in air trapped in ice core, assuming that δ 15 N is only affected by gravitational fractionation in the ﬁrn column. Since the ﬁrn densiﬁcation process is


Introduction
Antarctic ice cores have provided outstanding records of past changes in climate and atmospheric composition (e.g. Jouzel et al., 2007;Loulergue et al., 2008;Lüthi et al., 2008;Schilt et al., 2010). However, a precise evaluation of the phase relationship between changes in local surface temperature and atmospheric composition remains challenged by the fact that air is trapped only at the bottom of the firn, a 60-120 m permeable layer below the surface where snow progressively densifies into ice. This leads to air becoming trapped at the bubble lock-in depth (LID) that is surrounded by ice as old as several hundred to up to ∼ 5500 yr in the case of central Antarctic sites, such as for the European Project for Ice Coring in Antarctica Dome C (EDC) site Fig. 1). The gas-ice offset is characterized by the ice age-gas age difference at a given depth, denoted age. Alternatively, it can be characterized by the depth difference in the ice core record between gas and ice of a given age, noted depth and given in actual physical metres. Constraining the firn structure is crucial to accurately estimate age and depth and reduce uncertainties in ice and gas chronologies, in particular for clarifying the exact timing between CO 2 concentration and Antarctic surface temperature during deglaciations (Fischer et al., 1999;Caillon et al., 2001;Monnin et al., 2001;Pedro et al., 2011;Shakun et al., 2012;Parrenin et al., 2013).
Firn densification models have been specifically developed to build ice core gas chronologies, which require estimating age or depth (e.g. Blunier et al., 2004;Bender et al., 2006;Loulergue et al., 2007). They assume a homogeneous snow material where the density profile, and thus the LID, is mainly dependent on the temporal scenarios of accumulation rate, surface temperature and surface density (Herron and Langway, 1980;Pimienta, 1987;Schwander et al., 1993;Arnaud et al., 2000;Goujon et al., 2003).
The isotopic composition of nitrogen (δ 15 N of N 2 , hereafter δ 15 N) in air trapped in ice core also provides information on past firn depth. It is modified by physical process such as thermal and gravitational fractionations. Large δ 15 N anomalies develop during episodes of rapid climatic changes, such as those recorded in Greenland ice cores, due to a temperature gradient developing in the entire firn (e.g. Severinghaus et al., 1998). While seasonal temperature gradients can occur in the top ∼ 10 m of the firn (e.g. , glacial-interglacial Antarctic surface temperature changes as inferred from water stable isotope records are too slow to create a vertical temperature gradient in the whole firn that would lead to a significant thermal δ 15 N anomaly in the trapped air bubbles. In other words, δ 15 N variations are here considered as only caused by gravitational fractionation, which leads to an enrichment of the trapped air in heavy isotopes proportional to the diffusive column height (DCH). In δ notation (given in ‰), the gravitational fractionation follows the barometric equation (Eq. 1) where m is the mass difference (kg mol −1 ; for the case of δ 15 N, it is the mass difference between 15 N and 14 N), g is the gravitation acceleration (m s −2 ), z the DCH (m), R the gas constant (J K −1 mol −1 ) and T the firn temperature (K). It results from Eq. (1) that the gravitational fractionation of δ 15 N in the firn is influenced directly by the mean firn temperature and by any factor that changes the DCH such as the firn temperature, the surface accumulation rate, the initial snow density and the firn permeability. The convective zone, in the upper part of the firn, is characterized by convective mixing that overwhelms molecular diffusion and prevents isotopic fractionation. Assuming that this convective zone is negligible, the DCH provides an estimate of the LID, i.e. the depth where the gas diffusion becomes negligible. Sites where firn air studies have been conducted so far are characterized by a convective zone spanning from 0 m to up to 20 m depth Severinghaus et al., 2010).
The "state-of-the-art" firn densification models have been evaluated against modern firn air δ 15 N observations spanning a range of mean annual temperatures at various Antarctic and Greenlandic sites (from −19 to −55.5 • C for surface temperature and from 2.2 to 140 cm water equivalent per year, water eq yr −1 , for the accumulation rate; Goujon et al., 2003;Landais et al., 2006). These models are also able to reproduce the glacial LID inferred from δ 15 N records from various Greenland ice cores (e.g. Landais et al., 2004;NorthGRIP community members, 2004;Huber et al., 2006), and from the Antarctic Byrd ice core (Sowers et al., 1992). Glacial climatic conditions at these sites are within the present range of surface parameters for which the models have been evaluated, e.g. LGM mean surface temperature of about −43 and −52 • C, and LGM mean accumulation rate of ∼ 5 and ∼ 6 cm water eq yr −1 for Byrd (Blunier et al., 1998) and NorthGRIP, respectively (Johnsen et al., 2001).
Accumulation rates are reduced during glacial periods. Considering a densification process at a constant speed, a smaller accumulation rate leads to a decrease of the LID. By contrast, colder conditions induce a slower densification process, leading to a firn deepening. In several Antarctic sites characterized by low accumulation rates (Vostok, EDC, EDML, Dome F), firn models predict that the LID should decrease from glacial to interglacial periods. Thus, firn models predict that the temperature effect dominates over the effect of accumulation rate on the LID evolution at the glacialinterglacial scale. This is the opposite to the LID evolution inferred from δ 15 N measurements (Fig. 2). This modeldata δ 15 N mismatch has been largely discussed by Caillon et al. (2001), Landais et al. (2006) and Dreyfus et al. (2010).
First, as firn models have been fitted onto observations under present-day climate, the extrapolation of their results outside the range of observations may be incorrect . Alternatively, Landais et al. (2006) proposed that the relationships between water stable isotopes, temperature and accumulation used to produce climatic scenarios to force firn models may be incorrect . With these two potential explanations, the model-data δ 15 N mismatch relies on the common assumption that the physics of firnification models is globally correct, and that firn model outputs can be reconciled with δ 15 N data after adjustments for the last deglaciation. The altitude, the distance from the coast, the mean annual surface 823 temperature, T 0 (°C) and the accumulation rate (in water equivalent per year), A 0 are 824 indicated for the ice core sites discussed in this study (Mulvaney et al., 2000;EPICA-825 community-members, 2004;EPICA community members, 2006;Loulergue et al., 2007;;826 Parrenin et al., 2007a;Buiron et al., 2011;Stenni et al. 2011). Note that the surface 827 accumulation rates of BI and JRI have been deduced from the ice flow model and adjusted to 828 find the best agreement with stratigraphic markers (see Appendix A1 for details). For JRI, A 0 829 is close to the value reported by Abram et al. (2011) equal to 62 ± 1.4 cm of water equivalent 830 per year and deduced from annual layer counting. 831 832 Fig. 1. Location of the Antarctic ice cores where δ 15 N measurements have been obtained for the last deglaciation. The altitude, the distance from the coast, the mean annual surface temperature -T 0 ( • C) -and the accumulation rate (in water equivalent per year) -A 0 -values are indicated for the ice core sites discussed in this study (Mulvaney et al., 2002;EPICA community members, 2004EPICA community members, , 2006Parrenin et al., 2007a;Buiron et al., 2011;Stenni et al., 2011). Note that the surface accumulation rates of BI and JRI have been deduced from the ice flow model and adjusted to find the best agreement with stratigraphic markers (see Appendix A for details). For JRI, A 0 is close to the value reported by Abram et al. (2011) equal to 62 ± 1.4 cm of water equivalent per year and deduced from annual layer counting. of the forcing scenarios and/or of the modelled influences of accumulation rate and temperature on the firn LID, especially for inland sites characterized by low temperatures and accumulation rates (here called Hypothesis A). Second, Caillon et al. (2001) and Dreyfus et al. (2010) suggested that the discrepancy between measured δ 15 N and modelled δ 15 N is not due to errors in the firn model or climate forcing scenario but rather to the presence of a deep convective zone under glacial conditions (here called Hypothesis B) linked to an increased firn permeability in periods of low accumulation rate (Courville et al., 2007). Indeed, the existence of a deep convective zone would reduce the measured δ 15 N Model of Herron and Langway (1980) (green diamonds); (ii) Model of Arnaud et al. 6 (2000)(blue diamonds); (iii) Model of Goujon et al. (2003)(red curve). Error bars on 7 simulated LID represent a 30% uncertainty on the past accumulation rate estimate. Similarly 8 with the grey area for the Goujon simulated curve. An alternative accumulation rate scenario 9 deduced from volcanic stratigraphic markers  is also used to force the 0 Goujon model (purple curve). Diffusive Column Height (DHC) deduced from δ 15 N 1 measurements (opened black diamond curve), Landais et al., 2006) accounting for a 20% 2 uncertainty on the temperature estimate (black dotted curves).  Goujon et al. (2003) (red curve). Error bars on simulated LID represent a 30 % uncertainty on the past accumulation rate estimate. Similarly, the grey area represents uncertainty for the Goujon simulated curve. An alternative accumulation rate scenario deduced from volcanic stratigraphic markers  is also used to force the Goujon model (purple curve). Diffusive column height (DCH) deduced from δ 15 N measurements (opened black diamond curve,  accounting for a 20 % uncertainty on the temperature estimate (black dashed curves). (b) EDML δD profile (Stenni et al., 2010). levels through the reduction of the diffusive zone, but not the modelled δ 15 N since the later is calculated in the firn model as a function of the modelled LID. Third, Hörshold et al. (2012) have demonstrated recently that the snow/ice impurity content (e.g. insoluble dust or Ca 2+ concentrations) may have a significant impact on the densification process (hereafter referred to as Hypothesis C), with a decrease in firn depth at increasing impurity levels. At the moment, no parameterization of this effect is available for implementation in firn densification models.
Here, we present published and new measurements of δ 15 N and simulations of firn densification over the last deglaciation for five Antarctic sites: Dome C (EPICA Dome C ice core, EDC), Kohnen Station (EPICA Dronning Maud Land ice core, EDML), Talos Dome (TALDICE ice core), Berkner Island (BI ice core) and James Ross Island (JRI ice core). These sites offer surface climatic conditions spanning a very large range of accumulation rates and temperatures. Each of these sites provides also a specific case due to intersite differences in latitude (and therefore insolation), elevation and distance to the nearest open ocean (Fig. 1). During glacial periods, the coastal or semi-coastal sites are expected to undergo surface temperature and accumulation rates that fall within the densification model empirical validity range. Each of these sites is also characterized by a specific magnitude of glacial-interglacial changes in local insoluble dust concentration (Ruth et al., 2008;Albani et al., 2012;Lambert et al., 2012), allowing us to test Hypothesis C. Using these new datasets, together with water isotope profiles and tests conducted with firn models, we investigate and discuss the different hypotheses presented above and their ability to explain the past firn structure dynamics for semi-coastal and coastal Antarctic sites.
In the following the analytical method for δ 15 N measurements is summarized (Sect. 2). Simulations of firn densification during the last glacial-interglacial transition are conducted for the five ice core sites and discussed (Sect. 3). The new JRI, BI and TALDICE δ 15 N profiles are described and compared with existing profiles from the EPICA ice cores (EDML and EDC) and firn modelling results (Sect. 4). The mechanisms governing past firn structure evolution in Antarctica are finally discussed (Sect. 5).

Measuring δ 15 N from trapped air in ice:
analytical procedure Here, we complement existing ice core δ 15 N data from EDC and EDML sites (Dreyfus et al., 2010;Landais et al., 2006) by additional measurements on the EDML ice core and new data measured on the recently drilled BI, TALDICE and JRI ice cores (Fig. 1). New air isotopic measurements were performed at the Laboratoire des Sciences du Climat et de l'Environnement between 2007 and 2011, during several measurement periods, using a melt-refreeze technique (Sowers et al., 1989;Landais et al., 2003) to extract fossil air from the ice (Table 1). Air samples were then analysed on a 10-collector Delta V Plus (ThermoElectron Corporation) isotope ratio mass spectrometer which allows simultaneous measurements of masses 28, 29, 30, 32, 33, 34, 36, 38, 40 and 44. Corrections for pressure imbalance and chemical interferences of CO 2 and δO 2 /N 2 were applied to improve the measurement precision following the procedure fully described in Severinghaus et al. (2001) and Landais et al. (2003). The analytical precision over a given measurement period is calculated as the pooled standard deviation of depth pairs  and is presented in Table 1. The pooled standard deviation for each dataset varies from 0.005 ‰ for the JRI dataset to up to 0.022 ‰ for the EDML dataset. It does not affect the following discussion because the amplitude of the δ 15 N variations considered here is much larger.

Firn densification models
We have used the sophisticated firnification model of Goujon et al. (2003), hereafter referred to as the Goujon model. This model has been built by implementation of the heat diffusion in the firn model of Arnaud et al. (2000), which relies on physical processes of pressure sintering to describe the relationships between density, surface temperature and accumulation rate, while some parameters have also been fitted onto density profiles. The Goujon model requires a depthage correspondence for thermal diffusion calculations. The calculation of the LID by the Goujon model requires setting a closed porosity percentage threshold to define when the air diffusion becomes negligible within the firn. The closed porosity is expressed as a function of the total porosity, and total gas measurements suggest that the closed porosity at the close-off is equal to 37 % of the total porosity (Goujon et al., 2003). However, firn air measurements at Vostok and Summit suggest that gas diffusion in the firn stops at lower ratios of total porosity. Based on those data, Goujon et al. (2003) define that the firn gas diffusion stops at a closed porosity ranging from 37 to 13 % of the total porosity. As the total porosity is calculated from the modelled firn density (Eqs. 9 and 10 in Goujon et al., 2003), the selection of a closed porosity threshold corresponds to the selection of a firn density for bubble close off, and thus strongly relies on a correct estimate of the density profile in the firn. No data are available to precisely assess the closed porosity threshold that, for each site, defines the LID. Sensitivity tests were conducted using the Goujon model with several closed porosities within the 37-13 % range, showing a very small impact on the modelled δ 15 N values (less than 0.009 ‰ on average, not shown). We thus arbitrarily used for all our sites the same version of the Goujon model set with a closed porosity percentage of 21 %, as defined for the Vostok site (Goujon et al., 2003). Based on the LID simulations, MODEL-δ 15 N is then deduced from Eq. (1). Table 2 provides information on the gas and ice timescales onto which our new results have been transferred. The official published chronologies are used for the EDC, EDML and TALDICE ice cores Parrenin et al., 2007a;Buiron et al., 2011). Preliminary timescales for the BI and JRI ice core have been derived from a glaciological approach (Parrenin et al., 2007b) and coupled to chronological constraints derived from comparison of ice (δD) and gas records (CO 2 , CH 4 , δ 18 O of O 2 ) to the well-dated ice cores. Details are given in the Appendix A.

Temperature and accumulation scenarios as input parameters
Surface climatic condition scenarios used to force the firn densification models are deduced following the procedure described by Parrenin et al. (2007a). The past surface temperature and accumulation rates are both estimated from the water isotopic records for each ice core. Detailed equations are given in Table 2.
A 0 (cm of water eq yr −1 ) and T 0 (K) are, respectively, the surface accumulation rate and temperature for the present taken for each site, as given in Fig. 1. δD ( δ 18 O) corresponds to the difference between δD (δ 18 O) at a given depth and the present-day value, δD 0 (δ 18 O 0 ). Note that water isotopic profiles are corrected for the influence of vapour source changes using the mean ocean δ 18 O (Bintanja et al., 2005;Parrenin et al., 2007a). α D and α O (K ‰ −1 ) represent the spatial slope of the present-day isotopic thermometer, while the parameter β (‰ −1 ) controls the glacial-interglacial amplitude of the accumulation rate change. Alternatively to the use of Eq. (4) at EDML , Buiron et al. (2011) calculated a synthetic δD record from the TALDICE δ 18 O ice data through the following equation: assuming no change in deuterium excess. For EDC, EDML and TALDICE, the same values of α and β parameters are used as previously optimized to construct their official ice and gas chronologies Parrenin et al., 2007a;Buiron et al., 2011). For the JRI, the α value determined by Abram et al. (2011) is 0.1563 K ‰ −1 . As BI does not benefit yet from a local estimate of α, the classical spatial slope of 0.01656 K ‰ −1 is used (Lorius and Merlivat, 1977). For both sites we used the β values that enable obtaining the best agreement with both the ice flow constraints and the available stratigraphic constraints (see Appendix A). Table 2 summarizes the respective values of α and β parameters for these scenarios. All these equations rely on the assumption that the isotope-temperature relationship observed today spatially (and driven by distillation processes) remains valid for past changes (e.g. . It implies that surface and condensation temperatures co-vary, and requires limited precipitation intermittency biases or changes in moisture source conditions (for temperature estimates); this assumption has been challenged for warmer than present-day conditions, based on one atmospheric model (Sime et al., 2009). The uncertainty associated 988 E. Capron et al.: Glacial-interglacial dynamics of Antarctic firn columns Table 2. Information about each ice core: description of available chronologies and methods, α and β parameters relating surface temperature and accumulation rate to water stable isotopes, and estimated LGM surface conditions (temperature and accumulation rates).
Ice core site (official chronology name) Chronology available

Modelled δ 15 N variations
For all sites, the LGM MODEL-δ 15 N mean level is higher than the early Holocene (EH) MODEL-δ 15 N mean level (red curves in Figs. 3-4). This is particularly obvious for the EDC and EDML sites, while the amplitude of the MODEL-δ 15 N variation from the LGM to the EH is relatively reduced at JRI, TALDICE and BI. A greater gravitational fractionation during glacial time results from a deeper LID modelled under colder conditions. It illustrates that the Goujon model not only predicts that (i) the surface temperature increase is the dominant factor controlling the LID evolution during such a large climatic transition, for the sites characterized by the lowest accumulation rate, but (ii) also that stronger competition with the effect of accumulation rate occurs for the coastal sites. The opposite influences of temperature and accumulation rate on firnification processes are illustrated for the TALDICE case by comparing two simulations: (i) an "Acc MODEL-δ 15 N" curve, which represents MODELδ 15 N simulated in response to accumulation changes only, and (ii) a "Temp MODEL-δ 15 N" curve, simulated when considering only the effect of temperature change (Fig. 4b).
While the two factors have clearly opposite effects when considered individually, the total MODEL-δ 15 N curve is not simply the average of the two δ 15 N simulations considering each single factor. This result is also valid for the EDC, EDML and BI cases and is due to non-linear interactions because the accumulation rate influence is different for different temperature levels, and vice versa (see Fig. B1 and Appendix B).
The decrease of MODEL-δ 15 N over the last deglaciation is not monotonic for the EDC, BI, EDML and TALDICE sites. Indeed, we observe that MODEL-δ 15 N is increasing in parallel to δD at the start of the deglaciation and during the warming after the ACR (phases 1 and 3, respectively, on Fig. 4). Sensitivity tests performed on all sites show that the MODEL-δ 15 N increase observed during the start of the deglaciation, i.e. between the LGM and the ACR, should be attributed to the effect of accumulation rate (phase 1, Figs. 4 and B1, left panels; Appendix B).

Comparing modelled δ 15 N profiles with new δ 15 N measurements
For all ice core sites, including JRI, we confirm the overall model-data δ 15 N mismatch over glacial-interglacial variations, which was previously reported for central Antarctic ice cores (Kawamura, 2000;Caillon et al., 2001;Dreyfus et al., 2010;Figs. 3 and 4). The fact that this model-data mismatch is also depicted at JRI is a surprise because JRI surface climatic conditions (−14 • C of annual mean temperature, snow accumulation rate of 62 cm water eq yr −1 ; Abram et al., 2011) are warmer than for Greenland deep ice core sites, where firn models perform well for glacial-interglacial variations. Despite a general model-data mismatch on the absolute δ 15 N value at LGM, we still note relatively good agreement in the JRI ice core (Mulvaney et al., 2012). As a result, we cannot discuss the MODEL-δ 15 N 853

along the deglaciation and we focus only on the mean MODEL-δ 15 N levels for LGM (Last 854
Glacial Maximum) and EH (Early Holocene) climatic conditions. 855  Mulvaney et al., 2012), MODEL-δ 15 N (red, this study) and DATA-δ 15 N (blue, this study) over the time interval 7-30 ka. Note that the water stable isotope variation suggests an unrealistically fast deglaciation compared to all other Antarctic records, related to an unconformity present in the early deglacial interval in the JRI ice core (Mulvaney et al., 2012). As a result we cannot discuss the MODEL-δ 15 N along the deglaciation and we focus only on the mean MODEL-δ 15 N levels for LGM (Last Glacial Maximum) and EH (Early Holocene) climatic conditions. between MODEL-δ 15 N and DATA-δ 15 N trends over the two warming phases of LGM-ACR and ACR-EH (phases 1 and 3 on Fig. 4) recorded at BI, TALDICE, EDML and EDC, apart for phase 3 at BI. At TALDICE the agreement is even better since both trends and absolute values of MODEL-δ 15 N and DATA-δ 15 N are coherent over the three different phases of the deglaciation (Fig. 4b).
Coherent trends in the MODEL-δ 15 N and DATA-δ 15 N show that the MODEL-δ 15 N captures well some of the DATA-δ 15 N variations. Still, the DATA-δ 15 N at most of the sites shows a stronger variability of the DATA-δ 15 N than MODEL-δ 15 N. For example, at TALDICE, DATA-δ 15 N data increase by 0.080 ‰ over the LGM-ACR warming while MODEL-δ 15 N model increases by 0.030 ‰. In particular, the high-resolution measurements performed on the BI ice core reveal the largest millennial-scale variations so far measured in an Antarctic δ 15 N profile (Fig. 4c). These variations represent true variations in the DCH thickness as (i) they are significantly larger than the analytical error (less than 0.015 ‰) and (ii) each rapid increase/decrease is defined by several consecutive measurements. DATA-δ 15 N also exhibits significant variations which cannot be linked to any large variations in the water stable isotope profile, e.g. at EDML, DATA-δ 15 N decreases by 0.073 ‰ corresponding to a DCH thinning of ∼ 17 m in ∼ 1 ka at 19 ka and the fastest δ 15 N variation occurs at 11.2 ka with a δ 15 N increase of 0.090 ‰ in 170 yr (equivalent to a ∼ 20 m DCH increase).

Summary
This model-data comparison leads to three main conclusions: -During phases of the deglaciation with a significant increase in accumulation rate like the LGM to ACR period, the MODEL-δ 15 N trends derived from firn mod-elling are consistent with the DATA-δ 15 N measured for most Antarctic sites.
-Simulations predict significantly higher glacial MODEL-δ 15 N levels than the measured ones. The model-data mismatch is the strongest for sites characterized by a low accumulation rate.
-Larger levels of δ 15 N variability are depicted by measurements than simulations.

What controls glacial-interglacial changes in firn structure?
In the light of our new measurements and simulations, we now assess the three hypotheses given in the introduction to explain the observed model-data δ 15 N mismatch at EDC, EDML and BI. Hypotheses A (relationships between firn depth and accumulation or temperature) and B (convective zone) assume that the physics of the firnification model is generally appropriate, while Hypothesis C assumes that the mismatch is due to the effect of snow impurity content on densification, which is not implemented in firn models. Note that assessing the validity of the physics of the firn model relies on the comparison between DATA-δ 15 N and MODELδ 15 N in the absence of any convective zone (Hypothesis A), as summarized in the previous section. We include this information in our final discussion (Sect. 5.3).

Investigating the presence of a deep glacial convective zone at EDML
In order to disentangle between Hypotheses A and B for the EDML site, we compare the depth difference observed along the ice core record between two synchronous events recorded in the ice phase and the gas phase, respectively, seen as the  depth, obtained through three different methods (details are given in the Appendix C): -We deduce an estimate of the depth evolution based on the DATA-δ 15 N measurements. For that purpose we translate the EDML DATA-δ 15 N into a firn LID, hence assuming (1) that no significant convective zone affected the firn in the past, (2) that δ 15 N only reflects the gravitational settling, and (3) that the difference between the LID and the COD does not significantly affect the gas repartition during bubble close-off. Then, we account for the thinning due to ice flow and translate the firn equivalent DCH thickness to an ice equivalent DCH thickness (equation given in the Appendix C) to obtain the EDML DATA-δ 15 N-based depth (Fig. 5).
-We obtain a modelled depth estimate by combining the LID estimate from the Goujon firn model  together with an estimate of the thinning function under the assumption that there is no significant convective zone.
-We infer two independent empirical estimates of the depth during the Laschamp event (41.2 ± 1.6 ka on GICC05; Svensson et al., 2007) from the 10 Be records in the ice phase and from the CH 4 records in the gas phase (see Loulergue et al., 2007 and Appendix C for more details). Our empirical depth estimates are equal to 22.0 ± 2.5 m and 25.2 ± 2.5 m at about 1368-1407 m depth.
We observe that the modelled depth and the δ 15 N-based depth estimates are generally compatible within the uncertainty range ranging from 6.5 and 3 m over the depth interval 540-1410 m of the DATA-δ 15 N; Fig. 5). Over the depth interval 960-1110 m corresponding to the LGM time period, the uncertainty on the depth estimates translates into an uncertainty of about 10 m on the firn thickness. The average difference between the two depth estimates could suggest the presence of a convective zone of about 12 m deep during the LGM, which could explain why modelled depth (similarly, the LID) is larger than DATA-δ 15 N-based depth (similarly, the DCH). However, considering the associated uncertainties, we cannot reject the hypothesis of the absence of a deep convective zone either.
Thus, it is impossible to draw a firm conclusion on the presence and thickness of the convective zone at EDML during the LGM but, fortunately, we can provide a more precise evaluation of the "deep convective zone" hypothesis over the Laschamp event. Indeed, we compare the two empirical depth estimates with the depth estimate based on δ 15 N measurements equal to 26.8 ± 3 m over the 1360-1400 m depth interval. This value is slightly larger than the empirical depth estimates from gas CH 4 and ice 10 Be matching. The existence of a convective zone would lead to a δ 15 Nbased depth smaller than the modelled or empirically derived depth, which is opposite to our observation.
For now we cannot provide additional independent constraints on the evolution of depth along the EDML ice core such as proposed by Parrenin et al. (2013) for the EDC ice core. While our new results do not allow us to firmly state the thickness of the convective zone during the LGM, they show the absence of a convective zone at EDML at the time of the Laschamp event. Note that Parrenin et al. (2012) had ruled out a large glacial convective zone at EDC, using a similar depth-based approach.

A dust influence on firnification?
We test Hypothesis C by analysing the phase relationship between DATA-δ 15 N variations and changes in ice core dust concentration available from the EDML, TALDICE, BI and EDC ice cores on a depth scale. If the main control on density evolution is the impurity content, we could expect to observe large changes in DATA-δ 15 N at depths where large changes are recorded in markers of impurity content (such as Ca 2+ or the insoluble dust concentrations). For EDML, TALDICE, BI and EDC, we cannot observe any systematic visual link between changes in records of ice impurity content and DATA-δ 15 N variations on a depth scale (Fig. 4). In particular, the BI DATA-δ 15 N profile presents a large millennial-scale parallel to a regular decrease in dust concentration during the deglaciation. For EDML and EDC a rather clear anti-phase relationship is observed; however, it is not easy to separate any effect of impurity content from the impact of parallel changes in surface temperature and accumulation rates on the firn structure because glacialinterglacial changes in dust concentrations often strongly covary with Antarctic climate changes (e.g. Lambert et al., 2012). For time periods where large variations of DATAδ 15 N are measured without any concomitant variability in δD, no significant change in the impurity content is recorded either. Moreover, the fact that glacial and interglacial DATAδ 15 N levels measured on the TALDICE and BI ice cores are 992 E. Capron et al.: Glacial-interglacial dynamics of Antarctic firn columns approximately equal disfavours the dust hypothesis as the latter suggests a smaller LID during glacial time when impurity concentrations in snow are higher.
These new observations do not favour Hypothesis C as the major explanation for firn LID changes over the deglaciation. However, our study is limited by the temporal resolution of the available datasets and by the parallel variations observed between the impurity content and the climate variables. We are therefore unable to separate their effects on densification. Our analysis is strongly limited by the lack of understanding of the physics behind the impurity effect. In particular, the link between impurity content and density kinetics may not be linear, and thus a strict visual correlation may not be necessarily expected. Previous studies have highlighted a major role of dust on the modification of ice microstructure (through the pinning of grain boundaries, Durand et al., 2006). Thus, further investigations are required based on future high-resolution glacial dust concentration (or Ca 2+ ) and δ 15 N records, and also from an improved quantitative understanding of the links between dust concentrations, grain growth and metamorphism and densification processes.

Synthesis
Our study suggests that the physics of the firnification model is at least partly correct (Hypothesis A), but some processes controlling δ 15 N variations are still missing. We propose that the remaining mismatch between modelled and measured δ 15 N can be attributed to the following causes: 1. The process of firn deepening in response to deglacial accumulation rate increase is underestimated in the firnification model. The densification might be a more time-controlled phenomenon than a pressure-controlled phenomenon. Indeed, if the densification was only time controlled, age would be constant through time and the LID would be proportional to accumulation.
2. Inaccurate scenarios for past accumulation evolution are used to force the firn model and methods to estimate past accumulation rates need to be revised. Indeed, the water isotope-based approach to infer past accumulation rate might not be well suited for semi-coastal and coastal regions, where the atmospheric moisture content is probably not only controlled by local temperature but also by changes in cyclonic activity, changes in precipitation intermittency, moisture source conditions and distillation paths at synoptic and seasonal scales (van Ommen et al., 2004;Monnin et al., 2004;Landais et al., 2006;Masson-Delmotte et al., 2008). In particular, the site of BI has encountered very particular site-specific climatic and glaciological changes (R. Mulvaney, personal communication, 2013). Past changes of local ice sheet topography could have had significant impacts on (i) atmospheric circulation and elevation-accumulation-temperature-water stable iso-tope relationships, and (ii) ice flow, layer thinning and inferred accumulation rates. At EDML, an accumulation rate scenario inferred from volcanic signature matching with the EDC ice core produces accumulation variations during the glacial period that are not linked to any variations in the water isotopic profiles . Further investigations are required in particular to test whether the increase in Antarctic accumulation rates is underestimated over the deglaciation, especially from the ACR to the EH, when the disagreement is the largest between modelled and measured δ 15 N.
3. The heterogeneous behaviour of the firn structure evolution over the last deglaciation from one site to the other is also likely to result from strong competition and/or compensation between several of the discussed mechanisms which is specific to each studied site.

Conclusions and perspectives
In this paper we have presented new measurements and simulations of air δ 15 N profiles from several Antarctic ice cores spanning the last deglaciation. First, our new δ 15 N measurements highlight a heterogeneous behaviour of the firn structure evolution over the last deglaciation. In particular, TALDICE glacial δ 15 N-data values are similar to interglacial δ 15 N-data values and this undermines the hypothesis for a significant impact of the snow impurity content on the firn structure for this site. At BI we measure strong millennial-scale δ 15 N-data variations during climatic intervals associated with relatively flat δD, revealing that processes independent from the water isotopes affect the firn structure at this site. Moreover, our new results enable us to rule out the hypothesis of a large glacial convective zone as the single explanation for the model-data δ 15 N mismatch observed at EDML. Still, direct constraints on the extent of past convective zone, especially during the LGM (e.g. Severinghaus et al., 2006), are necessary not only in order to strengthen confidence in our conclusion for the EDML case but also to assess this hypothesis for other Antarctic sites.
Second, our δ 15 N model and data syntheses show that complex competition between the opposite impacts of changes in surface temperature and accumulation rate is at play during the last deglaciation in Antarctic firn. We suggest that the role of temperature in firnification process may have been overestimated in past studies, while the role of accumulation rate should be revised in current firn models. These new results also highlight the importance of using accurate past surface accumulation rate estimates to force firnification models. The processes that could induce deviations from simple relationships between accumulation, temperature and precipitation isotopic composition require more in-depth studies (e.g. Sime et al., 2013), and future highresolution chemical tracer profiles should help constraining past changes of the accumulation rate independently from water isotopic profiles, especially for coastal sites such as BI.
Overall, the temporal evolution of the firn structure is likely to result from a site-specific complex interaction between several of the discussed mechanisms, explaining why current firn densification models do not correctly resolve all the processes controlling δ 15 N variations. Further studies are necessary to be able to separate the effects of surface temperature, accumulation and impurity content on firn densification. New firn air sampling and ice cores drilled in West Antarctica (Fletcher Promontory, WAIS) will allow future investigations on the current and past firn structure in coastal and semi-coastal regions.

Establishing a chronology for the BI and JRI ice cores
Glaciological chronologies have been derived for BI and JRI following an approach similar to that presented in Parrenin et al. (2007b). It consists of an accumulation model and an ice flow model. The accumulation is assumed to be exponentially related to the isotopic content of the ice following Eq. (4) given in the main manuscript. The ice flow model is a simplified pseudo-steady state model (Parrenin and Hindmarsh, 2007) -that is to say that the geometry (bedrock and surface elevation) and the ratio melting/accumulation are assumed constant in time. For each ice core, the free parameters in the model including A 0 and β from Eq. (4) have been adjusted so that the resulting timescale is in good agreement with age markers obtained by comparison of ice and gas records to other well-dated palaeorecords. Stratigraphic markers were derived from matching gas records (δ 18 O atm , CH 4 and CO 2 , unpublished data) with gas records from EPICA Dome C (EDC), Byrd and Vostok on the EDC3 chronology (Le Floch et al., 2007;Parrenin et al., 2007a;F. Parrenin, personal communication, 2013). The Goujon model has then been forced with surface temperature and accumulation scenarios to estimate age, allowing production of the gas age scale.
It is beyond the scope of this paper to discuss in details the preliminary JRI and BI age scales. Sensitivity tests have shown that the underlying dating uncertainties do not affect the results discussed in the main manuscript.

Sensitivity tests with the Goujon Model (2003)
Sensitivity tests have been run for the EDML, TALDICE, EDC and BI sites with the Goujon firn densification model to examine the respective effects of surface temperature and accumulation rate variations over the last deglaciation and their influence on the LID and thus the evolution of MODELδ 15 N.
First, for each site we have compared two simulated curves (Fig. B1, left panels): -An "Acc MODEL-δ 15 N" curve representing MODELδ 15 N simulated in response to accumulation changes only. For that purpose the Goujon model is forced with a surface temperature scenario fixed to the present surface temperature (given on Fig. 1) and the original scenario of accumulation rate deduced from Eq. (4) with the β value given in Table 1.
-A "Temp MODEL-δ 15 N" curve representing MODELδ 15 N simulated in response to temperature changes only. For that purpose, the Goujon model is forced with an accumulation rate scenario fixed to the present surface accumulation rate (given on Fig. 1) and the original scenario of surface temperature deduced from Eqs. (2) and (3) with the α value given in Table 1.
For the four sites, we systematically observe the opposite influence of surface temperature and accumulation rate on firnification processes. Over the deglaciation the effect of an increase of accumulation rate only leads to larger δ 15 N values during the EH than during the LGM, while the effect of an increase in surface temperatures leads to smaller δ 15 N values during the EH than during the LGM. The total MODEL-δ 15 N curve is not simply the average of the two δ 15 N simulations considering each single factor. Non-linear interactions occur as the accumulation rate influence is different for different temperature levels, and vice versa. Second, to better investigate the complex interaction between the two, we have performed a second sensitivity test. For each site we have run the Goujon model forced by inputs parameters deduced from water isotopic profiles but with slightly different values for the coefficients α and β to convert them into the past surface temperature and accumulation rate through Eqs. (3) and (4) The MODEL-δ 15 N pattern observed at the start of the deglaciation (first half of phase 1) can be explained by the corresponding increases in accumulation rate at all sites (Appendix B; Fig. B1). However, depending on the warming amplitude at the start of the deglaciation and the associated increase in accumulation rate, our sensitivity tests suggest that during the second half of phase, either the MODEL-δ 15 N increases (i.e. all tested scenarios for BI, Fig. A1, right panels), or the MODEL-δ 15 N decrease (i.e. all tested scenarios for EDML and Acc high MODEL-δ 15 N curve for EDC, Fig. B1, right panels). These results illustrate again the complex interaction between the effect of surface temperature and the effect of accumulation rate on the LID.
For the four scenarios tested on the four sites, we systematically observe that the Goujon model predicts a decrease of MODEL-δ 15 N during the Antarctic Cold Reversal, followed by an increase of MODEL-δ 15 N toward the early Holocene (Fig. B1, right panel). It illustrates that a warming phase of smaller amplitude than a glacial-interglacial transition could lead to deepening in the firn.

Constraints on the EDML depth
To derive depth estimated from DATA-δ 15 N, we deduce the diffusive column height (DCH) from δ 15 N using the barometric equation (Eq. 1) and we convert then the DCH to depth through the following equation: In this equation we account for the thinning (t) due to ice flow by multiplying by the appropriate thinning factor from an ice flow model. The coefficient 0.7 represents the ratio of column-averaged firn density to ice density, which is required to translate the firn equivalent DCH thickness to an ice equivalent DCH thickness. We adopt a 5 % uncertainty to account for variations with respect to firn density profiles as a function of temperature and accumulation rate and varying ice density (Blunier et al., 2004). We also used two different estimates of the thinning factor: one from the EDML glaciological model of Huybrechts et al., 2007) and one from the new AICC2012 chronology (Bazin et al., 2012), and we consider a 10 % uncertainty linked to this parameter. Our new DATA-δ 15 N translate into a depth of 26.8 ± 3 m for the 1360-1400 m depth interval. To deduce depth from the Goujon model, we used Eq. (4) with the LID estimate from the Goujon model , taking into account similar uncertainties as previously described above.
Finally, we have revised the two depth estimates that Loulergue et al. (2007) deduced over the Laschamp event (41.2 ± 1.6 ka on GICC05; Svensson et al., 2007) using independent matching of the gas (CH 4 ) and ice ( 10 Be) records of the EDML ice core with the NorthGRIP gas and ice records. They estimate two empirical depth of 21.4 ± 4.6 m and 23.1 ± 4.6 m at about 1368-1407 m depth. Most of the estimated uncertainty is due to the complicated identification of tie points between the EDML CH 4 record and the corresponding NorthGRIP δ 18 O ice record (see Table 1 and Fig. 3b of their paper for more details). In order to reduce the uncertainty of the Loulergue et al. (2007) empirical depth estimates, we use an objective method based on the match protocol (Lisiecki and Lisiecki, 2002) to give the best matching between EDML and NorthGRIP CH 4 records. Our revised match between EDML and NorthGRIP CH 4 records is displayed on Fig. C2. After identifying the CH 4 change concomitant with the 10 Be peak on the GICC05 timescale, the corresponding depth (following the notation of Loulergue et al., 2007) can directly be read on the CH 4 record on the EDML depth scale from the correspondence between NorthGRIP and EDML CH 4 ). We thus obtain two EDML depths corresponding to the depths of CH 4 changes concomitant with the age of the two 10 Be peaks: 1390.4 m (instead of 1389.8 m in Loulergue et al., 2007) and 1408.5 m (instead of 1406.4 m in Loulergue et al., 2007). The two revised empirical depth deduced from our approach are slightly larger than the original estimates by Loulergue et al. (2007): 22.0 ± 2.5 m and 25.2 ± 2.5 m instead of 21.4 ± 4.6 m and 23.1 ± 4.6 m, respectively. The uncertainty is linked to the mean resolution of the EDML CH 4 record (2.5 m) and to the rate of CH 4 change (Fig. C1). We see that the CH 4 record undergoes a clear minimum at 1408.5 m corresponding to the second 10 Be peak for which we estimate the depth of 25.2 ± 2.5 m. Therefore, we consider this estimate as the most robust. These estimates are larger than the depth deduced from EDML firn modelling based on two different accumulation rate histories (modelled depth of about 22.9 m for Scenario 1 and 21.2 m for Scenario 4, as defined by Loulergue et al., 2007).

E. Capron et al.: Glacial-interglacial dynamics of Antarctic firn columns
Our results show that the modelled depth and the δ 15 Nbased depth estimates are compatible in general within the considered uncertainty range; however, we cannot draw a firm conclusion on the presence and thickness of the convective zone at EDML during the LGM. Still, we can provide a more precise evaluation of the "deep convective zone" hypothesis over the Laschamp event by comparing the two empirical depth estimates with the depth estimate based on δ 15 N measurements equal to 26.8 ± 3 m over the 1360-1400 m depth interval. This value is slightly larger than the empirical depth estimates from gas CH 4 and ice 10 Be matching. The existence of a convective zone would lead to a δ 15 N-based depth smaller than the modelled or empirically derived depth, which is the opposite of our observation. Thus, our new results show the absence of a convective zone at EDML at the time of the Laschamp event.