A first chronology for the East Greenland Ice-core Project (EGRIP) over the Holocene and last glacial termination

This paper provides the first chronology for the deep ice core from the East Greenland Ice-core Project (EGRIP) over the Holocene and the late last glacial period. We rely mainly on volcanic events and common peak patterns recorded by dielectric profiling (DEP) and electrical conductivity measurement (ECM) for the synchronization between the EGRIP, North Greenland Eemian Ice Drilling (NEEM) and North Greenland Ice Core Project (NGRIP) ice cores in Greenland. We transfer the annual-layer-counted Greenland Ice Core Chronology 2005 (GICC05) from the NGRIP core to the EGRIP ice core by means of 381 match points, typically spaced less than 50 years apart. The NEEM ice core has previously been dated in a similar way and is only included to support the match-point identification. We name our EGRIP timescale GICC05-EGRIP-1. Over the uppermost 1383.84 m, we establish a depth–age relationship dating back to 14 967 years b2k (years before the year 2000 CE). Tephra horizons provide an independent validation of our match points. In addition, we compare the ratio of the annual layer thickness between ice cores in between the match points to assess our results in view of the different ice-flow patterns and accumulation regimes of the different periods and geographical regions. For the next years, this initial timescale will be the basis for climatic reconstructions from EGRIP high-resolution proxy data sets, e.g. stable water isotopes, chemical impurity or dust records.

1 Introduction

The dating of an ice core establishes the depth–age relationship to derive a chronology of past climatic conditions from the measured proxy parameters. The proxy parameters reflect past atmospheric conditions and biogeochemical events along the core. Concerning the ice sheet as a whole, the depth–age relation is needed to map the ice sheet's internal architecture to interpret and understand the climatic evolution and the behaviour of ice streams (MacGregor et al., 2015). This is a particular focus of the East Greenland Ice-core Project (EGRIP). The drill site has been chosen close to the onset of the Northeast Greenland Ice Stream (NEGIS) (see Fig. 1), which is the largest ice stream of the Greenland ice sheet (Joughin et al., 2010, 2018). A main objective of the EGRIP project is to study the dynamics of the ice flow in the NEGIS ice stream by analysing the ice core's rheology and its relation to the deformation of the ice.

Figure 1 Locations of deep ice-core drill sites: EGRIP, NEEM, NGRIP, GRIP, GISP2, DYE-3 and Camp Century in Greenland, and close-up of the EGRIP drill site inside the Northeast Greenland Ice Stream (NEGIS). Colours show surface flow velocities from satellite data (Joughin et al., 2018).

The main objective of this work is to facilitate analysis of the data from the core by transferring the GICC05 timescale – which has already been transferred from the North Greenland Ice Core Project (NGRIP) to the Greenland Ice Core Project (GRIP), Greenland Ice Sheet Project 2 (GISP2) and North Greenland Eemian Ice Drilling (NEEM) deep ice cores (Seierstad et al., 2014; Rasmussen et al., 2013) – to EGRIP by aligning features in the DEP and electrical conductivity measurement (ECM) data sets of the EGRIP, NGRIP and NEEM ice cores. We establish a timescale for the time period of the Holocene and the last glacial termination. We stop specifically at 15 kyr b2k, as the density of the match points in the period of Holocene, GS-1 (often called the Younger Dryas) and GI-1 (or the Bølling–Allerød) is much higher than in the glacial period (and especially in the Last Glacial Maximum; LGM). We present this timescale to enable climate studies while work on a revised layer-counted timescale is ongoing.

The GICC05modelext timescale was transferred from NGRIP to the NEEM ice core by matching 787 match points of mainly volcanic origin identified in the ECM and dielectric profiling (DEP) records and – where available – verified by tephra horizons (Rasmussen et al., 2013). To apply this approach to the EGRIP core, we have profiled the upper 1383.84 m of the EGRIP core using ECM and DEP in the field during the 2017, 2018 and 2019 field seasons. We rely mainly on volcanic events as reflected in the common peak pattern in the DEP and ECM records for the synchronization between the EGRIP ice core and the NGRIP1 and NGRIP2 cores. The NEEM ice core is included in order to support match-point identification, while the GICC05 ages are transferred from NGRIP to EGRIP. Three identified tephra horizons independently verify the correct match of the ice cores.

2 Data and methods

2.1 GICC05

The annual-layer-counted Greenland Ice Core Chronology 2005 (GICC05) is derived from measurements of stable water isotopes in the DYE-3, GRIP and NGRIP (see Fig. 1) ice cores for the period back to 7.9 kyr b2k (Vinther et al., 2006) and high-resolution measurements of chemical impurities, conductivity of the ice, and visual stratigraphy from the GRIP and NGRIP ice cores for the period between 7.9 and 14.7 kyr b2k (Rasmussen et al., 2006). For the period from 14.7 to 42 kyr b2k, the dating of the cores is based on annual layer counting in the visual stratigraphy, the electrical conductivity profiles and a set of chemical impurity data (Andersen et al., 2006). The timescales are compared to timescales of different other climate archives at suitable tie points, e.g. marine sediment cores (Svensson et al., 2006). For the NGRIP core, the GICC05 timescale has been extended even further into the glacial period, back to 60 kyr b2k by annual layer counting (Svensson et al., 2008) and ice-flow modelling (Wolff et al., 2010). For the older parts (Wolff et al., 2010), the NGRIP ss09sea06bm model timescale, shifted to younger ages by 705 years, has been spliced onto the end of the GICC05 timescale, thereby forming the so-called GICC05modelext chronology. The GICC05modelext was also applied to the central Greenland GRIP and GISP2 cores by more than 900 marker points and verification with 24 tephra horizons (Seierstad et al., 2014). In summary, the GICC05modelext timescale is the consistent reference frame for the entirety of Greenland deep cores.

Previous studies assessed the differences between independent timescales of Holocene palaeoclimate records. Adolphi and Muscheler (2016) indicated that the GICC05 counting error underestimates the total uncertainty in some parts of the Holocene based on the comparison between the radiocarbon dating calibration curve (Reimer et al., 2013, IntCal13) and (Svensson et al., 2008, GICC05), and the work was extended in Adolphi et al. (2018). The objective of this work, however, is to extend GICC05 to the EGRIP core to allow parallel analysis of the records, and we thus refrain from a further discussion of the absolute accuracy of GICC05 here.

2.2 Ice-core data sets over the Holocene and last glacial termination

2.2.1 EGRIP

