The WAIS Divide deep ice core WD2014 chronology - Part 1

. The West Antarctic Ice Sheet Divide (WAIS Divide, WD) ice core is a newly drilled, high-accumulation deep ice core that provides Antarctic climate records of the past ∼ 68 ka at unprecedented temporal resolution. The upper 2850 m (back to 31.2 ka BP) have been dated using annual-layer counting. Here we present a chronology for the deep part of the core (67.8–31.2 ka BP), which is based on stratigraphic matching to annual-layer-counted Greenland ice cores using globally well-mixed atmospheric methane. We calculate the WD gas age–ice age difference ( (cid:49) age) using a combination of ﬁrn densiﬁcation modeling, ice-ﬂow modeling, and a data set of δ 15 N-N 2 , a proxy for past ﬁrn column thickness. The largest (cid:49) age at WD occurs during the Last Glacial Maximum, and is 525 ± 120 years. Internally consistent solutions can be found only when assuming little to no inﬂuence of impurity content on densiﬁcation rates, contrary to a recently proposed hypothesis. We synchronize the WD chronology to a linearly scaled version of the layer-counted Greenland Ice Core Chronology (GICC05), which brings the age of Dansgaard–Oeschger (DO) events into agreement with the U / Th absolutely dated Hulu Cave speleothem record. The small (cid:49) age at WD provides valuable opportunities to investigate the timing of atmospheric greenhouse gas variations relative to Antarctic climate, as well as the interhemispheric phasing of the “bipolar seesaw”.


154
C. Buizert et al.: The WAIS Divide chronology -Part 1 (Gow et al., 1997), and geothermal heat flow (Dahl-Jensen et al., 1998). Having a reliable ice core chronology (i.e., an age-depth relationship) is paramount for the interpretation of the climate records and comparison to marine and terrestrial paleoclimate archives.
The West Antarctic Ice Sheet Divide (WAIS Divide, WD) ice core (79.48 • S, 112.11 • W; 1766 m above sea level; −30 • C present-day mean annual temperature) was drilled and recovered to 3404 m depth (WAIS Divide Project Members, 2013). Drilling was stopped 50 m above the estimated bedrock depth to prevent contamination of the basal hydrology. Due to high accumulation rates of 22 cm ice a −1 at present and ∼ 10 cm ice a −1 during the Last Glacial Maximum (LGM), the WD core delivers climate records of unprecedented temporal resolution (Steig et al., 2013;Sigl et al., 2013) as well as gas records that are only minimally affected by diffusive smoothing in the firn column (Mischler et al., 2009;Mitchell et al., 2011Marcott et al., 2014). The combination of high accumulation rates and basal melting at the WD site results in ice near the bed that is relatively young (∼ 68 ka) compared to cores drilled in central East Antarctica.
In WD, annual layers can be identified reliably for the upper 2850 m of the core, reaching back to 31.2 ka BP (thousands of years before present, with present defined as 1950 CE). Below 2850 m depth an alternative dating strategy is needed. Several methods have been employed previously at other deep ice core sites. First, orbital tuning via δO 2 / N 2 has been applied successfully to several Antarctic cores (Bender, 2002;Kawamura et al., 2007). However, an age span of only ∼ 3 precessional cycles in WD, in combination with the low signal-to-noise ratio of δO 2 / N 2 data, makes this technique unsuitable for WD. The uncertainty in the orbital tuning is about one-fourth of a precessional cycle (∼ 5 ka), making it a relatively low-resolution dating tool. Second, in Greenland, ice-flow modeling has been used to extend layer-counted chronologies (e.g., Johnsen et al., 2001;Wolff et al., 2010). This method requires assumptions about past accumulation rates, ice flow, and ice sheet elevation. Particularly for the oldest WD ice, the resulting uncertainty would be substantial. Third, several radiometric techniques have been proposed to date ancient ice. Radiocarbon ( 14 C) dating of atmospheric CO 2 trapped in the ice is unsuitable as it suffers from in situ cosmogenic production in the firn (Lal et al., 1990), and the oldest WD ice dates beyond the reach of 14 C dating. Other absolute (radiometric) dating techniques, such as recoil 234 U dating (Aciego et al., 2011), 81 Kr dating (Buizert et al., 2014a), or atmospheric 40 Ar buildup (Bender et al., 2008), currently suffer from uncertainties that are too large (≥ 20 ka) to make them applicable at WD.
Instead, at WD we use stratigraphic matching to welldated Greenland ice cores using globally well-mixed atmospheric methane (CH 4 ) mixing ratios (Blunier et al., 1998;Blunier and Brook, 2001;Blunier et al., 2007;Petrenko et al., 2006;EPICA Community Members, 2006;Capron et al., 2010). This method is particularly suited to WD because of the small gas age-ice age difference ( age, Sect. 3) and the high-resolution, high-precision CH 4 record available (Sect. 2.1). The method has three main sources of uncertainty: (i) the age uncertainty in the records one synchronizes to, (ii) age of the ice core being dated, and (iii) the interpolation scheme used in between the CH 4 tie points. We present several improvements over previous work that reduce and quantify these uncertainties: (i) we combine the layercounted Greenland Ice Core Chronology (GICC05) and a recently refined version of the U / Th-dated Hulu speleothem record (Edwards et al., 2015;Reimer et al., 2013;Southon et al., 2012) to obtain a more accurate estimate of the (absolute) ages of abrupt ; (ii) we combine firn densification modeling, iceflow modeling, a new WD δ 15 N-N 2 data set that spans the entire core, and a Monte Carlo sensitivity study to obtain a reliable age estimate (Sect. 3); and (iii) we compare four different interpolation schemes to obtain an objective estimate of the interpolation uncertainty (Sect. 4.5).
This work is the first part in a series of two papers describing the WD2014 chronology for the WD core in detail. The second part describes the development of the annual layer count from both multi-parameter chemistry and electrical conductivity measurements. The WD2014 chronology is currently the recommended gas and ice timescale for the WD deep core, and as such it supersedes the previously published WDC06A-7 chronology (WAIS Divide Project Members, 2013).

