A first chronology for the NEEM ice core

Introduction Conclusions References


Introduction
The 2540 m long North Greenland Eemian Ice Drilling (NEEM) ice core was drilled during 2008-2010 through the ice sheet in Northwest Greenland (77.45 • N,51.07 • W, annual layer identification in the NEEM core may be undertaken either by the manual identification of annual layers as performed on data from the DYE-3, GRIP, and NGRIP cores in the construction of the Greenland Ice Core Chronology 2005, hereafter GICC05, , or by automated methods . 10 At this point, however, data for annual layer identification are not available. Thus, to provide a first time scale for the NEEM core, the GICC05 time scale and its model extension GICC05modelext is adopted as is from the NGRIP core using common reference horizons, or match points, as identified from the Electrical Conductivity Measurement (ECM), Dielectrical Profiling (DEP), and tephra records of the two cores.

The GICC05 and GICC05modelext time scales
Annual layers can be identified in all deep ice cores from Greenland, and due to differences in accumulation rates, ice thickness, temperature, ice flow conditions, and availability of data in sufficient resolution, different deep ice core records are best suited for annual layer identification in different age or depth ranges. For example, with a present-10 day annual accumulation of 0.56 m of ice equivalent, the DYE-3 ice core is optimal for dating the most recent millennia as the DYE-3 δ 18 O record contains a well-developed annual signal reflecting the seasonal temperature variation, while in return, the high accumulation rates lead to rapid thinning of layers at greater depths . In contrast to DYE-3, the present-day annual accumulation at the NGRIP site of 0.19 m 15 of ice equivalent means that the NGRIP δ 18 O record only marginally resolves the annual layers even at the top of the core. However, the combination of the availability of high-resolution Continuous Flow Analysis (CFA) impurity records and comparably low thinning rates in the glacial section due to bottom melting, permit the annual layers of more than 1 cm thickness to be identified back to approximately 60 ka b2k (Svensson 20 et al., 2008). The idea behind the design of the GICC05 time scale is to exploit these differences, creating a multi-core time scale based on layer counting, using those cores and records that optimally resolve the annual layers at a given time period. The resulting time scale is based on identification of annual layers in δ 18 O records in the top section, using 25 parallel data from DYE-3 (most recent 8.3 ka), GRIP (most recent 3.8 ka), and NGRIP (most recent 1.8 ka) linked together using common volcanic signatures (Vinther et al., 2971 between the different data series ensure maximum consistency. In the glacial section, annual layer thicknesses decrease further due to low glacial accumulation rates and flow-induced thinning and therefore only the NGRIP CFA, ECM and visual stratigraphy (VS) data can be used for annual layer identification (Andersen et al., 2006b;Svensson et al., 2006Svensson et al., , 2008. Beyond 60 ka b2k, the counted GICC05 time scale is extended with the modelled ss09sea06bm time scale (North Greenland Ice Core Project members, 2004) shifted 705 yr towards younger ages to ensure continuity (Wolff et al., 2010). The combined time scale has later been named GICC05modelext, and annual layer counting across Greenland Stadial 22 (GS-22) using new high-resolution CFA data largely confirms the validity of the GICC05modelext 15 estimates of the duration of GS-22 (Vallelonga et al., 2012).
With the aim of providing a common chronological framework for the Greenland cores, the GICC05 time scale has been transferred to the Renland and Agassiz ice cores (Vinther et al., 2008) and to the GRIP and GISP2 ice cores in the interval 32.45-11.7 ka b2k  and further extended to 48 ka b2k for GRIP in 20 Blockley et al. (2012) using an approach similar to the one applied here. Ongoing work will extend the GICC05modelext to the full length of the stratigraphically uncompromised sections of the GRIP and GISP2 cores (Seierstad et al., 2013), and GICC05 ages are also used in integrated inverse modelling efforts aiming at constructing a common chronological framework for Antarctic and Greenland ice cores (Lemieux-Dudon Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | increasingly marginal for resolving annual layers, and even when the GRIP CFA data becomes available in the section older than 7.9 ka, the robustness of the annual layer count is adversely influenced by the fact that the GRIP CFA set consists of only two independent parameters at moderate resolution. The section from about 10.3 to 7.5 ka b2k (after which the NGRIP CFA measurements become available) is the weakest sec-5 tion data-wise, and although this time period is partially overlapping with the brittle zone in the NEEM core, there is potential for a strengthening of the annual layer identification across this interval once NEEM data is available. Another potential addition to the GICC framework based on NEEM data is semi-automatic and objective annual layer identification using VS and CFA data, possibly in parallel between the NEEM and 10 NGRIP cores, e.g. using the method developed by Winstrup et al. (2012).
The uncertainty in GICC05 annual layer identification is given as the so-called Maximum Counting Error (MCE) derived from the number of uncertain annual layers, each counted as 1/2 ± 1/2 yr. Comparisons with other palaeoclimate records indicate that the GICC05 has little or no systematic bias and confirms that the MCE probably is 15 a conservative estimate of the total age uncertainty of the time scale (Svensson et al., 2008). In this work, we follow the same notation and report all GICC05 ages with corresponding MCE values. It should be stressed that the MCE is not a Gaussian uncertainty measure and that MCE values strictly speaking cannot be compared to Gaussian-style dating uncertainties of other records, but in the absence of a more appropriate un-20 certainty estimate, we recommend that the MCE is regarded as the 2σ uncertainty of GICC05 in cases where Gaussian uncertainties are needed.

Electrical Conductivity Measurements (ECM)
Electrical Conductivity Measurements (ECM) were conducted in the field during the NGRIP (Dahl-Jensen et al., 2002) and NEEM ice core campaigns. For both cores, the 25 ECM set-up and procedure was similar to the one described in Hammer (1980): first, a flat surface along the depth axis of the ice core is polished with a microtome knife to produce a level and freshly cleaned surface. The hand-held ECM instrument consists 2973 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of a pair of brass electrodes spaced ∼ 1 cm apart, between which a DC voltage difference of 1250 V is maintained. This device is moved along the ice in a steady fashion, while maintaining a constant pressure of the electrodes onto the ice surface. The electrical field causes H + -ions to be displaced resulting in a transient displacement current which is measured. Steady movement of the ECM instrument over the ice surface pre-5 vents the ice from being polarized leading to a reduced current over time (Hammer, 1980). Based on comparisons to high-resolution chemistry measurements, it has been established that the recorded ECM signal is almost entirely a response to the acidic components of the ice, even in the presence of large concentrations of neutral salt (Moore et al., 1992). Hence, ECM is a fast and non-destructive technique for providing 10 a high-resolution acidity profile for the ice core. The ECM technique is often used as one of the primary tools for establishing an event stratigraphy for volcanic eruptions, which are detectable as high-acidity layers in ice cores. ECM was recorded continuously along the NGRIP and NEEM ice cores. For NGRIP, the upper part (down to 1372 m) of the ECM profile was measured on the NGRIP1 core, whereas the lower part (1346-3085 m) is based on data from the NGRIP2 core. Inaccuracies in the logging of the NGRIP cores have given rise to a small offset (∼ 0.43 m) between corresponding volcanic events recorded in the two cores in the overlapping section, with NGRIP1 data being assigned the greater depth (Hvidberg et al., 2002).

Calibration
The ECM signal is temperature dependent and was calibrated to a standard temperature of 14 • C by using an Arrhenius law with an activation energy of 0.23 eV (Neftel et al., 1985). The exact conversion between ECM and ice acidity is ambiguous. Even for a specific 25 set of electrodes, the data display a large amount of scatter, and several calibration curves have been suggested based on measured H + ion concentrations. The initial calibration curve suggested by Hammer (1980) uses a square-law relationship, but 2974 Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the appropriate parameters in this calibration seem to be ice-core dependent (Moore et al., 1992). This is partly a result of different ambient conditions and operational differences (e.g. pressure applied to the electrodes) during measurements. Also, the exact procedures concerning storage and cleaning of the ice core surface may play a role. For the NGRIP and NEEM cores, many operators have been in charge of the 5 ECM measurements, and the ice core was subjected to varying conditions. Hence, the calibration is likely to be variable with depth, and thus an elusive target indeed. Here, the ECM current i (in µA) has been converted to ice acidity (in µequiv. H + kg −1 ) using the relationship [H + ] = 0.045 · i 1.73 as suggested by Hammer (1980).
For our purposes, however, a further fine-tuning of the ECM calibration is unneces-10 sary: the synchronization of the two cores does not rely on absolute values, but solely on the recognition of similar patterns of peaks in the two ice core acidity records.

Data quality
The quality of the ECM measurements is generally higher in ice than in the firn, where the lower density makes it difficult to maintain constant contact area between ice crys-15 tals and electrodes. This leads to generally lower currents and enhanced scatter in data from the uppermost part of the ice core.
To ensure the quality of the ECM data, the ECM profile was measured twice for each core section and the two profiles were subsequently checked for agreement. In case of major dissimilarities, measurements were repeated until reproducible results were 20 obtained and a representative run could be selected. In this way, the reproducibility of the ECM data was continuously checked.
The manually-operated ECM device is fast and simple, and it allows the scientific investigator to take advantage of the best-preserved section of the ice core surface by manoeuvring the electrodes around cracks or other areas of poor-quality ice.
Breaks across the entire core diameter cannot be avoided, and can be recognized as abrupt dips in the ECM profile. The location of breaks was logged as part of the ECM Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | measurements, and data from these sections have subsequently been removed from the resulting ECM data sets. A disadvantage of the handheld electrode approach is that the depth assignment is only accurate to 1-2 cm, which is a consequence of the flexibility in the electrode setup together with the uncertainty of the registration of the movement of the electrodes 5 along the length of the core. The depth uncertainty of 1-2 cm of the ECM data makes a detailed comparison with e.g. high-resolution chemistry measurements of the ice core challenging. Data have subsequently been interpolated to a depth resolution of 1 cm for NGRIP1 and NEEM data and 1 mm for NGRIP2 data, with the different resolution reflecting changes in the data post-measurement treatment.
Despite the rather low accuracy of the depth assignment, the high resolution of the ECM records is important for the synchronization of the ice cores, as it permits robust pattern recognition of features found in both cores.

Dielectrical Profiling (DEP)
The dielectrical stratigraphy of the NEEM core was recorded in the field with the Dielec- 15 tric Profiling (DEP) technique (Wilhelms et al., 1998). Ice core sections of 1.65-2.2 m were usually analysed no more than one day after drilling. The brittle zone ice could not be handled directly after drilling and was stored and processed during the subsequent season, but all the data used here comes from below the brittle zone. Approximately every day a free air measurement was recorded for data calibration. 20 For each measurement the cable stray capacitance and conductance -as determined in the corresponding daily free air measurements -were subtracted from the respective capacitance and conductance records. After multiplying the conductance with the vacuum permittivity, permittivity and conductivity are calculated by division with the free air capacitance. 25 The relative permittivity increases during the initial densification in the upper part of the core and remains just above 3 in deeper core sections. Therefore, too low permittivity is a good control for sections with bad core quality influenced e.g. by breaks or 2976 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | missing pieces. Data from these sections are not trustworthy, and a corrected data set is created by removing sections where the permittivity drops below a chosen threshold (see Fig. 1). The NEEM DEP data set is provided in processed form without correction for the temperature variation in the science trench or in the core-buffer. Therefore, care should 5 be taken when interpreting the absolute levels between sections. For the identification of patterns of peaks between different cores, however, this limitation is not important.
(A short section on any NGRIP DEP data processing (focusing on any differences in measurement and/or processing relative to NEEM) will be added here but was not available in time for the IPICS special issue deadline.)

Tephra horizons
Glacial ice from both NEEM and NGRIP ice cores is currently being intensively investigated with the aim of developing a high-resolution tephrochronological framework, which is able to link the cores together and to provide a robust template for correlating the ice-core records with other sediment records containing tephras (Lowe et al., 2008;15 Blockley et al., 2012). In this work we focus on the identification of cryptotephra horizons, i.e. horizons that contain a low concentration of volcanic glass particles and are invisible to the naked eye. A sampling strategy was devised based on three criteria: intervals corresponding to age ranges of known tephra layers not yet identified in the ice cores.
In the sampled sections, ice strips with a cross section of approximately 2 cm 2 were cut from the 55 cm long archive segments of the core stored at the University of 25 Copenhagen. The sample preparation methodology followed that outlined in Davies 2977 Introduction  (2008,2010) and , and involves centrifugation of the melted ice sample and preparation of the remaining particulate material onto ground glass slides. The samples were then embedded in epoxy resin and high-magnification light microscopy was undertaken to identify the presence of glass-shard particles. Samples found to contain volcanic material were then ground and polished using different 5 grades of diamond suspension (9, 6 and 1 µm) to produce a polished surface suitable for analysing single shards by electron microprobe. Electron-probe microanalysis (EPMA) of the identified glass shards took place at the Tephrochronology Analytical Unit based at the University of Edinburgh. A Cameca SX-100 electron microprobe with five vertical wavelength dispersive spectrometers was 10 used to analyse oxide values for 10 major and minor elements within individual glass shards. The operating conditions followed those outlined by Hayward (2012).
Fourteen tephras common to the NGRIP and NEEM ice cores are identified and are outlined in Table 1. These are Icelandic in origin with 7 of tholeiitic basaltic composition and 7 of transitional alkali basaltic composition (Fig. 2). Similarity coefficient values for 15 all identified tephra tie-points are high (> 0.96), supporting the proposed correlations (Table 1). Similarly, the statistical distance values do not exceed the critical value and are not considered to be statistically different.

Methods and results
In the logging and analysis of ice cores, the core length is the fundamental param-20 eter. The core is retrieved from a slightly inclined bore hole, and inaccuracies in the logged length accumulate along the core, especially in sections where drilling difficulties resulted in poor core quality; or in the brittle zone, where the core quality is often compromised by depressurisation of air bubbles within the ice. However, as all measurements on the core are made relative to the same logged depth scale, any such 25 inaccuracies are of no relevance to the synchronization of the ice cores or time scale transfer. A potential offset between true depth and logged depth could be of relevance when records from the ice core are compared to in-hole measurements or modelled quantities. We have no direct way of estimating the magnitude of this offset, but the NEEM drilling and logging procedures are similar to those employed at NGRIP where the drilling and logging of two parallel ice cores down to a depth of 1372 m made it possi-5 ble to evaluate the NGRIP2 logging precision to about 0.5 m at 1372 m depth (Hvidberg et al., 2002). As the estimated uncertainty on the logged depth therefore is small and varies slowly with depth, the effect on the comparison of true-depth ice flow modelling results and core-depth measurements is negligible. 10 To allow for the investigation of differences, including possible leads and lags, between the climate records of the two cores, the transfer of the time scale is based on reference horizons of non-climatic origin only. Ideally, only horizons that can be uniquely identified in both cores should be used. Extensive ongoing work on locating tephra in ice cores and geochemically characterizing these tephras is providing an ever-growing number 15 of such horizons that tie different ice cores together with high confidence (Mortensen et al., 2005;Narcisi et al., 2005Narcisi et al., , 2012Coulter et al., 2012). However, the number of available tephra horizons is still far from sufficient to provide the sole basis for the transfer of GICC05 to the NEEM core. Instead, the synchronization is based on the identification of common patterns in the continuous and highly resolved ECM and DEP 20 records with the tephra horizons providing an independent check of the synchronization.

Synchronization of the NEEM and NGRIP cores
Peaks in the ECM records originate from many sources with sulphuric acid of volcanic origin being a common source, while dips in the ECM record often are correlated to high NH + 4 concentrations, which are thought to be caused e.g. by episodes of high 25 biomass burning activity (Fuhrer et al., 1996;Taylor et al., 1996). Following the approach of Rasmussen et al. (2008), characteristic sequences of peaks and/or dips in the ECM signals can be used for robust climate-independent synchronization of ice 2979 Introduction cores, although individual peaks in themselves cannot be robustly matched between the cores. In the section of NEEM below 1757 m depth, DEP data was used to support the matching. The matching was performed by several investigators in parallel using the tool "Matchmaker" (available upon request from Sune O. Rasmussen). All sections were 5 matched independently by at least two investigators or groups of investigators. No significant differences between parallel matches were observed, although the number of match points differed markedly among investigators. In the second phase, a consensus set of match points was agreed upon by the investigators, excluding points only supported by one investigator. A total of 383 match points were identified between 10 the NEEM and NGRIP1 cores with an additional 423 match points tying together the NEEM and NGRIP2 cores (of which 19 are in the zone of overlap between NGRIP1 and NGRIP2 measurements, leading to a total of 787 match points). The match points are provided in the data file accompanying this paper and a typical sequence of ECM and DEP data with match points is shown on Fig. 3. 15 During both phases of match point identification, the Matchmaker tool is used to continuously evaluate whether the proposed match is glaciologically realistic, i.e. compatible with reasonable assumptions about the flow and accumulation regimes of the two core sites. As outlined in Rasmussen et al. (2008), the most useful diagnostics variable is the ratio r ab = (D a − D b )/(d a − d b ) of annual layer thicknesses between neighbouring 20 match points a and b found at depths D a and D b in the NEEM core and d a and d b in the NGRIP core. The values of r ab depend on both the relative accumulation rates of the cores and the amount of thinning induced by ice flow, and are expected to change only slowly (especially within a certain climatic state) as long as the spacing between a and b is large enough to eliminate the effect of short term accumulation variability, surface 25 redistribution, and uncertainties in the ECM depth scale (Andersen et al., 2006a).
The match points and the annual layer thickness ratios r ab are shown on Fig. 4  dust levels make the ice alkaline, thus muting the ECM acidity record (Ruth et al., 2003). As seen in Fig. 4, the higher NEEM accumulation leads to thicker annual layers in NEEM than NGRIP down to a depth of about 700 m, after which the increased flowinduced thinning at NEEM leads to NGRIP annual layers being thicker by a factor of 1.5-2 during most of the Glacial.

5
The (NGRIP depth, NEEM depth)-relation has a small kink at the 2087.698 m match point and clear kinks at the 2166.449 m and 2195.630 m match points. The first two kinks probably represent the boundaries between sections that have undergone different amounts of strain due to complex ice flow patterns that lead to overturned folds deeper in the ice (NEEM community members, 2013), but we believe that the ice 10 stratigraphy is still undisturbed. Below the 2195.630 m match point, the match between the cores is less robust, and there is a significant risk that the ice below this depth is stratigraphically disturbed. However, we do not have enough data of significant resolution at this point to shed light on this issue, and we therefore transfer GICC05modelext to NEEM down to a depth of 2203.597 m (or 108.2 ka b2k) below which the ice core 15 is known to be stratigraphically compromised, noting that the section below approx. 2150 m (93.6 ka b2k) is less certain and that the section below the 2195.630 m match point (107.0 ka b2k) must be considered tentative.

The relation between tephra horizons and the ECM synchronization
As described in Sect. 2.4 and Table 1, fourteen tephra horizons common to the NGRIP 20 and NEEM cores have been located. Of these, thirteen are fully consistent with the ECM-based synchronization.
The tephra particles are found in ice samples with a typical length of 15-20 cm, but their precise location can most often not be determined as the tephra layers are not visible. In contrast, the ECM record has a much finer depth resolution and precision. 25 More than half of the tephra horizons sit on top of (or very close to) one or more ECMbased match points, in which case only the more precise ECM-based match points are used for the time scale transfer. Five tephra horizons are located in areas without 2981 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | adjacent ECM-based match points and are included in the time scale transfer data set using the middle of the tephra sampling intervals (see Table 1 for details). While thirteen out of fourteen tephra horizons are consistent with the ECM-based synchronization, one tephra horizon (no. 12 in Table 1) is not. It is part of a sequence of three closely spaced horizons (no. 11-13) found around the onset of GI-8 at around 5 38 ka b2k in an area where ice samples for tephra identification have been cut continuously. The major element results in Fig. 2 show good correspondence between the tephras, and specifically for horizon 12, the layer at NEEM 1756.90-1757.10 m shows very close affinity to the tephra deposits identified in two adjacent NGRIP samples (74 shards in 2065.45-2065.65 m and 700 shards in 2065.65-2065.80 m). We believe that the two finds are related to the same volcanic event, most likely from Grimsvötn in Iceland. We find it likely that either the tephra horizon sits right in the cut that separates the two samples (2065.65 m) or that dispersion of the tephra shards has taken place, either due to processes related to the transport from the source or due to wind-driven reworking of the snow. With the exception of a few outliers in the upper sample, the ma- 15 jor element data are tightly clustered reflecting a distinct homogenous population and the geochemical signatures plot within the same envelope as horizon no. 11. Thus, two geochemically very similar horizons are present and based on the tephra analysis alone both horizons would be regarded as well-defined volcanic reference horizons.
However, from a glaciological point of view, the three tephra horizons present an in-20 consistency: the relative spacing of the horizons is different between cores and horizon 12 is offset by 0.5-1 m relative to ECM-based synchronization. If all three tephra horizons represent synchronous horizons, the NGRIP/NEEM accumulation ratio must have dropped sharply by 30 % for around a century and then increased to 120 % of the normal value for another 1-2 centuries in order to generate the apparent depth offset. Such 25 a scenario is highly unlikely given the slow variations in accumulation ratio that are usually observed. The high degree of similarity between the ECM records and the δ 18 O curves across the GI-8 onset support the conclusion that the ECM synchronization is robust and that tephra horizon no. 12 is unlikely to represent a time-synchronous event, and is therefore excluded from the time scale transfer match point data set. It is possible that the same volcanic source has produced several similar tephra layers within this period (similar to the closer spaced samples 14a and 14b in NGRIP), although subtle geochemical variations often (but not always) allow discrimination be-5 tween such events . However, as both cores have been continuously sampled across this interval (2050-2070 m NGRIP and 1740-1760 m NEEM), this explanation would require that one eruption only led to detectable tephra deposition at the NEEM site while the other event has only been registered at the NGRIP site. Although indeed possible due to the variable spatial and temporal distributions of 10 volcanic products found in Greenland ice (Coulter et al., 2012), we find this explanation somewhat speculative, but fail to be able to identify obvious alternative explanations.
Regardless of the cause, further investigations such as re-sampling of the tephra, ongoing investigations of the tephra content in GRIP and DYE-3 cores during this period, and re-measurement of the ECM record will hopefully reveal how tephra layers 15 across this section correlate.

Transfer of the GICC05modelext time scale from NGRIP to NEEM
The transfer of the time scale is based on the ECM and DEP-based match points and the five tephra horizons located away from ECM-based match points. Using the match points' NGRIP depths and the NGRIP GICC05modelext (depth, age)-relationship,

25
-Differences between the shape of peaks and inaccuracies in the depth registration of the ECM data set introduce synchronization uncertainty on the order of 2983 centimetres. The estimated synchronization uncertainty was estimated to 10 cm (1σ) by Rasmussen et al. (2008), and here we tentatively estimate its magnitude by calculating the effect on the (NEEM depth, NGRIP depth)-relation of removing every second match point. The results support the estimated synchronization uncertainty of 10 cm (1σ), leading to time scale transfer uncertainties ranging from 5 a few years to a maximum of a few decades at the deepest part of the record.
-Although we believe the set of match points to be robust, there is a risk that some sections have been erroneously matched up, leading to a larger systematic depth offset. This concern is particularly relevant for the NEEM ice below 2150 m (93.6 ka b2k).

10
As the MCE is typically 2 orders of magnitude larger than the matching uncertainty (when assuming no large systematic errors), we report GICC05modelext-NEEM-1 ages with the MCE uncertainty estimates only, but stress that observed phasing differences of up to a decade at the match point depths could be artefacts from the time-scale transfer. 15 As seen on Fig. 4, there is a section across parts of the brittle zone (845-1026 m NEEM depth) with only few match points and a section of alkaline ice of GS-2 (1511-1595 m NEEM depth) with no match points at all. In addition, many shorter periods (mainly in the stadials) do not contain match points. Even in the absence of match points we provide for convenience an interpolated (depth, age)-relation for each 0.55 m 20 depth increment, corresponding to the "bag" units used when cutting and packing the ice core. NEEM bag depths are first translated to NGRIP depths by linear interpolation using the (NEEM depth, NGRIP depth)-relation established by match points and tephra horizons. By using linear interpolation we assume a slowly varying ratio of annual layer thicknesses between the cores, which appears to be a reasonable assumption both 25 from our understanding of accumulation variability and from the smoothness of the (NEEM depth, NGRIP depth)-curve on Fig. 4. From the interpolated NGRIP depths, the ages are obtained from the NGRIP GICC05modelext (depth, age)-relation.

2984
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Note that obtaining the bag ages directly from the 792 NEEM (depth, age) match points by interpolating the age linearly between match points gives a different result. This approach corresponds to assuming constant annual layer thicknesses between the NEEM match points, which, especially in sections with match points in interstadials only, is an unrealistic assumption given the large variations in accumulation rates 5 across stadial-interstadial boundaries.

The precision of the interpolated time scale
As described in Sect. 3.2, the NEEM and NGRIP ice cores are precisely aligned at the match points and the accuracy of the NEEM time scale at the match points is thus essentially quantified by the MCE of GICC05. However, the interpolation applied to pro-10 vide a continuous NEEM time scale may introduce artificial offsets between the NEEM records and other records on the GICC05 time scale. These offsets are generally much smaller than the MCE and mainly matter when discussing the relative timing between events as observed in records from different cores.
Using linear interpolation between match point depths certainly is an approximation 15 and it is not straightforward to evaluate its uncertainty. However, the smooth shape and small curvature of the (NEEM depth, NGRIP depth)-curve seen on Fig. 4 illustrates that it is hard to justify more complex interpolation schemes. To test the influence of different interpolation schemes, interpolation was also performed using cubic spline interpolation. The change in interpolation method results in differences in the interpo- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | related to ice divide migration or lee-side effects near the GRIP and GISP2 drill sites, it cannot be ruled out that similar issues exist for the NGRIP-NEEM synchronization. Based on this and on visual inspection of the (NEEM depth, NGRIP depth)-curve, we estimate that the likely maximum depth offset caused by interpolation across the GS-2 match point gap is 1 m, corresponding to around 80 yr. However, the unusually strong 5 accumulation variability across GS-2 could potentially introduce an interpolation-based offset up to a few times larger than this. Similarly, across the ∼ 5 kyr long GS-18 and GS-19 sections where the synchronization is based on only two match points within GI-18, an interpolation-based offset of 1 m (here corresponding to around 200 yr) is conceivable.

Estimating the thinning function
Past accumulation rates can be inferred from the observed annual layer thicknesses in an ice core if the strain history at the drilling site is known. A thinning function, reflecting the accumulated vertical strain along the NEEM core, has been calculated using a Dansgaard-Johnsen (D-J) model (Dansgaard and Johnsen, 1969).
As discussed by NEEM community members (2013), the stratigraphy of the NEEM core has been compromised below 2203.6 m and, therefore, we do not attempt to model the strain in the deepest part of the core. A D-J model including basal melt and sliding (Buchardt, 2009) has been fitted to the (depth, age) relation of Sect. 3.2. The basal melt and sliding have no straightforward physical interpretation given the dis-20 turbed and folded stratigraphy but has been included in the model in order to properly represent the flow regime of the NEEM area as reflected in the upper 2203.6 m of core. The ice flow model is driven by a dynamical accumulation model, where the accumulation rate is calculated from the dated δ 18 O record through an exponential relation (Johnsen et al., 1995 ing a Monte Carlo analysis constrained by the observed (depth, age)-relation. Though match points are identified back to 108 ka b2k, only depth-age horizons younger than 85 ka b2k are used to constrain the flow parameters due to the accumulation ratio kinks discussed in Sect. 3.2 and shown in Fig. 4, which cannot be adequately described using a D-J model. As the ice is disturbed and folded below 2203.6 m, the layers immedi- sible with this simple model to obtain a good correspondence between modelled and observed depth-age horizons in the Holocene. The modelled depths are too shallow for match points younger than ∼ 8 ka and too deep for Holocene points older than this.
To further investigate this discrepancy, we look at the average annual layer thicknesses between the match points presented in Sect. 3.1. These are shown in blue in Fig. 5 15 together with the annual layer thicknesses calculated for a constant Holocene accumulation rate of 0.221 m yr −1 using a thinning function calculated from the Monte Carlo estimated parameters (Fig. 5, red curve). A prominent kink is seen in the observed annual layer thicknesses around 1200 m depth corresponding to an age of approximately 7.9 ka b2k. This kink cannot be simulated with the simple D-J model, so in order to 20 obtain a reasonable strain history for the ice above the kink, a Monte Carlo analysis constrained only by match points younger than 7.9 ka b2k is carried out. This leads to a better agreement between observed and modelled annual layer thicknesses for the period back to 7.9 ka b2k (Fig. 5, green curve). The NGRIP and NEEM annual layer thicknesses show a parallel behaviour with Early

25
Holocene layers being relatively thin compared to what is expected from a D-J model. Also, annual layer thicknesses from GISP2 seem to be somewhat thinner than expected from the model of Alley et al. (1993) between 11.7 and 9 ka b2k, whereas no such Introduction shift around 10-8 ka b2k is found in the Holocene annual layer thicknesses from GRIP . To further pursue the idea of a different Early Holocene accumulation pattern, we compare observed mean annual layer thicknesses between match points with layer thicknesses derived by the fitted D-J model (black curve in Fig. 6). A large discrepancy 5 is seen in the Early Holocene, across GI-1, and during the late Glacial. If the δ 18 Oinferred accumulation rates are reduced by 30 % during the Glacial and the reduction in accumulation rate is linearly phased out from 30 % at 11.7 ka b2k to zero at 7.9 ka b2k, the match between the Monte-Carlo-tuned D-J model and observed annual layer thickness improves dramatically (green curve in Fig. 6). 10 The origin of the changed accumulation-δ 18 O relationship in the Early Holocene is unclear, but possible explanations include (i) ice ridge movements causing the need for the flow parameters in the D-J model to be adjusted, (ii) circulation and accumulation pattern changes caused by changes in Greenland ice sheet height and geometry and the gradually diminishing effects of the Laurentide ice sheet, and (iii) changes in δ 18 O 15 caused by elevation changes and upstream effects.
The reconstructed accumulation history with the above-mentioned forced accumulation reduction in the Glacial and Early Holocene is presented in Fig. 9.

Establishing the gas chronology by estimating ∆age
Air bubbles are trapped at the firn-ice transition, leading to an ice age-gas age offset 20 that is commonly referred to as ∆age (Schwander and Stauffer, 1984). To obtain a gas chronology for the NEEM core we modelled the evolution of ∆age in the past using a coupled firn densification-heat diffusion model, with additional gas phase data to constrain the reconstruction (Fig. 7). We use δ 15 N of atmospheric nitrogen (N 2 ) where available (Guillevic et al., 2012) as well as methane (CH 4 ) data obtained through a com- 25 bination of discrete sampling (Mitchell et al., 2011;Rosen et al., 2013) and continuous flow analysis Schupbach et al., 2009;Chappellaz et al., 2013) The δ 15 N provides a strong constraint on both the timing and magnitude of abrupt temperature changes over Greenland through the imprint of thermal isotopic fractionation (Leuenberger et al., 1999;Severinghaus et al., 1998) and on past firn column thickness through the imprint of gravitational isotopic fractionation (Sowers et al., 1992). CH 4 variations are in phase with Greenland temperature, with CH 4 lagging temperature 5 by 0-70 yr Landais et al., 2004;Severinghaus and Brook, 1999;Vallelonga et al., 2012). Thus, we obtain additional timing constraints by assuming the midpoint of the CH 4 transition to slightly lag the midpoint in the δ 18 O transition at the abrupt onset of Greenland interstadials (GI); for simplicity we use a constant lag of 30 yr for most transitions, with the exception of GI-10 and 11 where available δ 15 N data in-10 dicate no lag. In several cases abrupt CH 4 and δ 18 O features are observed within an interstadial, providing additional timing constraints. The CH 4 tie points for GIs 5-11 are indicated by white circles in Fig. 7. For the gas chronology presented here, δ 15 N data were available for the last glacial termination (∼ 15-11 ka b2k) and for GIs 8-11 and 21; reliable CH 4 constraints were 15 available for the glacial termination and the GI-2 through GI-22. The ∆depth constraints from the gas data were converted to ∆age constraints through the NEEM ice chronology as presented in Sect. 3.2. A final constraint is the modern ∆age of 182 yr as derived from firn air measurements (Buizert et al., 2012).
Due to the many modelling uncertainties, the data-derived timing constraints must 20 be considered more accurate than the modelled ∆age. Because of the strong variability in Greenland temperature and accumulation during the glacial period, ∆age can change by hundreds of years in between the tie points, precluding the use of simple interpolation schemes. The philosophy behind our ∆age reconstruction is to use densification models to obtain a realistic, physically meaningful interpolation between the 25 data-derived tie points. Past changes in accumulation rate and surface temperature are the main drivers of ∆age variability, but both input functions are poorly known back in time. We adjust the accumulation (A) and temperature (T ) input to make the model fit the timing constraints.
2989 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | We use two different firn densification models in this study, with supporting information given in the Appendix. The first is the model by Goujon et al. (2003), which is a dynamical version of the 1-D densification model by Arnaud et al. (2000). The heat diffusion in the ice is calculated across the entire ice depth using a simplified version of the model from Ritz (1989). The model uses a time step of one year for both densifica-5 tion and heat diffusion calculation. The spatial resolution in the top 150 m is 0.25 m; it decreases with depth to 25 m at bedrock. The lock-in depth (LID), below which the air is no longer exchanging with overlaying layers, is defined as the depth where the ratio of closed to total porosity reaches 0.13. This LID definition based on the porosity of the ice is in agreement with the LID density threshold proposed by Schwander et al. (1997). 10 A constant convective zone of 2 m is used. This model is validated for present-day conditions in Greenland and Antarctica Arnaud et al., 2000).
The second model is a dynamical version of the Herron-Langway model (D-HL) (Herron and Langway, 1980). The model has 0.5 m spatial resolution along the depth axis, with the model domain extending down to 1000 m. The heat diffusion model 15 uses implicit Crank-Nicolson time stepping, a zero gradient boundary condition at 1000 m depth and an advective heat flux based on the 1-D ice flow model presented in Sect. 3.3. The densification and heat diffusion models use 1.0 and 0.2 yr time steps, respectively, with isochrones tracked downwards every fifth year. Lock-in density is set at 14 kg m −3 below the mean close-off density from Martinerie et al. (1994) as 20 suggested by Schwander et al. (1997); lock-in gas ages are calculated using the parameterization of Buizert et al. (2013). Following recent work by Hörhold et al. (2012), the model includes an empirical softening of ice that scales with the logarithm of the Ca 2+ concentration, where we used Ca 2+ data from Greenland summit (Fuhrer et al., 1993 The initial estimate for T (t) is based on the δ 18 O palaeo-thermometer, using a sensitivity of 0.562 % K −1 and a sea-water δ 18 O correction (Waelbroeck et al., 2002). The initial estimate for the accumulation is based on a combination of straincorrected annual layer thicknesses (similar to the red curve of Fig. 9), and a relationship with δ 18 O for the deeper core where the strain correction becomes unreliable: pared to an analytical precision on the order of 0.007 % . Consequently we use the D-HL model result for periods where we have sufficient timing constraints; while during two periods without constraints (108-90 and 23-17 ka b2k), output from both models is averaged. The combined ∆age curve is shown in Fig. 8 together with the ∆age constraints from CH 4 . The ∆age uncertainty σ ∆age is estimated in two ways. First, at the 20 CH 4 tie points we set σ ∆age equal to the uncertainty in that tie point (the error bars shown on Fig. 8) plus the absolute value of the model-tie point mismatch. In between the tie points we set σ ∆age to the linearly interpolated uncertainty of the two adjacent tie points plus 0.05 · ∆t, with ∆t the distance to the nearest tie point in years; after the last tie point at 88 ka b2k we keep σ ∆age constant. As a second estimate we use 25 the absolute value of the Goujon -D-HL model difference. For the glacial period (108-11.7 ka b2k) we use the larger of the above two estimates as the final ∆age uncertainty (shaded area). In the Holocene we tentatively set σ ∆age at 9 ka b2k to 30 yr, and again use linear interpolation to find σ ∆age at other times.

4 Summary and perspectives
The GICC05 time scale has been applied to the NEEM ice core by transfer from NGRIP using 787 horizons identified in the ECM and DEP records. Tephra layers confirm -with a single exception -this synchronization and add additional match points for the time scale transfer. The gas record has been dated by determining ∆age from firn modelling 5 with constraints from δ 15 N, CH 4 and modern day ∆age values, while the flow-induced strain history is determined from an ice flow model. The accumulation reconstructions from the ice flow and firn models show significant differences, in particular during the stadials in the middle part of the Glacial (50-25 ka b2k). The firn model reconstruction can be improved when δ 15 N data becomes available, while the model based results can be improved with improved ice flow models and additional knowledge on the evolution or the ice sheet. In addition, steps should be taken to treat ice flow and firn processes in an integrated way with the use of common accumulation and temperature input profiles. The presented time scale is intended as a first chronology for the NEEM ice core 15 facilitating the analysis of existing and forthcoming data sets from the NEEM core. When more tephra horizons and high-resolution CFA data from NEEM become available, a more detailed synchronization will be possible, hopefully also providing more match points within stadial periods and thereby increasing the precision of the time scale transfer. Future investigations will also show if NEEM data can be used to refine 20 the annual layer identification in ice cores from Greenland, e.g. in the period 10.3-7.9 ka b2k.

Data access
With the final version of this paper, the following data sets are made available at www. iceandclimate.dk/data and WDC Palaeo: Introduction -ECM data sets from NGRIP1 and NEEM in 1 cm resolution.
-NEEM DEP data from 1757 m depth downwards and NGRIP DEP data.
-The match points (including the five tephra horizons) used for the time scale transfer.
-The NEEM (depth, age)-relation in 0.55 m ("bag") resolution for ice and gas.

5
-∆age CH 4 tie points and the δ 15 N data used.
-Modelled total strain and accumulation in 0.55 m resolution.

Appendix
10 Figure A1 shows the temperature T and accumulation rates A used in the ∆age modelling. The densification models are essentially run as inverse models, where one is looking for the A and T that optimize the fit to the CH 4 tie points and δ 15 N data. The models start from an initial guess of A and T , given by the grey curves. A description of the initial curves is given in Sect.  Fig. 2) and based solely on the geochemical signatures it is uncertain which one is the NGRIP correlative. However, only the correlation to NGRIP 2077.90-2078.01 m is consistent with the ECM-based match. Similarity coefficient (SC) and statistical distance (SD) comparisons for tephra horizon correlations are presented. The similarity coefficient method is from Borchardt et al. (1972) and Hunt et al. (1995) and 7 major elements were used in the comparison. The statistical distance method is from Perkins et al. (1998Perkins et al. ( , 1995 and 10 major elements were used in the comparison. Critical value for testing the statistical distance values at the 99 % confidence interval is 23.21 (10 • of freedom 1669. 10-1669.25 1915.10-1915.50 1915.50-1915.63 Transitional  Mean NEEM annual layer thicknesses between the match points presented in this work (blue dots) and predicted Holocene annual layer thicknesses calculated from a Dansgaard-Johnsen model with constant accumulation rate of 0.221 m yr −1 and flow parameters determined from a Monte Carlo analysis constrained by observed depth-age horizons back to 85 ka b2k (red line). The green line shows the predicted annual layer thicknesses from a D-J model with flow parameter values determined from a Monte Carlo analysis constrained only by depthage horizons younger than 7.9 ka b2k. The transition from the Glacial to the Holocene (11.7 ka b2k) is marked with a horizontal grey line. High-frequency variability in the observed annual layer thicknesses caused by closely spaced match points in the top part of the core has been removed.  Fig. A1. Temperature (T , upper) and accumulation (A, lower) input into the densification models: results from the Goujon (blue) and dynamical Herron-Langway (orange) models are shown together with the grey curves, which show the initial T and A estimates as described in the text.