Here, we processed and analysed new DEP and ECM records and selected cryptotephra layers in the uppermost 1383.84 m of the EGRIP ice core. At the start of the drilling operation in 2016, the drilling site was located at 75^∘38^′ N and 35^∘60^′ W (see Fig. 1). The average annual accumulation rate is about 100 kg m⁻² yr⁻¹ (equivalent to 0.11 m yr⁻¹ of ice) for the period of 1607–2011 as determined from a firn core close to the main EGRIP drilling site (Vallelonga et al., 2014). Radar soundings suggest the ice thickness to exceed 2550 m, and traced radar layers from the NGRIP site suggest that the drill site preserves an undisturbed climatic record of at least 51 kyr (Vallelonga et al., 2014). The camp currently moves about 51 m to the north–northeast each year (Dahl-Jensen et al., 2019). Figure 2 presents an overview of the ice-core sections we used in this study. The EGRIP brittle zone is of better quality than the brittle ice from previous Greenland ice-core projects such as NEEM and NGRIP. For the EGRIP core, Fig. 3 presents a quality index on the basis of the ratio between validated and total measured DEP and ECM sample points. This quality index falls below 0.3 between 505 m (4220 years b2k) and 1210 m (11 163 years b2k) depth, which is consistent with the brittle zone between 550 and 1250 m according to the field season reports (Dahl-Jensen et al., 2019). The quality index calculated from the earlier released NGRIP and NEEM DEP data is presented in Appendix A (“Quality index for the NGRIP and NEEM ice cores”).

Figure 2 Overview of DEP and ECM that we used for the synchronization between ice cores over the Holocene and late last glacial periods.

Figure 3 Match points between EGRIP, NEEM (blue) and NGRIP (yellow) ice cores based on the DEP and ECM data sets. The core quality index Q is derived from the validated DEP (blue) and ECM (red) data, respectively.

2.2.2 NGRIP

The GICC05modelext timescale, as discussed in detail in Sect. 2.1 above, is well established for the NGRIP ice core. To fully exploit the potential of DEP records for matching, we processed unpublished DEP data from the NGRIP1 core for the upper part (down to 1298 m), and we used the NGRIP2 (below 1298 m) that was published with the NGRIP ECM data in Rasmussen et al. (2013). The NGRIP1 and NGRIP2 cores have a depth offset of around 0.43 m between corresponding events in the overlapping section (Rasmussen et al., 2013).

2.2.3 NEEM

The firn core NEEM-2008-S1 originates from the NEEM access hole of the main core, drilled during the 2008 field season to a depth of 103 m (Gfeller et al., 2014). We used only ECM data for the matching of the upper 100 m, as DEP was not measured on the access-hole core. Below this depth, both DEP and ECM were used to transfer the GICC05 timescale from the NGRIP to the NEEM core (Rasmussen et al., 2013). The shallow and deep cores overlap, forming a continuous record.

2.3 Field measurements and data processing

2.3.1 Dielectric profiling

DEP has been introduced as a system for rapid scanning of ice cores' electric permittivity and conductivity shortly after drilling (Moore and Paren, 1987; Wilhelms et al., 1998). The permittivity and conductivity of ice and firn are determined by their respective densities and conductivities (Wilhelms, 2005). The conductivity is related mainly to acidity, salt and ammonia concentrations of ice cores (Moore et al., 1992, 1994). The dielectric stratigraphy of the EGRIP, NEEM and NGRIP cores was recorded directly during the field seasons with the DEP device described by Wilhelms et al. (1998) (Fig. 4), in the discussion below referred to as “deep-core DEP”. DEP is the first measurement within the processing line directly on site. A few minutes before scanning, the core is moved from the core storage to the DEP table. Further along the processing line, the ice core is split into the different aliquots. For all three ice cores, DEP measurements were carried out on 1.65 m long sections (Fig. 4d).

Figure 4 Schematic of the DEP instrument.

All cores were scanned with the deep-core DEP device as described in Wilhelms (1996) and Wilhelms et al. (1998), and the entire scan of a core's section was calibrated with average values of the calibration measurements. For the processing of the NGRIP and NEEM cores, we have improved some features and adapted procedures that were developed for a different DEP device (Wilhelms, 2000), not yet described elsewhere in the application for the processing of data recorded with the deep-core DEP device. The identification of peaks is sufficient for the discussion in this paper. However, the data released with this paper will also be valuable for investigations relying on well-calibrated material properties with absolute calibration, e.g. modelling of synthetic radargrams. To provide a comprehensive presentation of the basis of the transfer of the timescale here, we outline the relevant discussion to operating the DEP system, while the related discussion on precision of the measured material properties, which are of more relevance for later use of the calibrated data, is presented in Appendix B (“Calibration and corrections to the DEP data”).

Due to the drilling procedure and properties of the ice, ice cores can exhibit breaks or broken-off slices, or may in some instances (especially in the brittle zone) be fragmented (Fig. 4b). The missing pieces and free surfaces with possibly high conductivity have the potential to introduce artefacts into the DEP record. These are clearly identifiable in the permittivity record by dropping spikes. For the validation of the data, any drop in permittivity below a certain threshold (illustrated by the red line in Fig. 5) identifies a spike to be rejected, where the segment to be rejected is extended from the spike to where the signal approaches the average of the permittivity record. The automated procedure as described in Rasmussen et al. (2013) (Sect. 2.3) is much faster and more consistent in between the three different cores and has proven to be superior to any approach based on a hand-written protocol, which depends on the judgement of the operator when identifying intervals of a bad core. As the permittivity is very sensitive to bad core quality and the conductivity is much less prone to bad core quality, the outlined validation procedure leads to a robustly validated conductivity record.

Figure 5 Example of DEP data processing (removing core breaks from the raw data). Raw data are plotted in grey and validated data are in black. We define insufficient data by permittivity drops below the threshold (indicated by the red line) and reject it during validation. Permittivity is presented in panel (a) and conductivity in (b). The insert shows the details for a short section (250–260 m).

The automated validation straight forwardly leads to a definition of a core quality index $Q\in [\mathrm{0},\mathrm{1}]$ by calculating the total length of validated core sections divided by the standard DEP and ECM run lengths of 1.65 m for all three cores. For the EGRIP core, the core quality indices as derived from the DEP and ECM records, respectively, are presented in Fig. 3.

Figure 6 Example of the DEP data gap length and relative precision of ECM vs. DEP depth assignment. The DEP measurement with a wrong depth in grey and after shifting depth in black. There is a gap of 35 cm (1300.4–1300.75 m) between two DEP measurements. Blue arrows show the early measurements before the starting position of the core. Corresponding ECM data are in red.

When processing the EGRIP core, over a certain period, the operators erroneously did not reset the starting position of the scanning electrode of the DEP device. This is clearly identifiable in the records and resulted in recording already ahead of the core's top depth (blue arrows in Fig. 6), then taking measurements over the (correctly set) length of the section but missing the corresponding length at the bottom end of the core section. The respective core sections have been shifted accordingly, but the missing end sections, which were not recorded, cannot be recovered. The reconstructed DEP record was compared and validated against the ECM record to assign the correct depths. The section at about 1285–1385 m was corrected in this way, where in total about 8.5 m of the 100 m were not measured. Furthermore, we relied more heavily on the ECM record than on the DEP record when matching peaks within sections with known problems. Figure 6 illustrates the corrections in the interval of 1299.10–1302.05 m, where a 35 cm data gap between two scanned sections cannot be reconstructed as it was not recorded due to the wrong positioning of the electrode.

2.3.2 Electrical conductivity measurements