Data description
Measurements of water stable isotopes. Water isotopic composition (δ 18 O and δD = δ 2 H) was measured at IsoLab, University of Washington. Procedures for the deep section of the core are identical to those used for the upper part of the core reported in WAIS Divide Project Members (2013) and Steig et al. (2013). Measurements were made at 0.25 to 0.5 m depth resolution using laser spectroscopy (Picarro L2120-i water isotope analyzer), and normalized to VSMOW-SLAP (Vienna Standard Mean Ocean Water -Standard Light Antarctic Precipitation). The precision of the measurements is better than 0.1 and 0.8 ‰ for δ 18 O and δD, respectively.
Measurements of CH 4 . Two CH 4 data sets were used for WD. The first is from discrete ice samples, and was measured jointly at Pennsylvania State University (0-68 ka, 0.5-2 m resolution) and Oregon State University (11.4-24.8 ka, 1-2 m resolution). Air was extracted from ∼ 50 g ice samples using a melt-refreeze technique, and analyzed on a standard gas chromatograph equipped with a flame-ionization detector. Corrections for solubility, blank size and gravitational enrichment are applied (Mitchell et al., 2011;WAIS Divide Project Members, 2013). The second data set is a continuous CH 4 record measured by coupling a laser spectrometer to a continuous flow analysis setup (Stowasser et al., 2012;Rhodes et al., 2013;Chappellaz et al., 2013), and was measured jointly by Oregon State University and the Desert Research Institute (Rhodes et al., 2015). The continuous data set is used to identify the abrupt DO transitions, as it provides better temporal resolution and analytical precision. Both records are reported on the NOAA04 scale (Dlugokencky et al., 2005). Analytical precision in the CH 4 data (2σ pooled standard deviation) is around 3.2 and 14 ppb for the discrete data from Oregon State University and Pennsylvania State University, respectively, and 3 to 8 ppb for the continuous CH 4 data, depending on the analyzer used (Rhodes et al., 2015); the 14 ppb stated for the PSU discrete data may be an overestimation, as depth-adjacent (rather than true replicate) samples were used in the analysis.
Measurements of δ 15 N. Atmospheric N 2 isotopic composition (δ 15 N) was measured at Scripps Institution of Oceanography, University of California. Air was extracted from ∼ 12 gram ice samples using a melt-refreeze technique, and collected in stainless steel tubes at liquid-He temperature. δ 15 N was analyzed using conventional dual-inlet isotope ratio mass spectrometry (IRMS) on a Thermo Finnigan Delta V mass spectrometer. Results are normalized to La Jolla (California, USA) air, and routine analytical corrections are applied (Sowers et al., 1989;Petrenko et al., 2006;Severinghaus et al., 2009). Duplicates were not run for most δ 15 N data in this study, but the pooled standard deviations of Holocene WD δ 15 N data sets with duplicate analyses are 0.003 ‰ (Orsi, 2013). We conservatively adopt an analytical uncertainty of 0.005 ‰ for this data set to allow for other sources of error.
Measurements of [Ca]. Ca concentrations in the ice were measured at the Ultra Trace Chemistry Laboratory at the Desert Research Institute via continuous flow analysis. Longitudinal samples of ice (approximately 100 cm × 3.3 cm × 3.3 cm) were melted continuously on a melter head that divides the meltwater into three parallel streams. Elemental measurements were made on meltwater from the innermost part of the core with ultra-pure nitric acid added to the melt stream immediately after the melter head; potentially contaminated water from the outer part of the ice is discarded. Elemental analysis of the innermost meltwater stream is performed in parallel on two inductively coupled plasma mass spectrometers (ICPMS), each measuring a different set of elements; some elements were analyzed on both. The dual ICPMS setup allows for measurement of a broad range of 30 elements and data quality control (McConnell et al., 2002(McConnell et al., , 2007. Precision of the Ca measurements in WD glacial ice is estimated to be ±3 %, with a lower detection limit of 0.15 ng g −1 . Continuous Ca and CH 4 measurements are done on the same ice, and are exactly co-registered in depth.

Firn densification model description
Air exchange with the overlying atmosphere keeps the interstitial air in the porous firn layer younger than the surrounding ice matrix, resulting in an age difference between polar ice and the gas bubbles it contains, commonly referred to as age (Schwander and Stauffer, 1984). Here we use a coupled firn-densification-heat-diffusion model to calculate age back in time (Barnola et al., 1991;Goujon et al., 2003;Schwander et al., 1997;Rasmussen et al., 2013), constrained by measurements of δ 15 N of N 2 , a proxy for past firn column thickness (Sowers et al., 1992). The model is based on a dynamical description of the Herron-Langway model formulated in terms of overburden load (Herron and Langway, 1980), which is solved in a Lagrangian reference frame. This model has been applied previously to the Greenland NEEM, NGRIP, and GISP2 cores Seierstad et al., 2015;Buizert et al., 2014b), where it gives a good agreement to the Goujon densification model Goujon et al., 2003). The model allows for the inclusion of softening of firn in response to impurity loading (Horhold et al., 2012), following the mathematical description of Freitag et al. (2013a). The equations governing the model densification rates are given in Appendix A.
The model uses a 2-year time step and 0.5 m depth resolution down to 1000 m, the lower model boundary. A thick model domain is needed because of the long thermal memory of the ice sheet. At WD, downward advection of cold surface ice is strong due to the relatively high accumulation rates, and the geothermal gradient does not penetrate the firn column (Cuffey and Paterson, 2010). We further use a lockin density that equals the mean close-off density (Martinerie et al., 1994) minus 17.5 kg m −3 (as in Blunier and Schwander, 2000) and an empirical parameterization of lock-in gas age based on firn air measurements from 10 sites .
We furthermore use the steady-state version of the Herron-Langway model (Herron and Langway, 1980) in performing sensitivity studies (Sect. 3.2) and the dynamical Arnaud model (Arnaud et al., 2000;Goujon et al., 2003) to validate our age solution.