For the EGRIP core, we recorded ECM profiles with the technique described by Hammer (1980) directly in the field. The ECM signal is related to acidity concentrations of ice cores, even with high concentrations of neutral salt (Moore et al., 1992). NGRIP (Dahl-Jensen et al., 2002) and NEEM (Rasmussen et al., 2013) were measured using similar equipment during the respective processing campaigns in the field. For each measurement, the hand-dragged ECM instrument was moved along the depth axis of the ice-core sections' microtome-polished surface (three-bag sections, about 1.65 m long). In order to calibrate the ECM data, as described in Rasmussen et al. (2013), the ice temperature was measured for each run. In addition, to ensure the ECM quality, we repeated our measurement at least twice for each core section and checked the profiles for the best quality of the measurements. Also, the core-break positions were registered along with measurement by moving the electrodes of the ECM instrument to the respective break position after the core scan and registering the position in the data file. During the processing, these recorded break marks were used to trim off artefacts and produce the final ECM data set. Data from each day were calibrated using independent measurements of the physical dimensions of the ECM measurement setup. The first and last few millimetres of recorded data are affected by the proximity to the end of the core and were removed. Areas with dips in the signal around logged core breaks were also muted during processing. Details on the acquisition and processing of the ECM record are laid out in Appendix C (“Details on the ECM procedures”).

The ECM current, i (in µA), was converted to ice acidity (in µequiv. H⁺ kg⁻¹) by using the relationship $\left[{\mathrm{H}}^{+}\right]=\mathrm{0.045}\times i^{\mathrm{1.73}}$, as suggested by Hammer (1980). Even though conversion from current to acidity and calibration curves have been shown to be ice-core dependent, the matching and synchronization of the ice cores are independent of the absolute values of the calibrated ECM signal as they rely on recognition of similar patterns and peaks in the acidity records (Rasmussen et al., 2013).

The quality of the processed data was checked by comparing independently processed ECM data by three investigators. No major disagreements were found when comparing, and one set of data was agreed on for further use in matching.

2.3.3 Tephra horizons

The EGRIP core was continuously sampled for tephra analysis in the field to maximize the probability for detection of volcanic ash deposits, particularly invisible deposits that exhibit low concentrations and/or small grain sizes, known as cryptotephra (Davies, 2015). An aliquot of ice was prepared from the outer curved edge of each 55 cm ice-core section and subsampled at 11 cm resolution, providing approximately 30 mL of meltwater per sample. Samples from sections with either significant peaks (Fig. 4f) in the DEP and ECM or visible layers were separated and prioritized for screening, whereby centrifuged samples were evaporated onto frosted-glass microscope slides and covered in epoxy resin to enable scanning by high-magnification light microscopy, as described by Cook et al. (2018). Electron probe microanalysis (EPMA) by wavelength dispersive spectrometry (WDS) was performed to determine the major element composition of individual grains in each deposit (Hayward, 2012), and EPMA measurements were performed using a Cameca SX100 electron probe microanalyser at the Tephrochronology Analytical Unit, University of Edinburgh. This system has five wave dispersive (WD) spectrometers and was calibrated daily using internal calibration standards as described by Hayward (2012). With optimized instrument settings for the analysis of small cryptotephra grains (<20 µm diameter), EGRIP samples were analysed with either a 5 or 3 µm beam diameter. Secondary standards were analysed to capture instrument drift. The geochemical composition of each layer was compared to deposits in NGRIP and positive matches were used to establish the independent tie points in between cores listed in Table 1. Major element biplots (Fig. 7) show graphical correlations for each NGRIP–EGRIP tephra match point, and these are supported by similarity coefficient (Borchardt et al., 1972) and statistical distance (D²) (Perkins et al., 1995, 1998) tests.

(Mortensen et al., 2005)

Table 1 Geochemical matches between EGRIP and NGRIP were supported by the similarity coefficient (SC) test of Borchardt et al. (1972) and statistical distance (D²) test of Perkins et al. (1995, 1998). Here, we provide SC and D² values for major elements (normalized to 100 %) where five major elements (with > 1 % wt) were used for SC calculations for sample pairs with rhyolitic composition and seven elements were used for sample pairs with basaltic composition. Values >0.95 suggest products are from the same volcanic source. For D², seven major elements were used for the comparisons (with >0.01 % wt). The value for testing the statistical distance values at the 99 % confidence interval is 18.48 (7 degrees of freedom).

Figure 7 Element–element biplots showing geochemical matches between EGRIP and NGRIP samples with the exception of EGRIP at 894.41–894.52 m (EGRIP bag 1627), which is shown here with Mazama data from Jensen et al. (2019). Geochemical data are normalized to 100 % (anhydrous basis), and analyses with totals below 94 %wt were excluded. Panels (a–d) are for EGRIP bags 177 and 1627; panels (e–h) are for EGRIP bag 2094.

2.4 Synchronization of dielectric profiling and electrical conductivity measurement records of EGRIP, NGRIP and NEEM

Patterns in the DEP records of NGRIP, NEEM and EGRIP were initially matched by one investigator. The same cores' ECM data were matched separately and independently by three different investigators. Both matches are mainly based on clearly identifiable volcanic peaks and also synchronous patterns of other events (see Fig. 8), which not necessarily need to be of volcanic origin but are assumed to reflect synchronous events. Based on these independent matches, the four investigators identified consistent and reliable common patterns that are represented in the ECM and/or the DEP records from NGRIP and at least one of the other ice cores. For the confirmation of match points, all records of all three cores were loaded into the Matchmaker tool (Rasmussen et al., 2013) and assessed jointly by all four investigators in the different display options featured by the software. The Matchmaker tool allows easy identification of wrong match points via interactive plots and online evaluation of the match. To validate match points, we plot the depths of the common match points D_i (in EGRIP or NEEM) against d_i (NGRIP). The slope of each of these (depth, depth) curves is the annual layer thickness ratio of the two cores, $r_i=\frac{D_{i+\mathrm{1}}-D_i}{d_{i+\mathrm{1}}-d_i}$. Points which deviate from the (depth, depth) curve or create jumps in r are easily recognized and checked again. We only expect significant abrupt changes in r at times where the climate (and thus the relative accumulation rates) shifts due to changes in climate conditions (Rasmussen et al., 2006, 2013; Seierstad et al., 2014; Winski et al., 2019), while the different ice-flow patterns at the cores' sites only lead to slow changes in r.

Figure 8 Example of the ECM and DEP data match between the NGRIP2 (a), NEEM (b) and EGRIP (c) cores. The match points are marked by grey bands.

Short-term accumulation variability due to both climatic factors and wind-driven redistribution of snow on the surface can lead to relatively large variations in the ratio of layer thicknesses between different cores, especially when match points are only a few years apart. To reduce short-term accumulation-rate variability in the final timescale, we re-evaluated intervals with large variability in annual-layer-thickness ratios and removed match points that were too closely spaced. The final minimum distance between match points is 0.22 m (1206.45–1206.67 m), corresponding to around 3 years. Overall, the match points are reasonably evenly distributed throughout the entire ice core, and the maximum distance between neighbouring match points is 26.6 m (490.06–516.67 m), corresponding to a time interval of 224 years.

2.5 Transfer of the GICC05 timescale to the EGRIP ice core