Temperature reconstruction and ice-flow model
Our temperature reconstruction (Fig. 1a) is based on water δD, a proxy for local vapor condensation temperature, calibrated using a measured borehole temperature profile (following Cuffey et al., 1995;Cuffey and Clow, 1997) and, for the last 31.2 ka, adjusted iteratively to satisfy constraints on firn thickness provided by δ 15 N and by the observed layer thickness λ(z). Using δ 18 O rather than δD in the temperature reconstruction leads to differences that are negligibly small. This borehole temperature calibration approach is possible at WD because the large ice thickness and relatively high accumulation rates help to preserve a memory of past  (c) δ 15 N data (black dots) with densification model output (green). (d) age calculated using the densification model (orange) and using the Parrenin depth method (black) with constant 4 m thick convective zone and no correction for thermal δ 15 N fractionation. (e) Modeled thinning function from ice-flow model (solid), and a simple Nye strain model for comparison (dashed); the Nye thinning function, which has a uniform strain rate as a function of depth, is given as f λ (z) = (H − z)/H with H the ice sheet thickness (Cuffey and Paterson, 2010, p. 616). temperatures in the ice sheet. A coupled 1-D ice-flow-heatdiffusion model converts surface T (t) into a depth profile for comparison to measured borehole temperatures. The 1-D ice-flow model calculates the vertical ice motion, taking into account the surface snow accumulation, the variation in density with depth, and a prescribed history of ice thickness. Vertical motion is calculated by integrating a depth profile of strain rate and adding a rate of basal melt. As in the model of Dansgaard and Johnsen (1969), the strain rate maintains a uniform value between the surface and a depth equal to 80 % of the ice thickness, and then varies linearly to some value at the base of the ice. This basal value is defined by the "basal stretching parameter" f b , the ratio of strain rate at the base to strain rate in the upper 80 % of the ice column. The basal ice is melting, so part of the ice motion likely occurs as sliding. The along-flow gradient in such sliding is unknown and thus so too is the parameter f b . We overcome this problem by making both the current ice thickness and the basal melt rate free parameters when optimizing models with respect to measured borehole temperatures. Because the basal melt rate and f b affect the vertical velocities in similar fashion, the optimization constrains a combination of melt rate and f b that is tightly constrained by the measured temperatures. Thus we find that varying f b through a large range, from 0.1 to 1.5, changes the reconstructed LGM temperature by less than 0.2 • C. Effects of the prescribed ice-thickness history are likewise minor; assuming a 150 m thickness increase from the LGM to 15 ka changes the reconstructed LGM temperature by less than 0.2 • C compared to a constant thickness. Note that the 1-D flow model used here is simpler than the one used by Cuffey and Clow (1997) in that it does not attempt to calculate changes in the shape of the strain rate profile; the unknown basal sliding motion at the WD site negates the usefulness of such an exercise.
One For the orange curve we have used the Nye thinning function (Fig. 1e); the final ice-flow model optimizes the agreement between A(t) obtained from λ(z) and from the inverse firn modeling approach. The WD2014 chronology uses the orange A init (t) scenario. (b) A(t) found in the inverse firn modeling approach using both A init (t) scenarios; color coding as in panel (a). The function ξ(t) is found as follows. We use control points at 1500-year intervals (blue dots); the algorithm has the freedom to change the value of ξ(t) at each of these points. In between the control points, ξ(t) is found via linear interpolation. (c) Modeled age using both A init (t) scenarios; color coding as in panel (a).
strain is represented by the thinning function f λ (z) (Cuffey and Paterson, 2010), the ratio of annual-layer thickness at depth in the ice sheet to its original ice-equivalent thickness at the surface when deposited. The modeled thinning function is shown in Fig. 1e (solid line). In the deep part of the ice sheet, f λ (z) becomes increasingly uncertain as the unknown basal melt rate and f b become the dominant controls. Here we optimize the model by comparing accumulation rates derived from f λ (z) with those implied by a firn densification model and the measured δ 15 N of N 2 (Sect. 3.1). While this has little effect on the temperature history reconstruction, it provides an important constraint on calculated basal melt rate, an interesting quantity for ice dynamics studies. Our analysis of basal melt rates and further details of the temperature optimization process and 1-D flow modeling will be presented elsewhere.
3 The gas age-ice age difference ( age)

The WD2014 age reconstruction
The firn densification forward model uses past surface temperature T (t) and accumulation A(t) as model forcings, and provides age(t) and δ 15 N(t) as model output.
For the past 31.2 ka, WD has an annual-layer-counted chronology; for this period the annual-layer thickness λ(z) provides a constraint on past accumulation rates via λ(z) = A(z) × f λ (z). WD accumulation reconstructed from λ(z) is plotted in black in Fig. 1b.
Prior to 31.2 ka we have no such constraint on A(t), and an alternative approach is needed. We use the densification model as an inverse model, where we ask the model to find the A(t) history that minimizes the root-mean-square (rms) deviation between measured and modeled δ 15 N, given the T (t) forcing. The δ 15 N data and model fit are shown in Fig. 1c, the A(t) history that optimizes the δ 15 N fit is shown in Fig. 1b (red), and the modeled age is shown in Fig. 1c (orange). The optimal A(t) history is estimated in two steps. First, we make an initial estimate A init (t) for the past accumulation history. Second, we adjust the A(t) forcing by applying a smooth perturbation ξ(t) such that ; an automated algorithm is used to find the curve ξ(t) that optimizes the model fit to the δ 15 N data. For the last 31.2 ka we obtain a good agreement between A obtained from λ(z) and the modeled f λ (z) (Fig. 1b, black) and A obtained from the inverse method (red). The solution we present here is therefore fully internally consistent, i.e., the A and T histories used in the firn densification modeling are the same as those used in the ice-flow modeling, and they provide a good fit to both the δ 15 N data and borehole temperature data. WD does not suffer from the δ 15 N modeldata mismatch that is commonly observed for East Antarctic cores during the glacial period (Landais et al., 2006;Capron et al., 2013).
We base our A init values on λ(z) for the past 31.2 ka; prior to that we use the common assumption that A follows δ 18 O (i.e., Clausius-Clapeyron scaling); the fit to the δ 15 N data is optimized for To test the validity of the Clausius-Clapeyron assumption, we additionally run the scenario A init (t) = 0.22 m a −1 (i.e., constant accumulation at present-day level). The A(t) and age reconstructed under both A init scenarios are similar at multimillennial timescales (Fig. 2). In the layer-counted interval (< 31.2 ka BP), A obtained from λ(z) and δ 18 O is significantly coherent at all timescales longer than 3000 years, but not at higher frequencies. This is equivalent to the variability resolved in the A init (t) = 0.22 m a −1 scenario above. We conclude that the WD δ 15 N data support the idea that A follows δ 18 O on multi-millennial timescales. However, there may not be a strong relationship at timescales less than a few thousand years, as is clear from the abrupt A increase around 12 ka seen in λ(z) that is not reflected in δ 18 O ( Fig. 1a and b). For consistency between the upper and deeper part of the core we use the age values obtained with the inverse densification model for the entire core.
Recently, another δ 15 N-based approach has been suggested that uses depth, rather than age, in reconstructing gas chronologies (Parrenin et al., 2012). This method removes the dependence on T (t) and replaces this with a dependence on the thinning function f λ (z). Note that this method is very successful in the upper part of an ice core, where f λ (z) is well constrained, but not very reliable near the base, where f λ (z) is highly uncertain. Therefore, the firn densification modeling approach should be considered to be more reliable at WD during marine isotope stages (MIS) 2 through 4. Results from the depth method are plotted in black in Fig. 1c, and generally show good agreement with the firn modeling approach. A notable exception is the 60-65 ka interval, where the depth method overestimates the age due to the fact that we have to compress λ(z) strongly in order to fit age constraints derived from DO 18 (Sect. 4.5).
Last, we want to point out that the δ 15 N data support an early warming at WD, as reported recently (WAIS Divide Project Members, 2013). WD δ 15 N starts to decrease around 20.5 ka BP, suggesting a thinning of the firn column. The λ(z) (as derived from the layer count) shows that accumulation did not change until 18 ka BP, at which point it started to increase (which would act to increase the firn thickness). The most plausible explanation for the δ 15 N decrease around 20.5 ka BP is therefore an early onset of West Antarctic deglacial warming, in agreement with increasing δ 18 O around that time. The warming enhances the densification rate of polar firn, thereby decreasing its thickness (e.g., Herron and Langway, 1980).