The procedure of transferring the timescale is similar to the approach described in Rasmussen et al. (2013). Note that we hereby assume that the ratio of annual layer thicknesses is constant between the match points of EGRIP and NGRIP. For each 0.55 m EGRIP depth segment (the so-called bag), we obtain the equivalent NGRIP depth by linear interpolation between the depths of the match points D_i in EGRIP and d_i in NGRIP. We then assign a GICC05 age from the annually resolved GICC05 timescale for NGRIP (Vinther et al., 2006; Rasmussen et al., 2006).

The EGRIP timescale inherits the maximum counting error (MCE) from the GICC05 timescale. Our match covers the time period back to 14 967 years b2k where the associated MCE is 196 years. The inaccuracies in the depth registration were estimated by Rasmussen et al. (2013) to 10 cm (1σ). For the joint assignment of DEP and ECM patterns, we repeat the assessment by Rasmussen et al. (2013) for the (EGRIP depth, NGRIP depth) relation by computing the difference δ_i between each EGRIP match point and the linear interpolated depth derived from the neighbouring points ${\mathit{\delta }}_i=\left(D_{i+\mathrm{1}}-D_{i-\mathrm{1}}\right)/\left(d_{i+\mathrm{1}}-d_{i-\mathrm{1}}\right)\times \left(d_i-d_{i-\mathrm{1}}\right)+D_{i-\mathrm{1}}-D_i$, thus merging all match points between EGRIP–NGRIP1 and EGRIP–NGRIP2 into one data set of 377 of the originally 381 match points for further statistical analysis. The difference in number occurs as δ is not defined for respective start and end points of both respective sequences. The statistical analysis of δ in Appendix D yields a standard deviation of 0.043 m for the depth assignment of a match point. As the annual layer thickness typically exceeds 0.04 m in the time period considered here, one expects an additional uncertainty for the peak assignment on the order of 1 year.

Larger errors would occur in the case of erroneously matched sections as discussed in Rasmussen et al. (2013), but wrong matches are even more unlikely here than in previous work, as three instead of two cores were matched. Rasmussen et al. (2013) also point out that the risk of erroneously matched sections is particularly relevant for the older part of the core, where the discussion here covers a section with comparably plentiful match points.

2.5.1 Precision and accuracy of the timescale transfer

The central mode of the (depth, depth) differences δ as defined above follows a Gaussian normal distribution with a standard deviation of 0.043 m (see Appendix D). This demonstrates that the synchronization ties two cores together at the match points with high precision. Besides the central Gaussian normal distribution, the statistical analysis of δ identifies an overlaid second Gaussian normal distribution with a standard deviation of 0.19 m. We interpret this distribution as stemming from curvature of the (depth, depth) curve and as indicative of the average amount of detail, that each point contributes to the description of this curvature. As this difference between the actual match point and the linear interpolation between the neighbouring points depends on the variable curvature of the (depth, depth) curve, which reflects both accumulation conditions and the evolution of glaciological conditions at both coring sites, the values of δ will generally be time correlated and cannot be expected to be randomly distributed with depth. We thus estimate that the EGRIP timescale may have time-correlated uncertainties relative to NGRIP of up to a handful of years related to changes in relative accumulation variability and ice-flow conditions not captured by the match points. This uncertainty will be largest where the closest match point is relatively far away and near climatic shifts where the accumulation changed abruptly and not necessarily by the same ratio at different ice-coring sites.

Now we assess the combined uncertainties. The GICC05 timescale t(D) inherits the associated maximum counting error from GICC05, and given the analysis of δ above, we conclude that, at the match points, the timescale is precise relative to NGRIP within about 1 year (1σ). However, when we want to know the age at an arbitrary EGRIP depth, additional uncertainties apply due to the interpolation between the match points. There are two dominant sources. As discussed above, variations in relative accumulation rates and ice flow may add up to a handful of years of additional uncertainty relative to GICC05, but there is also a contribution from the choice of interpolation scheme in between the match points. The difference introduced by the choice between the most widely used linear and cubic spline interpolation schemes (Press et al., 1992) is about an order of magnitude larger than the above-mentioned random uncertainty associated with the identification of the match points (see Appendix E).

We maintain linear interpolation for the timescale transfer despite the fact that the slope of the (depth, depth) curve changes instantaneously at the match points. While changes in this slope may in reality occur on many scales due to the intermittency of precipitation, wind-driven redistribution of snow and relative changes in accumulation rates at the sites, the most significant of these changes are likely to happen at times of climate change rather than at the arbitrary depths of the match points. Considering this, one could consider an interpolation scheme where the change of curvature is distributed over the entire curve and/or concentrated at times of climate changes as derived from the proxies of the ice core. However, we believe that the advantage of obtaining a smoother (depth, depth) curve does not compare favourably to the additional assumptions needed and the added complexity of the timescale transfer.

3 Results and discussion

3.1 Synchronization of the EGRIP, NEEM and NGRIP cores

A total of 257 match points between the EGRIP, NEEM and NGRIP1 ice cores and 124 match points between EGRIP, NEEM and NGRIP2 (total of 381 match points) with an additional three tephra horizons were identified. Figure 8 shows an example section of ECM/DEP records matched between ice cores. The match points between the ice cores are shown in Fig. 3. In the process of combining match points from all investigators, some match points were removed due to differences in the peak shapes between DEP and ECM data or when there were too many match points very close to each other. There are fewer match points in the interval of 600–1100 m due to the brittle zone, which in particular influences the NEEM and NGRIP1 cores. The ECM and DEP do not follow each other closely in the 1245–1283 m interval because of the alkaline nature of the ice associated with stadial conditions in EGRIP. This is due to high dust levels neutralizing the acidity of the ice (Ruth et al., 2003; Rasmussen et al., 2013).

3.2 Tephra horizons identified for the chronology

Three tephra horizons have been located in EGRIP (Table 1). The locations of these horizons were consistent with the DEP- and ECM-based synchronization. The tephra horizons thus provide an independent validation for our match points. In addition, ongoing tephra investigations will likely provide additional points for synchronization between ice cores in intervals without DEP and ECM match points.

3.3 GICC05-EGRIP-1

As described in Sect. 2.5, the GICC05 (depth, age) timescale was transferred from the NGRIP to the EGRIP ice core based on 381 match points. We name the relationship between depth and age for EGRIP over the Holocene and early last glacial periods GICC05-EGRIP-1 and present the average annual layer thickness between the match points in Fig. 9.

Figure 9 Average EGRIP annual layer thicknesses (dark blue line, left y axis) between the match points. The EGRIP depth–age relationship (right y axis) with match points (red dots) and the tephra horizons (black squares).

Figure 10 The EGRIP–NGRIP (orange) and EGRIP–NEEM (blue) annual-layer thickness ratio (left axis) calculated between neighbouring match points.