age sensitivity study
Besides A and T there are several model parameters that have the potential to influence the model outcome; these are the convective zone (CZ) thickness (Sowers et al., 1992;Kawamura et al., 2006), surface density (ρ 0 ), and sensitivity to ice impurity content. In this section we evaluate the sensitivity of the model output to all of these parameters. We performed 1000 model runs in which the model parameters were randomly perturbed. The spread in age model results is used to calculate the WD2014 age uncertainty. Convective zone thickness. In the WD2014 model run (Sect. 3.1) we use a constant 3.5 m CZ, corresponding to the present-day situation (Battle et al., 2011). In the sensitivity study we vary the CZ by one of two methods: (1) we let the CZ be constant in time; its thickness is set by drawing from a Gaussian distribution with 3.5 m mean and 3.5 m 2σ width (i.e., 95 % probability of drawing a value in the 0-7 m range).
(2) We let the CZ be a function of accumulation rate (Dreyfus et al., 2010), CZ = 3.5 + k × (A-0.22); we draw k from a Gaussian distribution with mean of −10 and a 2σ width of 40 (at an LGM A of 10 cm a −1 this gives a CZ of 0-10 m thickness). In both methods, whenever CZ values are selected that are smaller than 0 m, the CZ thickness is set to 0 m instead. For each of the 1000 model runs in the sensitivity study we randomly selected either of the two methods.
Surface density. In the WD2014 model run we use past surface densities (ρ 0 ) as given by the parameterization of Kaspers et al. (2004). In the sensitivity study we add a constant offset to the Kaspers values, the magnitude of which is drawn from a Gaussian distribution of zero mean and a 2σ width of 60 kg m −3 . This range corresponds to the full range of observed ρ 0 variability in Kaspers et al. (2004).
Past temperatures. Model temperature forcing is constrained by δD and measured borehole temperatures. There is, however, a range to the solutions allowed by the borehole temperature and ice-flow model; here we use the upper and lower extremes of this range, determined by Monte Carlo analysis using uncertainties of input variables. The scenarios were chosen to provide the maximum T range for the glacial period rather than for the Holocene, because we are interested in the uncertainty in the methane synchronization (68-31.2 ka BP). In the sensitivity study we use T (t) = T optimal (t) + κ × T (t), where T optimal is the forcing used in the WD2014 model run (Fig. 1a), T (t) is half the difference between the maximum-T and minimum-T scenarios, and κ is drawn from a Gaussian distribution of zero mean and unit 2σ width (giving 95 % probability that T (t) is within the extreme range identified from the borehole, Fig. 3a).
δ 15 N uncertainty. We conservatively adopt an analytical uncertainty of 0.005 ‰ for this data set; in addition, the interpretation of δ 15 N in terms of firn thickness is subject to further uncertainty due to irregular firn layering and the stochastic nature of bubble trapping, as was observed for other atmospheric gases such as CH 4 (Etheridge et al., 1992;Rhodes et al., 2013). For each run of the sensitivity study, we therefore perturb each of the individual δ 15 N data points by adding an offset that is drawn from a Gaussian distribution of zero mean and a 2σ width of 0.015 ‰.
Impurity-enhanced densification. Following recent work we include the possibility that increased glacial impurity loading could have enhanced densification rates (Horhold et al., 2012;Freitag et al., 2013a). We use the mathematical formulation of Freitag et al. (2013a), in which the activation energy of the sintering process is a function of the Ca concentration in the firn. The value of β, the sensitivity to Ca, is drawn from a Gaussian distribution with 0.0015 mean and a 2σ width of 0.0015. The topic of impurity-enhanced densification is discussed in detail in Sect. 3.3.
The A and age scenarios found in the sensitivity study are shown in Fig. 3b and c, respectively. The shaded areas in Fig. 3b and c give the total range of solutions, as well as the ±2σ and ±1σ confidence intervals. Note that the total range of solutions will depend on the number of model runs (here 1000) but that the position of the ±2σ and ±1σ envelopes will not. To investigate the distribution of values, we include histograms of age at 20 kyr intervals ( Fig. 3d-f). Based on the sensitivity study, we estimate the WD age to be 521 ± 120 years (2σ uncertainty) at the LGM (∼ 20 ka BP). The age value of 351 ± 73 years at 40 ka BP gives a representative age for MIS 3; age at 60 ka BP is 262 ± 50 years. Additionally, we have repeated our age reconstruction using the firn densification physics described by Arnaud et al. (2000) rather than the Herron-Langway description used so far; the Arnaud model provides the physical basis for the commonly used firn densification model of Goujon et al. (2003). More details on the implementation of the Arnaud model are given in Appendix A. age found using the Arnaud model is plotted in red in Fig. 3c. Averaged over the entire core, age found with the Arnaud model is 19 years (about 7 %) smaller than age from the Herron-Langway model. The root-mean-square (rms) difference between both solutions is 35 years, corresponding to 0.63 times the 2σ uncertainty found in the sensitivity study. Both solutions are thus found to be in good agreement. The Herron-Langway approach is preferred because the internally consistent solution of temperature, accumulation, and ice flow associated with it is in better agreement with borehole temperature data than the solutions associated with the Arnaud model. Furthermore, the Herron-Langway model is more successful in simulating the magnitude of the δ 15 N response to the accumulation anomaly at 12 ka (not shown), suggesting it has a more realistic sensitivity to accumulation variability. We use Ca sensitivities β = 0 (red) through β = 1 × 10 −2 (blue), in steps of 2.5 × 10 −3 (shades of deep purple). Black curves give A and age from ice-flow modeling and λ(z).