We synchronized the records of the ice cores back to 14.96 kyr b2k which corresponds to EGRIP depth of 1383.84 m, NEEM depth of 1493.29 m and NGRIP2 depth of 1611.98 m. Along with this publication, we release a timescale for each 0.55 m section (“bag”). For each EGRIP bag depth, the corresponding NGRIP depth was found by linear interpolation between the match points, and the GICC05 age was then determined from the GICC05 timescale for NGRIP. The maximal uncertainty resulting from the choice of interpolation scheme is assessed in detail (see Appendix E1) and is about 4 years. The relatively smooth (depth, depth) relation of EGRIP–NGRIP and EGRIP–NEEM (see Fig. 3) shows that the ratios of annual layer thicknesses between cores do not vary noticeably between match points. Figure 10 shows that EGRIP has thinner annual layers than both NEEM and NGRIP ice cores in the upper parts of the cores as also expected from the lower surface accumulation. Ice found in the EGRIP core originates from snow that was accumulating upstream, and accumulation rates increase upstream as the flow line approaches GRIP and NGRIP, where present-day accumulation is about twice that at EGRIP (Vallelonga et al., 2014; Riverman et al., 2019; Karlsson et al., 2020). Surprisingly, annual layer thicknesses in EGRIP remain almost constant back to 8 kyr b2k (Fig. 9), while the layer thicknesses in large parts of the Holocene part of the NGRIP and NEEM cores thin linearly due to ice flow. We believe that it is a coincidence that the combined effects of the increasing upstream accumulation and flow-induced thinning at EGRIP balance out for the last 8 kyr. Despite the lower accumulation at EGRIP, annual layers in EGRIP eventually get thicker than the annual layers in the NEEM and NGRIP ice cores. Below an EGRIP depth of around 700 m, annual layers in EGRIP are thicker than the layers from the same period in the NEEM core, and similarly below 1000 m, EGRIP annual layers are thicker than those in NGRIP (Fig. 10).

There are some gaps in the EGRIP ice-core record due to the brittle zone. However, the smoothness of the depth vs. depth plot in Fig. 3 and the annual layer thickness ratio in Fig. 10 robustly support our timescale based on the match points.

3.4 Error analysis for the timescale transfer

We statistically treated the “leave-one-out analysis” (δ), which was already used in a previous paper (Rasmussen et al., 2013) to estimate the uncertainty in match points. Binning the values revealed a robust mode that refers the statistical error of the match-point assignment. The tails of the distribution relate to more curved parts of the (depth, depth) curve defined by the match points. Cubic spline interpolation provides a similar estimate. However, when matching depths, we do not know the true curve nor how big the influence from the interpolation scheme is. We demonstrated that the two fundamentally different interpolation schemes of linear and cubic spline interpolation give a consistent result. The δ values give qualitatively a similar pattern to the high-resolution difference of the interpolation schemes and have demonstrated their suitability to estimate the uncertainty of the interpolation scheme. They deviate more than the statistical uncertainty of the match-point assignment, but this is less than 10-fold of the statistical error for peak identification, and for the timescale here it introduces an error of up to 4 years, while the inherited maximum counting error of the timescale increases from about 1 to 200 years over the matched record.

4 Conclusions

We have established the initial chronology for the EGRIP deep ice core in Greenland which encompasses the Holocene and late glacial periods. We have established the depth–age relation for the upper ∼1383.84 m of the core back to approximately 14.96 kyr b2k based on the GICC05 timescale and labelled it GICC05-EGRIP-1. After field measurements and processing of the ice-core data, we relied on the DEP and ECM records for the synchronization, using 381 match points between EGRIP, NEEM and NGRIP ice cores. The identification of tephra match points between the EGRIP and NGRIP cores provides an independent tool for validating this synchronization. We used the ratio of annual layer thickness between ice cores as a tool to evaluate our match points. This first timescale can help to interpret, design sampling strategies for and improve the understanding of the forthcoming EGRIP data sets.

Appendix A: Quality index for the NGRIP and NEEM ice cores

For the NEEM and NGRIP ice cores, we calculated similar quality indices to those provided for EGRIP above. They are presented together in Fig. A1.

Figure A1 Quality indices for the EGRIP, NGRIP and NEEM ice cores.

Appendix B: Calibration and corrections to the DEP data

As a correction in the few percent range, we correct the offset, introduced by the changing stray admittance due to the varying cable geometry due to its movement during the measurement, by subtracting the course of free-air measurements from the respective measurement of a core section along the DEP device when processing the data. This is possible as the stray admittance is connected in electric series and is thus additive. This measure reduces coherent noise on the record, which is, e.g. important when using the records for spectral analysis.

For the processing of the two NGRIP cores, reproducibility was ensured by laying the cables out to move freely in the same way for all measurements in between two recorded free-air measurements. This was improved for the NEEM and EGRIP processing by placing the cables into cable channels that enforce repeatable deformation.

For the calibration of the DEP device, free-air measurements without ice were recorded frequently, usually at least twice daily before processing started and finished. The slight capacitance and conductance variation on the order of less than 4 fF and 500 pS, thus corresponding to relative permittivity changes of $\mathrm{4}\mathrm{fF}/\mathrm{63}\mathrm{fF}=\mathrm{0.06}$ and conductivity changes of $\left(\mathrm{500}\mathrm{pS}\right)/\left(\mathrm{63}\mathrm{fF}\right)\times \left(\mathrm{8.8542}\mathrm{pF}\right)=\mathrm{70}\mathrm{nS}$, along the DEP device is due to the unavoidable deformation of the cables (Fig. 4a) when moving the scanning electrode along the device (Fig. 4e). Compared to the properties of pure glacier ice (refer to Fig. 5), these variations are on the order of 2 % for the permittivity and 5 ‰ for the conductivity. Additionally, an offset of few nS of residual conductance may remain even after performing the correction routines of the LCR meter (inductance L, capacitance C, resistance R) bridge (Fig. 4c).

As a further correction in the few percent range, we developed procedures to determine the true absolute free-air capacitance of the DEP capacitor. This is relevant when determining calibrated absolute values of the material properties: permittivity and conductivity. In the original publications (Wilhelms, 1996; Wilhelms et al., 1998), the proper calibration of the device is cross checked by comparison with the theoretical capacitance value of 63.4 fF, where, for a precisely adjusted DEP bench, the free-air capacitance coincides within less than 2 fF. These small deviations from the theoretical value might well be due to mechanical tolerances like the electrode length in the range of a few tenths of a millimetre. However, for a slightly differently adjusted device (e.g. one with slightly more clearance to the core), the deviation from the ideal value of the free-air capacitance might be a few fF more. Besides the calibration uncertainty of the LCR meter, it might also include a component of cable stray capacitance, which is not identifiable in the LCR meter's automated correction procedure.

To even proceed from the correction of cable stray admittance variation along the course of the device by simple subtraction towards absolute precision, one needs to know the free-air conductance, which is expected to vanish for the empty device, and the true free-air capacitance of the capacitor, which needs a special approach as it cannot be measured directly and separately from other interfering capacitances.

In parallel to the NGRIP project, Wilhelms (2000) developed a calibration procedure for a custom DEP device with fixed electrodes, which was optimized for firn studies to establish the density-and-conductivity mixed permittivity (DECOMP) model (Wilhelms, 2005). For the EGRIP core processing, we transferred the principle of introducing concentric metal tubes to the deep-core DEP bench and upgraded it with a rack to mount to move the tubes of different diameters along with DEP electrodes and record the capacitance along the DEP device. The tube in the electric field increases the capacitance of the arrangement and Wilhelms (2000) derives the theory to calculate the effective relative permittivity of the setup. For the calibration, tubes with radii in approximately 10 mm increments between 0 and 70 mm represented effective permittivity standards ($\left|\mathit{\epsilon }\right|$) between 1 and 4 (Fig. B1). The result is a calibration curve which holds for the calculation of a consistent free-air capacitance for the correction of the DEP measurement of the EGRIP core. The free-air capacitance is the proportionality factor of the measured capacitance and the effective permittivity, i.e. the slope of the graph in Fig. B1. From this analysis, one derives the true free-air capacitance of the deep-core DEP device as adjusted when it was assembled.

Figure B1 Calibration curve of the measurements with artificial tubes and free-air measurement for DEP device.

Now, we know the true free-air capacitance (as determined from the just-outlined calibration) together with the true (vanishing) conductance for an empty device. The free-air empty capacitor measurements thus determine the (additive) stray conductance and capacitance, which is the difference of the measurement and the reference values. The derived stray conductance and capacitance profile along the course of the DEP device are then subtracted from the actual ice-core measurements. The calculation of both material properties involves a division with the free-air capacitance, where introduction of further errors is minimized by using the precise value from the calibration procedure as outlined above.

Precise permittivity and resulting precise free-air capacitance values mainly make a difference when, e.g. computing the wave propagation speed of radar waves while modelling synthetic radargrams (Eisen et al., 2006). For NGRIP and NEEM, we did not yet perform the calibration procedure. From assigning the two-way travel time of prominent radar reflectors to volcanic spikes in the core, one can determine the radar wave propagation speed and calculate the ice's permittivity. As the NGRIP and NEEM data sets will also be used for comparison with radar surveys later, we determined the free-air capacitance by averaging the measured capacitance over deep-core sections and dividing by the expected permittivity of ice of 3.15±0.1 which computes the free-air capacitance with 3 % relative error, which is only about a 2 fF absolute error for the free-air capacitance.

To sum up, all material properties' data sets we derived here are only subjected to a few percent absolute error. While the EGRIP record is calibrated independently, for the NGRIP and the NEEM cores, the permittivity, as determined from radar wave propagation, was used. The latter is a very minor restriction as this cannot be checked independently but is not of practical relevance.

Due to the varying temperature in processing area throughout the field seasons, the core was not processed at a consistent temperature and we do not have the temperature readings avail to provide consistently harmonized conductivity data. The missing temperature correction does not affect the use of conductivity peaks for synchronization purposes in between ice cores, which is relevant for the discussion here. When, e.g. deriving radar wave absorption coefficients from the presented conductivity record, one would have to be very cautious and have this limitation of the data in mind.

Appendix C: Details on the ECM procedures

The depth scale of the ECM profile was assigned based on the recorded movement of the electrodes interpolated between the logged top and bottom depth of individual ice-core sections. To investigate the quality of the depth assignment, a bag mark position analysis was carried out on the section below the brittle zone in EGRIP (∼1160–1760 m). Only ice-core sections with an undamaged core of 1.65 m length were included in the analysis. Each 1.65 m section contains the equivalent of three 0.55 m bags, and the true depths of the bag interfaces (“bag marks”), separating the first and second bags and the second and third bags, respectively, are known. During ECM measurements, these bag marks are logged (just as the break marks) but are not used for the processing. After processing of the ECM signal, the positions of logged bag marks were interpolated onto the same depth scale as the processed ECM signal, making it possible to compare the true depth of these marks to their depths in the processed data. The distances in depth between logged and expected positions of individual bag marks were calculated for all sections included in the analysis. It was found that the depth assignments of the bag marks were almost always accurate within 20 mm, with mean distance μ=8.3 mm and standard deviation σ=7.9 mm.

Appendix D: Statistical analysis of the (depth, depth) match with linear interpolation

Each match point's depth assignment has an uncertainty due to the varying peak form, which is caused by regional deposition differences and short-term accumulation variations. This peak assignment uncertainty dominates the distribution of δ in sections where the EGRIP–NGRIP1/2 depth curve is straight. However, varying conditions at the time and place of the snow deposition or different ice-flow patterns between the ice cores cause differently evolving annual layer thickness ratios $r_i=\frac{D_{i+\mathrm{1}}-D_i}{d_{i+\mathrm{1}}-d_i}$ and recognizable curvature of the (depth, depth) curve. For the distribution function of δ, we expect a normal distribution of the peak assignment overlaid by a distribution from the curvature of the (depth, depth) curve, referred to as “the refinement” in the following.

Figure D1 Histogram of the difference data with linear interpolation ${\mathit{\delta }}_i=(D_{i+\mathrm{1}}-D_{i-\mathrm{1}})/(d_{i+\mathrm{1}}-d_{i-\mathrm{1}})(d_i-d_{i-\mathrm{1}})+D_{i-\mathrm{1}}-D_i$ for the match points between EGRIP and NGRIP1/2.

For further statistical treatment, we bin the δ values. As $\mathit{\delta }\in (-\mathrm{0.42929}\mathrm{m},\mathrm{0.385968}\mathrm{m})$, we counted the occurrence N(i) of values for the i=0…20 intervals $(-(\mathrm{0.525}+i\times \mathrm{0.05})\mathrm{m},(-\mathrm{0.475}+i\times \mathrm{0.05}\left)\mathrm{m}\right)$ in between −0.525 and 0.525 m, and display the data in the following histogram (refer to Fig. D1). The standard deviation of the counts is $\mathit{\sigma }\left(N\right)=\sqrt{N}$. Overall, n=14 bins are occupied. The distribution of δ is presented in Fig. D1.

The Gaussian normal distribution for the peak assignment in the centre seems to be overlaid by a second Gaussian normal distribution representing the refinement of the depth scale. This is indicated by the quite wide tails of the distribution. Following the just-posed assumption that the refinement distance is statistically normal distributed $N_B\left(\mathit{\delta }\right)=B\exp \left(-(\mathit{\delta }-{\mathit{\mu }}_B{)}^{\mathrm{2}}/{\mathit{\sigma }}_B^{\mathrm{2}}/\mathrm{2}\right)/\sqrt{\mathrm{2}\mathit{\pi }}/{\mathit{\sigma }}_B$ and that the uncertainty of the peak assignment in smooth intervals is normal distributed $N_A\left(\mathit{\delta }\right)=A\exp \left(-(\mathit{\delta }-{\mathit{\mu }}_A{)}^{\mathrm{2}}/{\mathit{\sigma }}_A^{\mathrm{2}}/\mathrm{2}\right)/\sqrt{\mathrm{2}\mathit{\pi }}/{\mathit{\sigma }}_A$, we χ² fitted the sum $N\left(\mathit{\delta }\right)=N_A\left(\mathit{\delta }\right)+N_B\left(\mathit{\delta }\right)$. The weight of the χ² is the counting error $\mathit{\sigma }\left(N\right)=\sqrt{N}$. The fit of the c=6 independent parameters converged and yielded a χ²=4.7 for the $f=n-c=\mathrm{14}-\mathrm{6}=\mathrm{8}$ degrees of freedom, indicating that the fitted distribution is supported by the data $\left({\int }_{\mathrm{4.7}}^{\mathrm{\infty }}{\mathit{\chi }}^{\mathrm{2}}\left(\mathrm{8}\right)=\mathrm{0.79}\right)$.