Impurity softening of firn?
Recent work suggests a link between densification rates and impurity content (for which [Ca 2+ ] is used as a proxy) in polar firn (Horhold et al., 2012;Freitag et al., 2013a). Here we measured total [Ca] by ICP-MS, but at WD nearly all Ca is in the form of Ca 2+ . The influence of the impurity sensitivity β (see Eq. A6 in the Appendix) on age at WD is shown in Fig. 4. The sensitivity recommended by Freitag et al. (2013a) from investigating present-day firn packs is β = 1×10 −2 . We reconstructed A and age with the firn densification inverse model using five values of β ranging from β = 0 (red) to β = 1 × 10 −2 (blue) in steps of 2.5 × 10 −3 . Average [Ca] is around 0.8 ng g −1 in the early Holocene and around 9 ng g −1 in the LGM; a change to about an order of magnitude. Following Freitag et al. (2013a) we use the total [Ca] rather than non-sea-salt Ca. If densification rates are sensitive to impurity loading (large β, blue curves), this results in increased firn compaction during the LGM. The densification model, which is trying to match the δ 15 N data, will compensate by increasing the A forcing, which in turn results in a decreas-ing age. Hence the model simulations with large β (blue) give a higher A and smaller age.
For the past 31.2 ka we have an independent A estimate from λ(z) that we can compare to the solutions from the firn model (Fig. 4, black curve). We also plotted age reconstructed via the depth method of Parrenin et al. (2012). Remarkably, we find consistent solutions only when using a Ca sensitivity β ≤ 2.5 × 10 −3 , i.e., less than one-quarter of the sensitivity suggested by Freitag et al. (2013a). The best fit to the independent LGM (25-20 ka BP) A and age estimates is obtained for β = 0. We conclude that WD does not provide any evidence for impurity (or, more specifically, Ca) enhancement of densification rates.
An important caveat is that our model uses 10-yearaverage [Ca] values, and therefore cannot resolve effects of interannual layering within the firn. Explicitly modeling the layering would require centimeter-scale resolution in the dynamical firn model, which is prohibitive from a computational point of view. Furthermore, [Ca] data at the required sub-annual resolution are difficult, if not impossible, to measure for the deepest part of the core, where λ(z) is below 1 cm a −1 . Increased firn layering and enhanced bulk densification affect the firn thickness in a similar manner; both lead to a shallower lock-in depth, and thereby a reduced δ 15 N. Therefore, in order to reconcile our WD results with the impurity hypothesis of Horhold et al. (2012), one would need to invoke a strong reduction in LGM firn layering relative to the present day to compensate for the impurity-driven increase in bulk densification rates. Recent work on the EDML core suggests that firn density layering may have been more pronounced during glacial times (Bendel et al., 2013); including firn layering is therefore likely to only exacerbate the problem.
Work on present-day firn has provided support for firn softening by impurity loading (Horhold et al., 2012;Freitag et al., 2013a, b). More work is needed to understand how densification rates are linked to impurity content in a mechanistic, rather than purely empirical, way. Perhaps such a microscopic description could provide an explanation why firn densification rates at WD, to first order, do not appear to be affected by order-of-magnitude variations in [Ca] loading. One possible explanation could be that densification rates are controlled by some parameter that co-varies with Ca in modern day firn yet does not change appreciably over glacial cycles (Fujita et al., 2014). 3. a reassessment of the dating between 2300 and 2800 m (15.3-29.5 ka) using ECM and dust particle measurements, with the ECM having increasing importance with depth; 4. an extension of the annual-layer dating between 2800 and 2850 m (29.5-31.2 ka) using ECM.
Details on the updated WD layer count and the layer counting methodology are presented in part 2 of the WD2014 papers.

Methane synchronization (31.2-68 ka)
For the deep part of the core where an annual-layer count is not available, we date WD by synchronization to well-dated Northern Hemisphere (NH) climate records of abrupt DO variability using the WD record of globally well-mixed CH 4 (Fig. 5). This process consists of several steps: 6. Repeat steps 3-5 iteratively until the depth-age relationship is stable within 1 year. At WD this happened after three iterations.
These steps are described in more detail in the following sections.

Establishing the midpoint in abrupt DO transitions
The procedure for determining the midpoint of the abrupt DO warming transitions is depicted in Fig. 6. where it is much more gradual than in the other records (possibly because calcite sampling was not perfectly perpendicular to the stalagmite isochrones, or because growth rates were variable in between the U / Th ages). the transitions we manually determine pre-event and postevent averages, as indicated by the orange lines. The averaging time is set to 150 and 50 years for stadial and interstadial periods, respectively; this difference in duration is used because (i) several of the interstadials are of short duration and (ii) Greenland δ 18 O is more variable during stadial climates, requiring longer averaging. For DO 16.1, the duration of the pre-event stadial baseline climate was shorter than 150 years, and the averaging time was reduced to 100 years (Fig. 6).
After determining the pre-and post-event averages, we use linear interpolation of the time series to find the time at which the variable of interest had completed 25, 50, and 75 % of the total transition (Fig. 6). We use the 50 % marker (red) as the midpoint of the transition, which is used in the methane synchronization. The 25 and 75 % markers (blue) are used as the ±1σ uncertainty estimate. In rare cases the time series contain inversions within the transitions that lead to ambiguity in the timing of the markers; for these events we find the markers using a monotonic spline fit to the data.
The midpoints of abrupt interstadial terminations were determined in the same fashion (WD CH 4 and NGRIP only). Tables 1 and 2 give the results for NH warming and NH cooling, respectively.