The χ² fit computes for the peak assignment distribution N_A(δ) a scaling factor $A=(\mathrm{16}\pm \mathrm{1})$ m, the shift from the centre ${\mathit{\mu }}_A=(\mathrm{0.003}\pm \mathrm{0.003})$ m and a standard deviation ${\mathit{\sigma }}_A=(\mathrm{0.043}\pm \mathrm{0.002})$ m. Similarly for the refinement distance N_B(δ): $B=(\mathrm{3}\pm \mathrm{0.7})$ m, ${\mathit{\mu }}_B=(-\mathrm{0.05}\pm \mathrm{0.04})$ m and ${\mathit{\sigma }}_B=(\mathrm{0.19}\pm \mathrm{0.04})$ m.

Solving N_A(δ_j)=N_B(δ_j) for the roots ${\mathit{\delta }}_{\mathrm{1}}=-\mathrm{0.11}$ m and δ₂=0.12 m defines the inner interval that is dominated by the peak assignment statistics. A total of 349 points are in the interval (δ₁,δ₂), and the direct statistical evaluation confirms σ_A=0.043 m (in the main paragraphs of the paper, we label this as statistical error for the peak assignment ΔD) and μ_A=0.003 m, where the skew (0.06) and the kurtosis (−0.05) are small and support normal distribution $\mathit{\delta }\in (-\mathrm{0.11},\mathrm{0.12}\mathrm{m})$. A Shapiro–Wilk test confirms normal distribution of the peak assignment errors, as W is 0.996 and the corresponding p value is 0.54.

As the refinement distance is overlaid by the peak assignment statistics in the centre of the distribution, none of the standard statistics are applicable, but the χ²-fitted N_B(δ_i) can be χ² tested for the refinement-distance dominated bins. χ²=4.0 over the bins $i=\mathrm{1},\mathrm{2},\mathrm{6},\mathrm{7},\mathrm{13},\mathrm{14},\mathrm{15},\mathrm{17},\mathrm{18}$ in the tail that at most are marginally influenced by N_A(δ_i). The n=9 bins together with initially c=3 fitted parameters compute $f=n-c=\mathrm{6}$ degrees of freedom, which supports the refinement distance being normal distributed $\left({\int }_{\mathrm{4.0}}^{\mathrm{\infty }}{\mathit{\chi }}^{\mathrm{2}}\left(\mathrm{8}\right)=\mathrm{0.68}\right)$.

Appendix E: Statistical analysis of the (depth, depth) match with cubic spline interpolation

Here, we quantify the difference between using linear interpolation and interpolation by cubic splines, which is a widely used scheme. The latter has the benefit of using smooth curves, such that the (depth, depth) curve and its derivatives are continuous, but as discussed above, several factors may cause the real (depth, depth) curve to be non-differentiable or even discontinuous, and we therefore maintain our practice of linear interpolation between the depths of the match points.

Analogous to the definition of δ, we define ${\mathrm{\Sigma }}_i=S(D_{\mathrm{1}},\mathrm{\dots }$, , $\mathrm{\dots },D_n)-S(D_{\mathrm{1}},\mathrm{\dots },D_i,\mathrm{\dots },D_n)$, where $S(D_{\mathrm{1}},\mathrm{\dots },D_i,\mathrm{\dots },D_n)$ is a cubic spline calculated for all match points and S(D₁,…, , …,D_n) a cubic spline calculated for all but the ith match point.

Analogous to the above analysis with linear interpolation, we expect that each match point's depth assignment has an uncertainty due to the varying peak form, which is caused by regional deposition differences and short-term accumulation variations. This peak assignment uncertainty determines the distribution of Σ in sections where the EGRIP–NGRIP1/2 depth curve is straight and no systematic glaciological differences occur. However, varying accumulation conditions at the time and position of snow deposition or different ice-flow patterns influencing the ice cores cause differently evolving annual layer thickness ratios $r_i=\frac{D_{i+\mathrm{1}}-D_i}{d_{i+\mathrm{1}}-d_i}$, which leads to curvature of the (depth, depth) curve. This curvature may not be captured if there are no match points at the relevant depth, and different interpolation schemes will make different predictions across intervals without match points. For the distribution function of Σ, we thus expect two contributions: one from the uncertainty of peak assignment and the other one from the refinement.

Figure E1 Histogram of the difference data with cubic spline interpolation ${\mathrm{\Sigma }}_i=S(D_{\mathrm{1}},\mathrm{\dots }$, , $\mathrm{\dots },D_n)-S(D_{\mathrm{1}},\mathrm{\dots },D_i,\mathrm{\dots },D_n)$ for the match points between EGRIP and NGRIP1/2.

For further statistical treatment, we bin the Σ values. As $\mathrm{\Sigma }\in (-\mathrm{0.674604},\mathrm{0.364446}\mathrm{m})$, we counted the occurrence N(i) of values for the i=1…21 intervals $(-(\mathrm{0.725}+i\times \mathrm{0.05})\mathrm{m},(-\mathrm{0.675}+i\times \mathrm{0.05}\left)\mathrm{m}\right)$ in between −0.675 and 0.375 m and display the data in the following histogram (refer to Fig. E1). The standard deviation of the counts is $\mathit{\sigma }\left(N\right)=\sqrt{N}$. Overall, n=17 bins are occupied. The distribution of Σ is presented in Fig. E1. The weight of the χ² is the counting error $\mathit{\sigma }\left(N\right)=\sqrt{N}$.

Like for the (depth, depth) differences for the linear interpolation δ, the distribution of Σ appears having too-wide tails to match a Gaussian distributions. We support this by trying to fit a single Gaussian normal distribution $N_G\left(\mathrm{\Sigma }\right)=G\exp \left(-(\mathrm{\Sigma }-{\mathit{\mu }}_G{)}^{\mathrm{2}}/{\mathit{\sigma }}_G^{\mathrm{2}}/\mathrm{2}\right)/\sqrt{\mathrm{2}\mathit{\pi }}/{\mathit{\sigma }}_G$ to the data. Even when treating the three values below Σ<0.325 as outliers and restricting the fit to the n=14 bins around 0, we minimize χ²=30.8. For the c=3 fitted constants, the degree of freedom is f=11. This indicates that the fitted distribution is not supported by the data $\left({\int }_{\mathrm{0}}^{\mathrm{30.8}}{\mathit{\chi }}^{\mathrm{2}}\left(\mathrm{11}\right)=\mathrm{0.99}\right)$.