Synchronizing WD to a NGRIP-Hulu hybrid chronology
Abrupt DO variability is expressed clearly in a great number of NH climate records (Voelker, 2002). For the purpose of methane synchronization, our interest is in highresolution records that express the abrupt DO events very clearly, and are furthermore exceptionally well dated. We here use a combination of two such NH records (Fig. 5), namely the Greenland NGRIP δ 18 O record (NGRIP community members, 2004), and a refined version of the Hulu Cave speleothem δ 18 O record (Edwards et al., 2015;Reimer et al., 2013;Southon et al., 2012) with improved resolution and additional dating constraints (see Wang et al., 2001, for the original, lower resolution Hulu δ 18 O record). The DO events are resolved most clearly in the NGRIP δ 18 O record, which is available at 20-year resolution. We use the GICC05modelext chronology for this core, which is based on annuallayer counting back to 60 ka BP and ice-flow modeling for ice older than 60 ka Svensson et al., 2006;Wolff et al., 2010). While annual-layer counting provides accurate relative ages (e.g., the duration of DO interstadials), it provides relatively inaccurate absolute ages due to the cumulative nature of counting uncertainty (Table 1). The refined Hulu δ 18 O record also shows the abrupt DO events in high temporal resolution (Fig. 6). The speleothem chronology is based on U / Th radiometric dating, providing much smaller uncertainty in the absolute ages than GICC05 (Table 1). The reason for selecting this record over other speleothem records is the large number of U / Th dates, the low detrital Th at the site, and the high sampling resolution of the δ 18 O record (Wang et al., 2001). In the Hulu data, as in other records of DO variability, the interstadial onsets are more pronounced and abrupt than their terminations. We therefore only use the timing of the former as age constraints, as they can be established more reliably. The onset of NH interstadial periods as expressed in Hulu δ 18 O is given in Table 1.
In both the NGRIP and Hulu Cave δ 18 O records we have determined the ages of the midpoints of the DO transitions ( Fig. 6; Table 1); a plot of their difference (Hulu age minus NGRIP age) is shown in Fig. 7, where the error bars denote the root sum square of the NGRIP and Hulu midpoint determination uncertainty (Sect. 4.3). The Hulu ages are systematically older than the NGRIP ages, and the age difference increases going further back in time. Note that the Hulu-NGRIP age difference is smaller than the stated GICC05 counting uncertainty (832 to 2573 years) but larger than the Hulu age uncertainty (92 to 366 years). A linear fit through these data, forced to intersect the origin, is given by   Table 1). The GICC05 ages are placed on the BP 1950 scale rather than the b2k scale (years prior to 2000 CE). 0.0063 × GICC05 age, suggesting that the GICC05 annuallayer count on average misses 6.3 out of every 1000 layers. Because of this observation we use a linearly scaled version of the GICC05 chronology (GICC05 × 1.0063) as the target chronology for methane synchronization. This approach has several advantages. First, it respects both the superior relative ages (i.e., interval durations) of GICC05, as well as the superior absolute ages of the Hulu chronology. Second, it is very simple to convert between the WD2014 and GICC05 chronologies (CH 4 -synchronized section of the chronology only); one simply needs to divide WD2014 ages by 1.0063 (and add 50 years to convert to the b2k reference date). Third, it still allows for direct synchronization of WD CH 4 to the NGRIP δ 18 O record, providing more tie points than direct synchronization to the Hulu record would. Note that the GICC05 × 1.0063 target chronology respects the Hulu age constraints in an average sense only; the age of individual events differs between Hulu and our target chronology by up to 180 years. Our approach therefore represents only a firstorder correction of a growing offset between GICC05 and Hulu; nonlinear temporal changes in the counted dating error may exist from one tie point to the next (Fleitmann et al., 2009).
The exercise of finding the transition midpoints and determining the GICC05-Hulu scaling factor was performed by two of the authors (J. P. Severinghaus and C. Buizert), independently of each other. The scaling factors obtained were 1.0063 and 1.0064, respectively, showing that, to first order, this result is insensitive to (subjective) judgment in identifying the transitions. The difference between the Hulu ages and 1.0063 × GICC05 ages are all within the stated Hulu 2σ dating error (Table 1). Consequently, our chronology is not in violation of any Hulu constraint as it respects the Hulu 2σ error at all of the tie points. In deriving the scaling we have assumed that the abrupt DO transitions observed in NGRIP and Hulu are simultaneous, which is not necessarily true. The variations in monsoon intensity represented by Hulu δ 18 O are commonly explained by meridional movement of the Intertropical Convergence Zone (ITCZ) and tropical rainfall belts (Wang et al., 2001(Wang et al., , 2006Kanner et al., 2012); modeling work suggests such atmospheric readjustments occur on decadal timescales in response to NH high-latitude forcing (Chiang and Bitz, 2005;Broccoli et al., 2006;Cvijanovic and Chiang, 2013). Moreover, CH 4 emission changes are nearsynchronous with Greenland δ 18 O variations, which they lag by only a few decades on average Baumgartner et al., 2014;Rosen et al., 2014). Since both CH 4 emissions and Hulu δ 18 O are closely linked to tropical hydrology, timing lags between NGRIP and Hulu are also expected to be on decadal timescales. The uncertainty in the NGRIP-Hulu phasing is therefore probably small (decadal) relative to the correction we apply (up to 400 years).
Rather than synchronizing WD CH 4 to Greenland CH 4 records, we have chosen to synchronize directly to NGRIP δ 18 O, which varies in phase with CH 4 (but with a nearly constant time lag). We let the midpoint in the CH 4 transitions lag the midpoint in the NGRIP δ 18 O transition by 25 years, as suggested by studies of Greenland δ 15 N-CH 4 phasing Baumgartner et al., 2014;Rasmussen et al., 2013;Kindler et al., 2014;Rosen et al., 2014). The rationale behind this approach is threefold. First, throughout MIS 3 the NGRIP δ 18 O record has both better precision and higher temporal resolution than any available Greenland CH 4 record Brook et al., 1996;Blunier et al., 2007). Second, the dating of Greenland gas records depends on the highly variable age function, which is not equally well constrained for all DO events (Schwander et al., 1997;Rasmussen et al., 2013). This reliance on Greenland age would introduce an additional source of uncertainty. The NGRIP δ 18 O record, on the other hand, is accurately dated through the GICC05 layer count. Third, Greenland CH 4 records are more strongly impacted by firn smoothing than the WD CH 4 record, because glacial accumulation is lower in Greenland (Greenland glacial age is about 2-3 times as high as WD age during that time). In summary, our approach circumvents the uncertainties associated with using Greenland CH 4 as an intermediary, or, to state this another way, the uncertainty in the phasing between CH 4 and Greenland δ 18 O is smaller than the uncertainty in the Greenland age.

Interpolation between age constraints
We can assign a gas age to each of the depths where an abrupt WD CH 4 transitions occurs; we do this for DO 4.1 through DO 18, i.e., the events prior to 31.2 ka BP (the onset of the WD layer count). The gas age we assign is equal to 1.0063 times the GICC05 age for the same event, with 25 years subtracted to account for the slight CH 4 lag behind Greenland δ 18 O. By adding age (Sect. 3) to this gas age we assign an ice age. These assigned ages are printed in boldface in Tables 1 and 2. To obtain a continuous depth-age relationship between these ice age constraints, we have to apply an interpolation strategy. This task amounts to estimating the annual-layer thickness λ(z) along the deep part of the core. The simplest approach is to assume a constant accumulation rate in between the age constraints; this is shown in Fig. 8b for the case where we use the age constraints from NH warming events only (black) or the age constraints from both NH warming and cooling events (red). The disadvantage of this approach is that it results in discontinuities in λ(z) (the first derivative of the depth-age relationship), which we consider highly unrealistic. A more realistic approach is therefore to assume that λ(z) is continuous and smooth ); Fig. 8b shows two scenarios in which we use a spline function to estimate λ(z), where again we have applied age constraints from NH warming events only (orange) or age constraints from both NH warming and cooling events (blue).
For comparison, past A obtained from the firn densification model (Sect. 3) is plotted in green (Fig. 8b). While the δ 15 N-based A follows the synchronization-based A estimates broadly, the millennial-scale details do not agree. We want to point out that this is not unexpected, since both methods have their imperfections. In particular, any errors in the (stretched) GICC05 age model or in our modeled thinning function f λ (z) will strongly impact the synchronization-based A estimates in Fig. 8b. The discrepancy is pronounced between 60 and 65 ka, where we have to strongly reduce λ(z) in order to fit the age constraint(s) from DO 18, while δ 15 N provides no evidence for low A during this interval.
For the WD2014 chronology we have applied the smooth λ(z) interpolation scheme using all age constraints (i.e., both NH warming and cooling events). The midpoint detection uncertainty is comparable for all events and systematically smaller at the start of interstadial periods than at the terminations (Tables 1 and 2). For short interstadials (e.g., DO 9) this leads to a large relative uncertainty in the event duration, and thereby a large uncertainty in the implied accumulation rates (Fig. 8b). We force the interpolation to fit all NH warming constraints perfectly, yet relax this requirement for NH cooling constraints to prevent large swings in λ(z) for the short-duration events. The WD2014 chronology fits the NH warming and NH cooling age constraints with a 0-and 16-year rms offset, respectively. Because the duration of (inter)stadial periods is well constrained in the layer-counted GICC05 chronology, using both the NH warming and NH cooling tie points results in a more robust chronology. The duration of (inter)stadial periods is 0.63 % longer in WD2014 than in GICC05, which is well within the stated GICC05 counting error of 5.4 % (31.2-60 ka interval). , 11, 153-173 annual-layer thickness scenarios, converted to an accumulation rate for comparison to the δ 15 N-based firn model reconstructions. The interpolation strategy is to use either constant accumulation rates between tie points ("constant") or a smoothly varying λ(z) ("smooth"); the age constraints used are either only the NH warming events ("warming"), or both the NH warming and cooling events ("all"). (c) Estimated 2σ uncertainties in the WD2014 chronology due to age, choice of interpolation scheme, midpoint detection, and the absolute age constraints used in the synchronization. Total absolute ice age uncertainty plotted in solid black; relative age uncertainty (i.e., with absolute age uncertainty in the Hulu-GICC05 master chronology withheld) plotted in dashed black.

Age uncertainty
The age uncertainty we assign to the deep part (> 2850 m) of the WD2014 chronology has four components.
The first source of uncertainty is the age calculation; we use the 2σ uncertainty obtained in the age sensitivity study (Sect. 3.2). The second source of uncertainty is the choice of interpolation scheme used to obtain a continuous chronology; here we use the standard deviation between the four different interpolation schemes of Fig. 8b as an uncertainty estimate. The third source of uncertainty is the difficulty in determining the timing of the abrupt events in the time series; we use the uncertainty in the midpoint evaluation (root sum square of WD CH 4 and NGRIP δ 18 O estimates). The last source of uncertainty is the age uncertainty in the hybrid NGRIP-Hulu chronology that we synchronize to. We use the stated Hulu age uncertainty plus 50 years to account for possible leads or lags in the NGRIP-Hulu δ 18 O phasing, plus the absolute value of the offset between the Hulu ages and the 1.0063 × GICC05 ages. For DO events where we do not have reliable Hulu age estimates (Table 1), we set the uncertainty to the Hulu age uncertainty of the nearest event, plus the uncertainty in the interval duration specified by the GICC05 layer count. For example, for DO 14 we do not have a reliable Hulu age estimate, and we use the Hulu age uncertainty of DO 16.2 (226 years) plus the uncertainty in the DO 14 to DO 16.2 interval duration on GICC05 (209 years), giving a total of 226 + 209 = 435 years. The uncertainties (2σ values) are plotted in Fig. 8c (log scale). We assume these four uncertainties to be independent, and use their root sum square as the total uncertainty estimate on the WD2014 ice age scale (Fig. 8c, black curve). Note that the fourth source of uncertainty is only relevant when considering absolute ages; when evaluating relative ages (e.g., between WD ice and WD gas phase, or between WD and NGRIP), this last contribution does not need to be considered. For the deepest WD ice (3404 m depth) we thus find an age of 67.7 ± 0.9 ka BP.

Discussion
While the WAIS Divide ice core does not extend as far back in time as deep cores from the East Antarctic Plateau, its relatively high temporal resolution (due to the high snow accumulation rate) makes it an ice core of great scientific value. WD accumulation rate during the LGM (∼ 10 cm a −1 ice equivalent) is still higher than the present-day accumulation rate at the EPICA (European Project for Ice Coring in Antarctica) Dronning Maud Land core (7 cm a −1 ), which is generally considered a high-accumulation core (EPICA Community Members, 2006). With 68 ka in 3404 m of core, the core average λ is 5 cm a −1 , at the onset of the last deglaciation (18 ka BP) λ is around 4 cm a −1 , and near the bed λ is around 0.4 cm a −1 . This high temporal resolution provides the opportunity for obtaining very detailed climatic records.
High accumulation rates also result in a small age. Figure 9 compares age between several Antarctic cores (note the logarithmic scale). age at WD is approximately onethird of the age at EPICA DML (EDML) and Talos Dome (TALDICE), and one-tenth of the age at EPICA Dome C (EDC), Vostok, and Dome Fuji. Because the uncertainty in the age (or depth) calculation is typically on the order of 20 %, a smaller age allows for a more precise interhemispheric synchronization with Greenland ice core records using CH 4 . The small WD age uncertainty during MIS 3 allows for investigation of the phasing of the bipolar seesaw (Stocker and Johnsen, 2003) at sub-centennial precision (WAIS Divide Project Members, 2015).
In comparing the shape of the age profiles, there are some interesting differences (Fig. 9). It is important to realize that not all the age histories shown were derived in the same way; WD and Dome Fuji age were derived using densification models, and the other four were derived using the depth approach (Parrenin et al., 2012) and a Bayesian inverse method that includes a wide range of age markers (Veres et al., 2013). We will therefore focus on comparing the WD and Dome Fuji results. age at WD shows more pronounced variability than at Dome Fuji, particularly during MIS 3. The reason is that the glacial firn pack at Dome Fuji is about 4000 years old, and consequently the firn column integrates over 4000 years of climate variability, thereby dampening the age response to millennial-scale climatic variability. At WD the glacial firn layer is only about 350 years old, and therefore the firn is in near equilibrium with the millennial-scale climate variations. This difference in response time is also obvious during the deglaciation, where WD age transitions from glacial to interglacial values between 18 and 14.5 ka BP, while Dome Fuji takes more time (18-10 ka BP). Surprisingly, EDML age does not show a strong deglacial age response, unlike all the other cores.
The relatively small age at WAIS Divide also allows for precise investigation of the relative timing of atmospheric greenhouse gas variations and Antarctic climate (Barnola et al., 1991;Pedro et al., 2012;Caillon et al., 2003;Parrenin et al., 2013;Ahn et al., 2012). Recent works suggest that during the last deglaciation the rise in atmospheric CO 2 lagged the onset of pan-Antarctic warming by approximately 0 to 400 years (Pedro et al., 2012;Parrenin et al., 2013). This Antarctic warming around 18 ka BP is presumably driven by the bipolar seesaw, as it coincides with a reduction in Atlantic overturning circulation strength as seen in North Atlantic sediment records (McManus et al., 2004). The WD age at 18 ka (gas age) is 515 ± 91 years (2σ ), much smaller than at central East Antarctic sites such as EPICA Dome C, where age is approximately 3850 ± 900 years (Veres et al., 2013, with the age uncertainty taken to be the difference between the gas age and ice age uncertainties). The precision with which one can determine the relative phasing of climatic (i.e., δ 18 O of ice) and atmospheric signals is set by the uncertainty in age (or equivalently, the uncertainty in depth). High-resolution WD records of CO 2 and CH 4 (Marcott et al., 2014) place the onset of the deglacial rise in the atmospheric mixing ratio of these greenhouse gases on the WD2014 chronology at 18 010 and 17 820 years BP, respectively. However, evaluating the relative phasing of CO 2 and Antarctic climate is complicated by the observation of asynchronous deglacial warming across the Antarctic continent (WAIS Divide Project Members, 2013). Attempts to capture the climate-CO 2 relationship in a single lead-lag value may be an oversimplification of deglacial climate dynamics.
An important next step will be to synchronize the WD chronology with other Antarctic cores via volcanic matching and other age markers (e.g., Severi et al., 2007;Sigl et al., 2014). Because of the annual-layer count and possibility of tight synchronization to Greenland ice cores, WD could contribute to an improved absolute dating of Antarctic cores, as well as improved cross-dating between cores. Such cross-dating could help inform the WD chronology as well, particularly in the deepest part of the core, where the ice is potentially highly strained, as suggested by the interpolation difficulties in the 60-65 ka interval (Fig. 8b). With a synchronized chronology, WD could improve the representation of West Antarctic climate in Antarctic ice core stacks (Pedro  Parrenin et al., 2013), and provide a more refined pan-Antarctic picture of the climate-CO 2 relationship.