The Gaussian normal distribution for the peak assignment in the centre seems – as for the δ – to be overlaid by a second Gaussian normal distribution representing the refinement of the depth scale. Following the just-posed assumption that the refinement distance is statistically normal distributed N_D(Σ) and that the uncertainty of the peak assignment in smooth intervals is normal distributed N_C(Σ), we tried to χ² fit $N\left(\mathrm{\Sigma }\right)=N_C\left(\mathrm{\Sigma }\right)+N_D\left(\mathrm{\Sigma }\right)$ to the entire data set with c=6 fitting parameters. For the n=17 bins of the entire data set Σ, we minimized χ²=48.7, which suggests the model does not describe the data (for details, refer to Fig. E1, where the figures for a similar treatment as for rejecting N_G(Σ) are provided).

By just treating the three values (out of 377 in total) below Σ<0.325 as outliers, which is justified, when we are mainly interested in assessing the central part of the distribution, which refers to the match-point assignment. Thus, repeating the fit for the bins $i=\mathrm{8},\mathrm{\dots },\mathrm{21}$ (n=14) and labelling the fitted function $N\left(\mathrm{\Sigma }\right)=N_E\left(\mathrm{\Sigma }\right)+N_F\left(\mathrm{\Sigma }\right)$ to clearly distinguish the result from the fit to the entire data set before, we are able to minimize χ²=3.09 for the $f=n-c=\mathrm{11}$ degrees of freedom, indicating that the fitted distribution is supported by the data $\left({\int }_{\mathrm{3.1}}^{\mathrm{\infty }}{\mathit{\chi }}^{\mathrm{2}}\left(\mathrm{11}\right)=\mathrm{0.93}\right)$.

The χ² fit computes for the peak assignment distribution N_E(Σ) a scaling factor $E=(\mathrm{12.6}\pm \mathrm{1.6})$ m, the shift from the centre ${\mathit{\mu }}_E=(-\mathrm{0.0003}\pm \mathrm{0.004})$ m and a standard deviation ${\mathit{\sigma }}_E=(\mathrm{0.044}\pm \mathrm{0.004})$ m. Similarly for the refinement distance N_F(Σ): $F=(\mathrm{6.0}\pm \mathrm{1.5})$ m, ${\mathit{\mu }}_F=(\mathrm{0.004}\pm \mathrm{0.013})$ m and ${\mathit{\sigma }}_F=(\mathrm{0.13}\pm \mathrm{0.02})$ m.

For the linear interpolation scheme (refer to Appendix D), the peak assignment contribution (N_A(δ) dominates the central bins. Only 4 % of the counts are attributed to the refinement (N_B(δ)), and we could define the interval where to perform a direct statistical analysis of the δ values for the peak assignment mode from the intersections of N_A(δ) and N_B(δ). For the cubic spline interpolation scheme, 16 % of the counts are contributing to the refinement N_F(Σ) in the central bins. Thus, the roots ${\mathrm{\Sigma }}_{\mathrm{1}}=-\mathrm{0.089}$ m and Σ₂=0.087 m of N_E(Σ_k)=N_F(Σ_k) do not define the complete inner interval where the peak assignment contributes, and it does not as clearly dominate the distribution of Σ. The peak assignment and the refinement modes do not separate as clearly as for the linear interpolation, and a Shapiro–Wilk test for the (Σ₁,Σ₂) interval fails. When extending the interval to $(-\mathrm{0.16},\mathrm{0.16}\mathrm{m})$ – which covers the contributions of the peak assignment mode well – the direct statistical analysis for the 348 points estimates a standard deviation of 0.056 m, skew of 0.06 and kurtosis of 0.21. This suggest a symmetrical distribution that is slightly higher with wider wings, which is consistent with the above observation of the refinement mode being recognizable in interval, and the standard deviation is overestimated. A Shapiro–Wilk test supports normal distribution in the interval $(-\mathrm{0.16},\mathrm{0.16}\mathrm{m})$, as W is 0.995 and the corresponding p value is 0.38. As we already excluded outliers for the analysis, there is no meaning in the statistical analysis of the refinement distance, and we assume Σ as a measure for the systematic deviation when using cubic spline interpolation.

The cubic spline scheme confirms the above-observed maximal error of the match points $\mathrm{\Delta }D=\mathrm{0.043}\mathrm{m}={\mathit{\sigma }}_A\mathit{\approxeq }{\mathit{\sigma }}_E=\mathrm{0.044}$. When restricting $\left|\mathit{\delta }\right|,\left|\mathrm{\Sigma }\right|\le \mathrm{0.375}$ m, then both δ and Σ have three match points exceeding this threshold and both distributions fit a profile with a standard deviation of σ_F=σ_B (without three outliers) =0.13 m for the refinement each match point contributes on average. This means that one of the interpolation schemes is not superior to the other, but comparing them illustrates the uncertainty associated with interpolation in between the match points.

Now, we calculate the systematic deviation between linear and cubic spline interpolation from data sets in 0.01 m resolution. ζ_i denotes the maximal absolute difference in the interval between the ith and (i+1)th match point, which is a direct measure of systematic differences due to the interpolation schemes.

Figure E2 The (EGRIP, NGRIP) match points, the (EGRIP depth, GICC05 time) scale, the (EGRIP depth, GICC05 MCE) curve, $\left|\mathit{\delta }\right|$, $\left|\mathrm{\Sigma }\right|$, $\left|\mathit{\zeta }\right|$, $\frac{\partial t_{\mathrm{GICC}\mathrm{05}}}{\partial D_{\mathrm{EGRIP}}}$, $\mathrm{\Delta }t\left(D\right)=\left|\frac{\partial t}{\partial D}\right|\sqrt{(\mathrm{\Delta }t{)}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}}$.

$\left|\mathit{\zeta }\right|$ and $\left|\mathit{\delta }\right|$ are both less than 0.4 m and exhibit a similar pattern, while $\left|\mathrm{\Sigma }\right|$ has less in common with both δ and ζ. $\left|\mathit{\zeta }\right|$ is a good measure for the interpolation uncertainty along the record as it is the direct comparison of two fundamentally different interpolation approaches (see Fig. E2).

For linear interpolation, the statistical error for the computed depth in between two match points is limited by the maximal error of the match points ΔD=0.043 m and the error of the interpolated depth D is therefore $\sqrt{(\mathrm{\Delta }D{)}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}}$. To propagate the depth error and estimate the additional error of the time match, we start from the highest-resolution published GICC05 dating of NGRIP with 2.5 and 5 cm depth resolution above and below 349.8 m, respectively (Vinther et al., 2006; Rasmussen et al., 2006), and linearly interpolate the EGRIP depth (D) onto the NGRIP depth (d) to get the timescale t(D) for EGRIP. We calculate $\frac{\partial t_{\mathrm{GICC}\mathrm{05}}}{\partial D_{\mathrm{EGRIP}}}$ in the high-resolution data set and sample it at the match points.

The matching error related to the timescale transfer $\mathrm{\Delta }t\left(D\right)=\left|\frac{\partial t}{\partial D}\right|\sqrt{(\mathrm{\Delta }t{)}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}}$ is maximally about 4 years; it exceeds the MCE on two occasions in the uppermost 200 m by 1 year and becomes increasingly smaller compared to the MCE for increasingly deeper parts of the record (see Fig. E2).