Summary and conclusions
We have presented a first chronology for the deep (> 2850 m) section of the WAIS Divide ice core, which is based on stratigraphic matching to Greenland ice cores using globally well-mixed methane. We use a dynamical firn densification model constrained by δ 15 N data to calculate past age, and find that age was smaller than 525 ± 120 years for all of the core. Using high-resolution WD records of atmospheric CH 4 , we synchronize WD directly to Greenland NGRIP δ 18 O for the abrupt onset and termination of each of the DO interstadials. To each event we assign an age corresponding to 1.0063 times its GICC05 age, which brings the ages in agreement with the high-resolution U / Th-dated Hulu speleothem record. The uncertainty in the final chronology is based on the uncertainties in (i) the age calculations, as evaluated with a sensitivity study; (ii) the interpolation strategy, as evaluated by comparing four different interpolation methods; (iii) determining the timing of events in the Due to the combination of a small age and a highresolution methane record, the WAIS Divide ice core can be synchronized more precisely to Greenland records than any other Antarctic core to date. This is important when investigating interhemispheric climate relationships such as the bipolar seesaw. The small WD age furthermore provides valuable opportunities for precise investigation of the relative phasing of atmospheric greenhouse gas variations and Antarctic climate.

Appendix A: Densification physics
The densification rates used in this work are based on the empirical steady-state model by Herron and Langway (1980) (the H-L model). We use the H-L model with minor modifications that allow it to be run dynamically (i.e., with timevariable T and A) and to include the softening effect of impurities following Freitag et al. (2013a). The H-L model divides the firn column in two stages, separated at the critical density ρ c = 550 kg m −3 , occurring at the critical depth z c .
For the upper firn (ρ ≤ ρ c , stage 1), the densification rates are given by with where E 1 = 10.16 kJ mol −1 is the activation energy for stage 1 and R is the universal gas constant. Because both the sinking velocity of deposited layers (w = dz/dt) and the densification rate scale linearly with A, the resulting density-depth profile ρ(z) in stage 1 becomes independent of A, and sensitive to T variations only. For the deeper firn (ρ > ρ c , stage 2), we use Eq. (4c) from Herron and Langway (1980), which was first derived by Sigfús J. Johnsen. This equation gives the densification rate in terms of overburden load, which allows the model to be run dynamically. The stage 2 densification rates are given by with where E 2 = 21.4 kJ mol −1 is the activation energy for stage 2 and σ z denotes the firn overburden load at a given depth in Mg m −2 : Note that we divide by 1000 to convert from kg m −3 to Mg m −3 , the units used by Herron and Langway (1980). We use the mathematical description by Freitag et al. (2013a) to include the hypothesized firn softening effect of impurities. In this approach an increasing Ca concentration, as a proxy for mineral dust content, lowers the activation energy of firn, thereby enhancing densification rates. This is tantamount to stating that dusty firn behaves as if it were "warmer" than its climatological temperature. The H-L activation energies of Eqs. (A2) and (A4) are modified by [Ca] in the following way: where E Ca and E HL are the Ca-modified and original H-L activation energies, respectively, [Ca] crit = 0.5 ng g −1 is the minimum concentration at which impurities affect densification, and α and β are calibration parameters. Whenever [Ca](z) < [Ca] crit , we set [Ca](z) = [Ca] crit . The parameter β sets the sensitivity to dust loading, and α is a normalization parameter that is included to account for the fact that the original H-L model was calibrated without the impurity effect. Consequently, if β > 0, one needs to compensate by setting α > 1 to preserve the original H-L calibration. The work by Freitag et al. (2013a) recommends β = 0.01 and α = 1.025 (which yields E Ca = E HL at [Ca] = 5.73 ng g −1 ).
Using the recommended value of α = 1.025 at WD provides a poor fit to observations of present-day firn density and close-off depth. The optimal fit to present-day WD observations is obtained using an activation energy equal to 1.007 × E HL ; this is in between the values suggested by Herron and Langway (1980) and Freitag et al. (2013a). In the experiment presented in Fig. 4 we changed the dust sensitivity β; it is clear that we need to simultaneously change α to keep the model well-calibrated to present-day conditions. Due to the fact that the mean late Holocene WD [Ca] is around 0.8 ng g −1 , we let α = 1.007/(1−β ln[0.8/0.5]) in the experiment of Fig. 4. This approach ensures that the present-day E Ca is invariant with β, and equals E Ca = 1.007×E HL . This means that whatever value we choose for β, we will obtain a good fit to the present-day age, δ 15 N, and A values that are well known from direct observations (Battle et al., 2011).
To validate the H-L model age simulations, we repeated the firn modeling using the densification physics of Arnaud et al. (2000), which is also the basis of the model by Goujon et al. (2003). Our implementation of the Arnaud model is based on the description in the latter paper, with one modification at the critical density that we outline here. In the Arnaud model, densification in the stage 1 follows the work of Alley (1987), and is given by with D the relative density D = ρ/ρ ice , P the overburden pressure, and γ a scaling factor used to make the densification rates continuous across the critical density D c . Stage 2 densification is given by with where a is the average contact area between the grains, Z is the coordination number, and E A is the activation energy (60 kJ mol −1 ). Arnaud densification rates for stage 3 (D ≥ 0.9) are describe elsewhere (Goujon et al., 2003;Arnaud et al., 2000). The difficulty in implementing this model is the following. The densification rates of Eqs. (A7-A8) exhibit a discontinuity at the critical density D = D c = 0.6 that cannot be remedied with the scaling factor γ . On approaching D c , densification rates given by Eq. (A7) go to zero (due to the inclusion of the term (1 − 5 3 D), while densification rates given by Eq. (A8) go to infinity because the contact area a equals zero at D = D c . Clearly neither equation gives a realistic result at D = D c . Therefore, in our implementation of the Arnaud model we use the H-L densification rates of Eq. (A1) instead of Eq. (A7) in stage 1. We take the onset of stage 2 to be the density at which Eqs. (A1) and (A8) intercept, thus avoiding the singularity in Eq. A8. This approach has the additional advantages of removing dependence on ad hoc scaling factor γ and introducing realistic temperature dependence for stage 1 densification. Because stage 1 spans just the top 10-20 % of the firn column, the modification has only a minor influence on the overall behavior of the Arnaud model. The Goujon model code avoids the singularity in Eq. (A8) by extending stage 1 to D c + ε (Anaïs Orsi, personal communication, 2014), a procedure not described in Goujon et al. (2003).
The Supplement related to this article is available online at doi:10.5194/cp-11-153-2015-supplement.