Implementation of counted layers for coherent ice core chronology

. A recent coherent chronology has been built for four Antarctic ice cores and the NorthGRIP (NGRIP) Greenland ice core (Antarctic Ice Core Chronology 2012, AICC2012) using a Bayesian approach for ice core dating (Datice). When building the AICC2012 chronology, and in order to prevent any confusion with ofﬁcial ice core chronology, the AICC2012 chronology for NGRIP was forced to ﬁt exactly the GICC05 chronology based on layer counting. However, such a strong tuning did not satisfy the hypothesis of independence of background parameters and observations for the NGRIP core, as required by Datice. We present here the implementation in Datice of a new type of markers that is better suited for constraints deduced from layer counting: the duration constraints. Estimating the global error on chronology due to such markers is not straightforward and implies some assumption on the correlation between individual counting errors for each interval of duration. We validate this new methodological implementation by conducting twin experiments and a posteriori diagnostics on the NGRIP ice core. Several sensitivity tests on marker sampling and correlation between counting errors were performed to provide some guidelines when using such a method for future dating experiments. Finally, using these markers for NGRIP in a ﬁve-core dating exercise with Datice leads to new chronologies that do not differ by more than 410 years from AICC2012 for Antarctic ice cores and 150 years from GICC05 for NGRIP over the last 60 000 years


Introduction
The reference timescale for Greenland ice cores, GICC05, has been obtained by layer counting back to 60 ka (thousands of years before 1950; Vinther et al., 2006;Rasmussen et al., 2006;Andersen et al., 2006;Svensson et al., 2008).This chronology is absolute, with an increasing associated uncertainty with depth, reaching more than 2.6 at 60 ka.Because this chronology is based on layer counting, the duration of events is rather precise, even for old ages, with an average relative counting error of 0 to 20 %.The uncertainty of the GICC05 age scale, however, cumulates the counting error through the maximum counting error (MCE).It maximizes this uncertainty since it assumes that the layer-counting error is fully correlated from one interval to another (Rasmussen et al., 2006).
This chronology has been used as a reference for many records of the North Atlantic region (Austin et al., 2012;Walker et al., 2012;Austin and Hibbert, 2012;Davies et al., 2012;Blockley et al., 2012b).It has also been used as a basis over the last 60 kyr for the recent construction of the coherent Published by Copernicus Publications on behalf of the European Geosciences Union.
Antarctic Ice Core Chronology (AICC2012) gathering one Greenland ice core (NorthGRIP -NGRIP) and four Antarctic ice cores (EPICA Dome C -EDC; EPICA Dronning Maud Land -EDML; Talos Dome ice core -TALDICE; and Vostok; Bazin et al., 2013;Veres et al., 2013).For the construction of AICC2012 with the Bayesian tool Datice (Lemieux-Dudon et al., 2009, 2010), we have imposed a 1σ deviation for NGRIP of 50 years maximum.Even if such a constraint is artificially too strong compared to the true uncertainty of GICC05, it forces a coherency within 5 years between the NGRIP AICC2012 chronology and GICC05.
Still, the strong tie of AICC2012 to GICC05 has led to some technical problems when optimizing the chronology with Datice.Three glaciological parameters are indeed optimized during this process: accumulation rate, ice thinning and lock-in depth (i.e. the depth at which air is trapped when snow is sufficiently compacted).The Bayesian approach requires starting with first-guess (background) scenarios for the three parameters.They are then modified within their imposed variance range so that the final chronology fits the absolute and relative age constraints for each ice core within error bars.
In practice, to force the NGRIP AICC2012 chronology to fit the GICC05 age scale, the modelled thinning function and accumulation rate of the GICC05 chronology (hereafter DJ-GICC05 scenarios; Vinther et al., 2006;Rasmussen et al., 2006;Andersen et al., 2006;Svensson et al., 2008) have been imposed as background scenarios in Datice.The variance associated with these background scenarios was set to be very small in order to prevent any deviation from the GICC05 timescale.
Even if the uncertainty of the GICC05 timescale is well constrained, this is not true for the DJ-GICC05 scenarios of thinning and accumulation.The thinning function is deduced from a simple Dansgaard-Johnsen (DJ) ice flow model (Dansgaard and Johnsen, 1969;Andersen et al., 2006) that has been parameterized to obtain the best match between the modelled and observed depth-age horizons in the ice cores.Then, the thinning function calculated with the DJ model is used together with the observed annual layer thicknesses to produce an accumulation rate history.No uncertainty value is associated with the reconstructions of thinning and accumulation rate in Greenland ice cores, but thinning and accumulation reconstructed from such 1-D ice flow models are only a first approximation (Cutler et al., 1995;Parrenin et al., 2004Parrenin et al., , 2007)).
Recently, studies combining air isotopic measurements (δ 15 N of N 2 ) with firnification models have suggested that, both in NGRIP and NEEM, the accumulation rates reconstructed from the GICC05 or ss09sea chronologies, through layer counting and the DJ flow model, were overestimated for the last glacial period (Huber et al., 2006;Guillevic et al., 2013;Kindler et al., 2014).Indeed, δ 15 N of N 2 of air trapped in an ice core indicates the depth and the amplitude of abrupt temperature changes in the gas phase through thermal fractionation.The depth difference between the same abrupt temperature changes recorded in the ice phase through ice δ 18 O increase/decrease and in the gas phase through a positive/negative δ 15 N peak is called delta-depth ( depth).Moreover, in the absence of any abrupt temperature change and convection at the top of the firn, the δ 15 N gives an indication of the past lock-in depth (LID) through gravitational fractionation.A firnification model including heat diffusion and mainly driven by temperature and accumulation rate can reproduce long-term and abrupt δ 15 N variations with depth for Greenland ice cores.The same would not be true in Antarctica, where a strong discrepancy between firnification models and data is observed (Landais et al., 2006).For NGRIP, it has been shown that the δ 15 N profile is best reproduced when the ss09sea accumulation rate is decreased by ∼ 20 % over the period 20 to 60 ka (Kindler et al., 2014;Huber et al., 2006).
It thus appears that the way NGRIP was implemented in the Datice tool for the AICC2012 chronology is not optimal.In addition to GICC05 chronological uncertainties that were not taken into account by construction, imposing the DJ-GICC05 accumulation rate and thinning scenarios with artificially reduced variances most probably led to incorrect output scenarios for these glaciological parameters.
In this paper, we propose an improvement of Datice to better implement the chronological uncertainties.Markers of duration are integrated in Datice with associated counting error.This allows for the strong constraints on thinning and accumulation rate to be relaxed.It also allows the NGRIP chronology to differ from GICC05 chronology within its error bars.
The outline of the manuscript is as follows.First, a methodological section presents and validates the improvements made to the Datice tool in order to integrate the duration constraints with their uncertainties.Then, we discuss different ways to implement the counting errors within the global chronological uncertainty.We also present some sensitivity experiments using the modified Datice tool for optimizing the sampling strategy and correlation between counting errors.Finally, we focus on how this new version of Datice modifies the NGRIP and the four Antarctic ice core chronologies compared to AICC2012.

Methodology
The purpose of the following section is to describe the modifications implemented in Datice (Lemieux-Dudon et al., 2009, 2010) to take into account the duration constraints.This type of marker enables one to constrain the duration of depth intervals along ice cores.This constraint is applied by feeding Datice with the beginning and end depths of the  Ide et al., 1997).Datice requires that the background age model and the palaeoobservations be independent from each other, since the cost function J is derived from the Bayes theorem: where P (Y |X) and P b (X) are the likelihood and prior probability distribution (Tarantola, 2005).
In practice, Datice is applied to several ice cores with a large set of palaeo-observations in order to calculate coherent chronologies for both the ice and gas phases.The chronologies are deduced from the scenarios of three glaciological parameters for each core indexed with k (Appendix Sect.A): (i) T k the total thinning function, (ii) A k the accumulation rate, and (iii) C k the lock-in depth in ice equivalent (LIDIE).To run a Datice experiment, palaeo-observations and the background parameters T b k (z), A b k (z) and C b k (z) must be provided with their respective uncertainties.The minimization of the cost function J enables one to refine the background by identifying correction functions τ k (z), α k (z) and γ k (z) at each depth level z k : From a particular set of correction functions, one can deduce a particular age model.Hereafter, a particular age model is written The Datice cost function formulation (Eq.5) relies on the following important statistical assumptions.In the prior probability distribution of Eq. (1), the parameters T k , A k and C k are supposed to be independent and log-normally distributed, with medians equal to , respectively.The prior probability distribution is further rewritten in terms of the correction functions (Eqs.2, 3, 4), to which we apply the change of variable X = log(X) in order to transform lognormal into normal probability distributions (pdf's;Tarantola, 2005).Using this change of variable, since observations of different types are supposed to be independent with either normal or log-normally distributed errors, the likelihood of Eq. ( 1) is itself a product of normal pdf's.Under these assumptions, the cost function J sums up quadratic terms (Eq.5).
Until now, observation Y could be of the following types: ice and gas age markers ( ia and ga ), delta-depth markers ( dd ), or ice and gas stratigraphic links ( is and gs ; Lemieux-Dudon et al., 2010;Buiron et al., 2011;Veres et al., 2013;Bazin et al., 2013).The application of the duration constraints ( ad ) leads to an additional term in the cost function (Eq.5), with special care to preserve the Datice hypothesis of no error correlation between (i) observations of different types, (ii) observations from different cores, or (iii) observations and background model scenarios: In Eq. ( 5), the first term measures the distance between the current age model X and the background scenarios X In this article, we wish to design Datice experiments with duration constraints derived from the GICC05 counted layer chronology.In Sect.2.3, we especially investigate the setting of the observation error covariance matrix R ad associated with these markers: where R ad ij accounts for the error covariance between the i th and j th pair of markers Y ad i and Y ad j .σ ad i and σ ad j are their standard deviations and ρ ad ij is their error correlation coefficient.

Validation of Datice developments: twin experiments
In this section, twin experiments are performed to test the incorporation of duration constraints within the Datice tool.
A twin experiment enables one to test any data assimilation system.It consists in the construction of some synthetic data and background by applying random perturbations of known statistical distribution to a given model scenario.The unperturbed model scenario is referred to as the "truth".The aim of this validation method is to rebuild the truth by running the data assimilation system on the perturbed data and background.
In our case, we designed a twin experiment based on 51 simulations where the Datice system is run with the NGRIP ice core alone.The duration constraints are the only type of observation included.The GICC05 age scale is considered as the "truth".We construct synthetic observations and backgrounds by applying random perturbations to the "true" GICC05 model scenario.The objective is to run the Datice system with the synthetic data as input and to rebuild GICC05 as accurately as possible.
Each twin experiment inputs are prepared using the following method.To build the 51 sets of synthetic markers of duration constraint Y ad , we first sample the "true" duration constraints Y ad, t from the GICC05 age scale every 100 years.
For this experiment, the markers Y ad, t represent the "truth" (superscript t) as extracted from the "true" age model GICC05.We then construct the observation error covariance matrix R ad (see Sect. 2.1) based on the MCE data under the assumption of full error correlation (see details in Sect.2.3).To provide the markers of duration Y ad that will effectively be applied in the simulations, the "true" markers Y ad, t are perturbed within their uncertainty range through random normal perturbations specified according to the observation error covariance matrix R ad : In the same way, the 51 background scenarios are built from the "true" GICC05 thinning function and accumulation rate by applying random perturbations, with the particularity that Background and analysed chronologies thinning function and accumulation rate are log-normally distributed as discussed in Sect.2.1.We construct the random lognormal perturbations δ α and δ τ on the basis of the background error covariance matrices B τ and B α , which are specified according to Bazin et al. (2013) with adapted values (see Table 4 and Appendix Sect.B): where B α and B τ are the first two diagonal and uncorrelated blocks of matrix B introduced in Eq. ( 5) of Sect.2.1.
To construct the perturbed background thinning function T b and accumulation rate A b , the δ α and δ τ vectors are applied as multiplicative factors to the "true" GICC05 thinning function T t i and accumulation rate A t i at each depth level z i (with index i running from 1 at the top of the core to n at the bottom): Figure 1 shows the large spread of the resulting perturbed background age scales (dashed lines) and a superimposition of the corresponding analysed age scales (orange lines).
Figure 2 shows the difference between the set of analysed chronologies minus GICC05 (upper panel) and the error of the analysed chronologies, i.e. σ a the a posteriori standard deviation calculated by Datice (lower panel).Histograms of the background and analysed chronologies are shown in Fig. 3 for the 1800 m depth level.Figures 1, 2 and 3 all show the convergence toward GICC05 even though Datice is fed with perturbed background scenarios and duration constraints.In the Datice system, the calculation of the analysed error σ a relies upon Clim. Past, 11, 959-978, 2015 www.clim-past.net/11/959/2015/ the assumption of normally distributed errors.This may be a strong assumption.However, the histogram of analysed chronologies (i.e.output chronologies) is rather symmetric and centred on GICC05 compared to the very asymmetric histogram of the perturbed background age scales (Fig. 3).
Moreover, at depth level 1800 m, 96 % of the 51 output chronologies are located inside a ±2 σ a envelope centred on GICC05 (Fig. 3).This result gives confidence in the methodology applied to calculate the analysed error (see Appendix Sect.E).A larger number of samples would be helpful to refine this analysis.
Finer diagnostics confirm the reliability of the Datice methodological developments.As investigated in Desroziers et al. (2009), several levels of a posteriori diagnostics can be applied on data assimilation system on the basis of ensemble of analyses conducted on ensemble of perturbed background and observations.The construction of our twin experiment appropriately relies on an ensemble of both perturbed background and observations.It consequently enables one to verify the first level of these diagnostics (Desroziers et al., 2009).It states that for weakly non-linear observation operators h (see Eq. ( 5), Sect.2.1), when both the B and R matrices are calibrated, averaging the values of the cost function at the optimum (when X a optimum is reached) must be equal to the number of observations p: In our twin experiment, we apply 633 duration constraints.The average of the cost function at optimum X a over our 51 simulations gives a value of 626 in accordance with Eq. ( 12).This is quite a fair result, and it validates our methodological development.One should note that we have applied perfectly calibrated background and observation error covariance matrices.Indeed, the background and observation errors specified in the cost function are exactly the B and R matrices that have been used to produce the synthetic backgrounds and observations based on the true scenario.In a more complex experiment, the B and R matrices are usually misspecified because the background and observation errors, b and o , are usually poorly known since the truth itself, X t , is the unknown (see Appendix Sect.D1).In such cases, the a posteriori diagnostics are applied to calibrate the error covariance matrices.In future work we wish to conduct such calibration on Datice experiments involving several ice cores.

Implementing layer-counting error (MCE)
Layer counting consists in identifying annual cycles on the basis of annual layer proxies recorded along the core.The identification of annual cycles is subjected to errors.In order to deal with uncertain annual layers and to derive a counting error estimate, GICC05 adopted the following statistical approach.If the i th cycle is identified as a certain annual cycle, the layer is counted as a full year with a zero error.Otherwise, for an uncertain i th annual cycle, the layer counts as half a year plus or minus half a year.For each cycle numbered with index i, the two following variables n i and σ i are introduced in order to record the layer detection and associ- n i ± σ i = 0.5 ± 0.5 yr for an uncertain layer.( 14) Ages along the ice core can thus be inferred by summing up the n i cycles.For instance between depths z q and z p , delimiting the start and end depths of the q th and p th individual cycles, respectively, the duration Y z q ,z p (in years) is written as The official GICC05 age scale provides depth levels and cumulative counting error corresponding to time windows of 20 years.The GICC05 error estimate is called the maximum counting error (MCE).The MCE sums up the error of the individual cycles (i.e.σ i ) over the corresponding time window: For our experiments, the objective is to apply the GICC05 measure of duration Y 20 yr z q ,z p as duration constraints in Datice simulations.Two questions arise at this stage: -Over which time window should we define our GICC05 markers of duration?Shall we apply markers of duration on 20 yr time window or choose another sampling rate (i.e.20-, 40-, 60-year time window)?
-How should we infer the associated error when applying different time windows?
Neither of these questions are trivial.They are closely interlinked through the existence of error correlation between annual layers and the assumptions inherent to the MCE construction.
The MCE can be expressed through formulation of the GICC05 counting process with two normal probability density functions (pdf's): (i) the pdf of annual cycles identified as certain with a 1-year mean and a variance that tends to zero and (ii) the pdf of annual cycles identified as uncertain with a mean and standard deviation both set to half a year.Under this formalism (Appendix Sect.C1), the calculation of the error z q ,z p on any counting measure Y z q ,z p is well documented, and the role played by the error correlation between annual cycles n i and n j becomes quite clear.If ρ ij records such correlation, the z q ,z p error is written as The z q ,z p error reaches a minimum value in the case of a null error correlation between any pair of cycles (i.e.ρ ij = 0): Conversely, the error reaches a maximum value when the error correlation between annual cycles is maximum (i.e. The MCE calculation is based upon Eq. () and is therefore an upper estimate of the error regarding the value of the correlation coefficient (but not regarding the assumptions on the error σ i ).In particular, the error correlations have an "infinite range" along the core (ρ ij does not decrease with the distance between the measured cycles) and pairs of annual cycles are fully correlated regardless of their respective position.Such description is not entirely realistic.Still, Rasmussen et al. ( 2006) have acknowledged that the assumption of full correlation of counting errors is not correct and stated that "recognizing that the counting errors in reality are neither uncorrelated nor fully correlated, we adopt the simple and conservative approach, summing up the uncertainties as if they were correlated" .In this study, the 1σ uncertainty of the GICC05 ice core is considered as half the MCE.Following this approach means that errors for duration constraints at 40, 60 or 80 years will be derived by summing up the GICC05 20-year-window MCE 2, 3 and 4 times, respectively, in the case of full correlation within the time windows associated with the chosen sampling rate (Appendix Sects.C3 and C2).
However, the final chronology error should not depend on the arbitrary choice of the sampling rate.The option has thus been included in the Datice approach to apply error correlation on a finite interval and avoid abrupt cut-off of error correlation between adjacent intervals.This development should permit sampling of the markers with a certain step and apply error correlations beyond this time interval.Indeed, in future chronologies constructions, the value of the error correlation may change along the core in relationship with changes of climatic periods.
For this formulation of error correlation on a finite range, the correlation coefficients ρ ad ij of the observation error covariance matrix R ad (Eq.6) are set according to a correlation function f that smoothly decreases with the distance between two duration constraints Y ad i and Y ad j : Clim.Past, 11, 959-978, 2015 www.clim-past.net/11/959/2015/ The shape of the function f is chosen as the product of a Gaussian and a triangular function: where L ad must be set in metres in order to adjust the width of the f function and therefore the scope of the error correlation: the larger L ad is, the higher the correlation between markers of duration.
With this new formulation of the error correlation, we can explore how both sampling and error correlation independently affect the final chronology and provide some guidelines for future constructions of ice core chronologies.

Tests and optimization of the Datice system to apply the GICC05 duration constraints
In this section, we extract several sets of duration constraints from GICC05, with different sampling and/or different assumptions regarding their associated errors.These inputs are used to conduct multiple Datice experiments and thus to investigate the sensitivity of the solution to sampling and error correlation assumptions.In the two following sections, experiments are run on the NGRIP core alone with only duration constraints.Details on the background settings are provided in Table 4.The marker errors are derived under either (i) the full correlation assumption AddMCE or (ii) the assumption of non-correlation beyond the 20-year-window SqrAddMCE (see Sect. 2.3 and Appendix Sect.C3).
In these experiments, a classical 1 m depth grid resolution is imposed, as in AICC2012.On such a depth grid, the annual layer thickness drops below 0.05 m yr −1 at some depth level so that the number of years in a 1 m layer becomes larger than 20 years.Datice cannot handle markers of duration that are sampled below the depth grid resolution.This technical issue prevents us from applying the GICC05 20-year-window markers and MCEs directly (Eqs. 16 and 17).In order to test fine sampling of the duration markers whilst avoiding sampling resolution below the depth grid, we implemented an adaptive sampling ranging from 40 to 140 years back to 60 ka.

Sampling and error correlation influence
To study the influence of the sampling, we run three experiments with markers sampled at three uniform rates (100, 200 and 300 years) as well as one experiment with the adaptive sampling between 40 and 140 years.
In these four experiments, the associated errors are derived from the 20-year-window MCE data under the AddMCE assumption of full error correlation between annual cycles over the length of the sampling interval.As discussed in Sect.2.3, error correlation and sampling are interlinked together.To investigate the error correlation influence, we run a second This assumption corresponds to the addition of squared MCE over 20-year windows to obtain the squared error over the sampling interval (40 to 140 years) used as input for Datice.Table 1 summarizes the experiment configurations.
Figures 4 and 5 show the different NGRIP simulations.As expected and discussed in Sect.2.3, the age solutions and their associated errors are sensitive to the sampling.For the four experiments run under the AddMCE assumption, we better reproduce the GICC05 details with finer sampling rates (Fig. 5).Still, finer sampling of duration constraints is not the main reason for the better agreement with GICC05.Indeed, as error correlations are cumulative under the Ad-dMCE assumption, the observation error largely increases with the length of the marker sampling window.Consequently, the strength of the constraint decreases, which deteriorates the convergence toward GICC05.The impact of observation error is also illustrated by Fig. 4 for the case of comparing the two adaptive sampling simulations, i.e. 40-140 yr_AddMCE and 40-140 yr_SqrAddMCE.The simulation run under the assumption of non-correlation beyond 20 years (SqrAddMCE) better converges toward GICC05, with a smaller associated error.As expected, the SqrAddMCE assumption strongly reduces the observation error at any depth along the core with respect to the full correlation assumption AddMCE.
Option SqrAddMCE may therefore be a way to relax the dependence of the analysed error to the sampling.However, as mentioned above, the abrupt loss of correlation at the boundaries of sampling interval may be questioned.At the junction of two duration constraints, neighbouring annual cycles from either side do not share any error correlation, while each of them correlates with much distant layers (as long as these layers are included in the same sampling interval).We actually rather expect error correlations to smoothly decrease with the distance between annual cycles.To circumvent this problem, we have designed an experiment called CorrCoeff_40-140 yr_AddMCE (shown in Figs. 4 and 5), which implements a correlation coefficient smoothly decreasing with the distance between markers.This implementation is discussed in the next section.We have demonstrated the sensitivity of the solution to the sampling and to the MCE assumptions applied to derive the observation error.However both issues are not fully decoupled in this first illustration.Hereafter, we investigate possible ways to study the error correlation independently of the sampling.

Finite range versus infinite range error correlation influence
In this section, we apply different correlation coefficients between duration constraints as implemented in Datice (Eqs.6 and 21).In the following experiments, we investigate two correlation configurations: (i) correlation coefficients with infinite depth range along the core (hereafter InfiniteRangeCorr) and (ii) correlation coefficient smoothly decreasing with the distance between markers (hereafter FiniteRangeCorr) (Eq. ).
In this set of experiments, due to large correlation coefficient values, the level of observation error is largely increased compared to the experiments of Sect.2.4.1.To operate with configurations where the minimization and solution are still strongly driven by the constraint of the markers of duration, the background errors have been exaggerated (Table 4).With such a configuration, the analysed error should tend toward the observation error (Appendix Sect.E1): where σ b , σ o and σ a are the background, observation and analysed errors, respectively.In a first set of experiments, we investigate the InfiniteRangeCorr option.The correlation coefficient ρ ad ij (Eq.21) is set to constant values ranging from 0.2 to 0.8.Such a configuration implies identical error correlations between markers separated by a large or a small distance, as the MCE formulation does for the GICC05 chronology.The MCE formulation implies a correlation coefficient of 1 all along the ice core.In the Datice approach, it is technically impossible to attribute a value of 1 to ρ ad ij due to the R ad matrix inversion in the cost function formulation (Eq.5).As a consequence, the Datice experiment run with a correlation coefficient of 0.8 is the closest analogue to the MCE formulation, and we expect the analysed error to closely approach the MCE.
Figure 6 shows the InfiniteRangeCorr experiments with (i) a comparison of the background and analysed chronologies with the reference chronology GICC05 as well as (ii) a comparison of the Datice analysed errors and the MCE.As expected, the analysed errors tend toward the MCE for higher correlation coefficients.The full convergence to the MCE values is, however, hampered since more error correlation between markers progressively rules out the hypothesis of Eq. ( 23): when the observation error becomes too large, the analysed error is also driven by the background error (Eq.E11).Analysed chronologies also show some predictable behaviour.When the correlation coefficient increases, the confidence in the duration constraints decreases and the analysed chronologies stay close to the background chronologies.Importantly, the reconstructed chronologies show an increasing bias relatively to GICC05 with increasing correlation coefficient.This problematic bias is closely related to the infinite depth range of the correlation coefficient.
We test hereafter the FiniteRangeCorr experiment with the finite depth range correlation coefficient.We ran simulations with five different types of sampling: (i) four uniform sampling rates (300, 200, 100 and 80 years) and (ii) the adaptive sampling between 40 and 140 years (CorrCoeff_40-140 yr_AddMCE experiment).An error correlation is applied between markers according to Eq. ( 21), and for all experiments we have set the correlation length L ad to 300 years (Eq. 22). Figure 7 shows a comparison between background, analysed chronologies and GICC05 as well as analysed errors.As expected and already shown in previous sections, we clearly observe that there is a good fit to GICC05 with a better resemblance for the highest sampling rate.Importantly, in comparison to the InfiniteRangeCorr experiments, the analysed chronologies are not biased relative to GICC05.Moreover, despite the different sampling rates, the analysed errors show very similar values contrary to tests presented in the previous section.This is due to the fact that the analysed error is mainly influenced by the correlation coefficient on a finite length, which more efficiently intercorrelates markers sampled on short time windows.
In summary, the tests presented in the two previous sections suggest some guidelines for future constructions of chronology using the duration constraints.The central problem is the definition of the error associated with annual layer counting and how this error is correlated with other layers' error.We showed that making different assumptions on the error correlation leads to significant difference in the final chronology and associated error.For experiments with Datice applied to several ice cores including NGRIP, if the objective is to preserve the NGRIP age scale, our recommendations are (i) to sample the duration constraints over small time windows (e.g. 100 years or apply an adaptive sampling rate), (ii) to use a small uncertainty for the observations (this is directly linked with a large or short range of correlation be- tween layer-counting errors), or (iii) to increase the NGRIP overall background error.

Application to five site experiments and comparison with AICC2012
After having validated the new developments for the implementation of duration constraints and possible error correlation, we show a first application of the new Datice tool to a five-ice-core experiment (NGRIP, EDC, EDML, Vostok, TALDICE).
To use Datice properly, the age constraints and the background scenarios need to be independent from each other.This was not the case when building AICC2012 for the NGRIP ice core.Here, the new development of Datice allows one to use scenarios for background accumulation rate and thinning function independently from the age constraints deduced from GICC05 for NGRIP.In this application, the thinning function is the same as for AICC2012, obtained from the 1-D DJ glaciological model adapted to NGRIP (Andersen et al., 2006).However, we have largely increased its associated variance to make it comparable with variances associated with the background thinning function of the other cores implemented in Datice.For the accumulation rate, we use the ss09sea accumulation rate based on the water isotopes record (Johnsen et al., 2001).The variance of NGRIP background accumulation rate is similar to the variances of background accumulation rates of the other ice cores in the building of AICC2012.The LIDIE background scenario in AICC2012 was built from a firnification model (Goujon et al., 2003) whose input parameters (temperature and accumulation rate) were roughly adjusted to be coherent with the mean δ 15 N values measured over the NGRIP ice core.It is thus independent of GICC05 and has been kept unchanged for our study.
Concerning the age constraint, the absolute age markers deduced from GICC05 were replaced by the duration constraints.The markers of duration are obtained from the GICC05 chronology with adaptive lengths of intervals between 40 and 140 years under the AddMCE assumption (full correlation between annual cycles) and a correlation length of 300 years.In order to constrain the relative gas chronology vs. the ice chronology, we use information derived from δ 15 N of air trapped in ice bubbles.New δ 15 N data on the NGRIP ice core have been published since the AICC2012 chronology (Kindler et al., 2014).In particular, these data allow for identification of depths of rapid temperature increases associated with the beginning of Greenland Interstadial (GI) 1 to 7 in the gas phase.The depth differences between peaks of δ 18 O ice and δ 15 N of a concomitant event recorded in the ice and the gas phases are thus used as delta-depth ( depth) constraints.With the new set of data from Kindler et al. (2014), we were thus able to deduce new depth markers that were not available for the construction of AICC2012 (Table 3).Their uncertainties depend on the resolution of measure- ments and the difference of depth estimates.Indeed, the depth can be estimated from the difference between midslopes of δ 18 O ice and δ 15 N increases or from the difference between the maxima of δ 15 N and δ 18 O ice .
Figure 8 compares the new Datice chronology produced here (NGRIP-free) to AICC2012 for the five sites between 35 and 48 ka.We emphasize that the NGRIP-free chronology discussed here should not be taken as a new official chronology.It is only a test for our methodological development.Moreover, the AICC2012 chronology has the strong advantage of being in exact agreement with the GICC05 chronology and hence facilitating the multi-archive comparison, taking GICC05 as a reference, as has already been done in many studies (INTIMATE project: Blockley et al., 2012a).When looking at the NGRIP ice records, the final NGRIP-free chronology does not differ from the GICC05 or AICC2012 chronologies by more than 150 years over the last 60 kyr (Fig. 8).
The Antarctic chronologies are not much modified compared to the AICC2012 chronologies.They all differ by less than 410 years from AICC2012 (Fig. 8), which is well within the uncertainties of these chronologies (400-1000 years over this period).The small differences between the NGRIP-free and AICC2012 chronologies mean that the relationship between Greenland and Antarctic climate discussed with AICC2012 for the millennial-scale variability of the last glacial period stays valid on NGRIP-free (Veres et al., 2013).We observe a classical see-saw pattern with Antarctic temperature increasing during the Greenland stadials, with a faster and shorter increase at EDML than at EDC (Fig. 8).

Conclusions
The Bayesian tool Datice used for the construction of coherent ice core chronology has been improved and now enables one to consider the duration of events as dating constraints.We validated this new methodological implementa-    (Stenni et al., 2011), EDML δ 18 O (EPICA Community Members, 2006Members, , 2010)), Vostok δD (Petit et al., 1999) and EDC δD (Jouzel et al., 2007) water isotopes on different coherent chronologies (AICC2012 in dark blue and NGRIPfree in light blue).The differences between the NGRIP-free and AICC2012 chronologies for each sites are represented by the black lines.
tion by conducting twin experiments and a posteriori diagnostics.
In comparison to age markers, duration constraints are more coherent with the building of chronologies based on layer counting where the absolute error, defined as the maximum counting error for GICC05, increases with depth due to cumulative effects.To account for the fact that the counting errors on duration constraints are neither fully correlated nor uncorrelated, we have also introduced the possibility to adjust correlation between duration errors with a correlation co-efficient that smoothly decreases with the distance between markers.
There is no objective way to choose the best representation of the correlation, and future dating experiments may propose different correlation coefficients for layer counting performed at different periods (glacial vs. interglacial times).We have thus presented here some sensitivity tests for the sampling and correlation of errors associated with duration constraints.These tests lead to general guidelines for future dating experiments including layer counting as absolute age constraints.For example, to best respect an ice core chronology based on layer counting, we would favour a highfrequency sampling of duration constraints with a correlation on a finite depth range.Finally, the comparison of AICC2012 with the chronology obtained over five polar sites using the improved Datice tool incorporating duration constraints and associated correlation of errors shows differences of less than 410 years over the last 60 kyr, well within the uncertainties associated with the AICC2012 chronology.Huge efforts in annual layer counting were produced in the recent years for ice core chronologies, in particular for the Western Antarctic Ice Sheet (WAIS) ice core (WAIS Divide Project Members, 2013) With the objective of better handling the MCE data, we make Gaussian assumptions and reformulate the GICC05 layer counting with two probability density functions (pdf's): -The duration of an annual cycle identified as certain is normally distributed with a 1-year average and a zero standard deviation.
-The duration of an annual cycle identified as uncertain is normally distributed with a mean and a standard deviation both set to half a year.
The counting variables n i and σ i are statistical parameters of the Gaussian distribution, i.e. the mean and standard deviation.
It should be noted that such a formulation may be questioned: (i) the annual layer counting has a discrete underlying nature, and one might rather prefer to introduce discrete random variables to handle it; (ii) the Gaussian pdf applies to continuous random variables ranging from −∞ to +∞, which is far from being the case here; and (iii) Gaussianassumption-based theorems are tricky to apply in the zerovariance limit assessed for annual layer identified as certain.
C2 Sampling markers of duration and amount of error correlation accounted for: general case Sampling markers of duration from the GICC05 layercounted chronology may lead to a different amount of error correlation between individual measures of annual cycles.Let us first sample the markers on a T -time window and get the two constraints, Y T z q ,z p and Y T z p ,z m , which measure the annual cycles in years in the two neighbouring depth intervals z q , z p and z p , z m along the core.The errors T z q ,z p and T z p ,z m associated with each marker are written as If we double the sampling rate (2T -time window markers), we get a single marker Y 2T z q ,z m (instead of the two constraints Y T z q ,z p and Y T z p ,z m ), which measures the annual cycles in years over the depth interval z q , z m .
The error 2T z q ,z m associated with 2T -time window marker Y 2T z q ,z m is now written as Rearranging Eq. (C4) in terms of the errors T z q ,z p and T z p ,z m In Eq. (C5), the last term corresponds to a part of the error accounted for in the 2T -window marker Y 2T z q ,z m that will never be accounted for in the case of the T -window markers Y T z q ,z p and Y T z p ,z m .It corresponds to error correlations between annual layers i and j that are separated by the longest distance as they are located in the [z q , z p ] depth interval for the first layer, and in the next interval [z p , z m ] for the second.The longer range correlation can only be accounted for with the larger sampling rate.This point is illustrated with Fig. C1.
It is worth noting that Eq. (C5) further simplifies with respect to the ρ ij correlation coefficients of the last term: -When the correlation coefficients are identically null, we get the sum of the squared errors: -When the correlation coefficients are set to 1, we get the squared sum of the errors: Yq,m ± Σq,m

Figure C1.
Error covariance matrix R associated with a duration constraint Y 2T z q ,z m sampled at the 2T years rate on a layer-counted chronology, e.g.GICC05.The matrix stores error information related to the measures of annual cycles on the depth interval [z q , z m ].The diagonal elements record the error variances σ i 2 associated with each identified annual cycle, while the non-diagonal elements store the error covariances, with especially the error correlation coefficient ρ ij between pairs of annual layers i and j .The error 2T z q ,z m associated with marker Y 2T z q ,z m takes into account the whole error correlations stored in the R matrix.If the measures of duration are instead sampled at the T sampling rate (i.e.half the previous rate), the marker of duration Y 2T z q ,z m splits into two markers: (i) Y T z q ,z p (in blue) and (ii) Y T z p ,z m (in brown).The error T z q ,z p associated with Y T z q ,z p will only account for the correlation of the upper diagonal block of R (dashed blue line around block).Symmetrically, the error T z p ,z m associated with Y T z p ,z m will only account for the correlation of the lower diagonal blocks of R (dashed brown line around block).Correlations of the non-diagonal blocks of R, which correlate annual layers i ∈ [z q , z p ] and j ∈ [z p , z m ], are only accounted for in the total error when applying the 2T sampling rate.
-Option one: we believe that the full error correlation assessed over the 20-year time window between annual layers cuts-off.Then, no correlation exists between the annual cycles included in the two separated but adjacent depth intervals z q , z p and z p , z m .Under this assumption, the theory shows that we must sum up the squared 20-year MCEs: -Option two: we believe that the full error correlation assessed over the 20-year time window between annual layers extends over the 40-year time window (which means over the depth interval z q , z m = z q , z p ∪ z p , z m ).In that case the theory shows that we must sum up the 20-year MCEs: From this simple illustration, it follows that markers of duration and errors sampled on GICC05 at different rates (i.e.40-60-80-100 years), derived by summing up the GICC05 20-year-window MCE, must be understood as very different inputs and different simulation outputs must be expected.
In Sect.2.3 of the main text, we refer to the AddMCE assumption (Eq.C9) when MCEs are added, while at the same time we refer to the AddSqrMCE assumption (Eq.C8) when the squared MCEs are added.To approximate the a posteriori error of the analysed chronology, the covariances of errors of X a are required.These covariances of errors are recorded in P a , the analysed error covariance matrix, which can be approximated: where B and R are the background and observation error covariance matrices, respectively Eq. (D1), and where H is the tangent linear observation operator (linearization of h at X a ).
Datice calculates the components of P a at each depth level on the basis of Eq. (E4).Importantly, P a operates a balance between the background and observation errors.The P e error covariance propagates to the analysed chronology a X a .
If we define E a as the a posteriori error of the analysed chronology, the corresponding analysed error covariance matrix a is by definition We intend to show how matrix a depends on matrix P a , and then on the error matrices B and R. We reiterate the steps to show this link, as described in Lemieux-Dudon et al. (2009).
One can first linearize the age scale of Eq. (A1) around X a Eq.D5): Inserting Eq. (E7) in Eq. (E5) enables one to approximate the a posteriori error E a : Importantly, the later approximation is valid if a represents sufficiently small perturbations, i.e. the correction functions X a must be close to the true scenario X t .Under this strong assumption, Eq. (E8) leads to Finally, from Eq. (E9), one can approximate the error matrix of the ice age a by applying the expected value operator to Eq. (E9) and by using Eq.(D6): a ∼ X X a T P a X X a .(E10) Datice applies Eq. (E10) to approximate the covariances of errors of the analysed chronology.This approximation especially requires that the optimum correction functions X a obtained after the minimization of the cost function remains sufficiently close to the true scenario X t .On the assumption of normally distributed errors, matrix a provides the standard deviation of the analysed age scale.The process to calculate the analysed error of the gas age scale is similar but relies on Eq. (A3).

E1 Balance between background and observation error and impact on the analysis
The variances of errors of the analysed chronology cumulate the error covariances recorded in matrix P a (Eq.E10).
The age solution and its error are therefore largely determined by the balance between observation and background errors (Eq.E4).To fix ideas, instead of matrices P a , R and B, let us suppose that we deal with the scalars σ a , σ o and σ b .With such a simplification, Eq. (E4) is written as According to the ratio between observation and background errors, there are two extreme configurations: -If σ o σ b , the minimization and the solution are strongly constrained by the observation, and the analysed error tends to be the observation error:

Figure 2 .
Figure 2. Twin experiments: 51 analysed chronologies, i.e. output chronologies from Datice.Top: comparison of the 51 analysed chronologies with GICC05.Bottom: the corresponding 51 analysed errors (red).The dashed black line represents the maximum counting error associated with GICC05 and considered as equivalent to a 2σ uncertainty.

-
Conversely, if σ b σ o , the background scenario dominates and the solution is close to the background.The analysed error tends to the background error: σ a ∼ σ b .(E13) www.clim-past.net/11/959/2015/Clim.Past, 11, 959-978, 2015 An intermediate ratio of background to observation error leads to an intermediate analysed solution and error.In the special case of an equal number of errors in the observation and background, i.e. σ = σ o ∼ σ b , the analysed error is writ--past.net/11/959/2015/ interval, its duration, the duration uncertainty, and optionally the error correlation between markers.Datice aims at obtaining the best age model scenario by formulating an optimization problem with a cost function that is accounted for by two main types of constraint: the palaeo-observations Y and a first-guess age model X b (referred to as the prior or background; Clim.Past, 11, 959-978, 2015www.clim-past.net/11/959/2015/

Table 1 .
Summary of the simulation configurations for the experiments of Sect.2.4.1.Sampling and error correlation influence.

Table 2 .
Summary of the simulation configurations for the experiments of Sect.2.4.2.Finite range versus infinite range error correlation influence.
and 200-year rates are shown with the blue to pink lines.The adaptive sampling between 40 and 140 years is shown in brown.The MCE assumption is AddMCE (full error correlation between cycles).Table 2 summarizes experiment configurations.

Table 4 .
Summary of the simulation configurations.
. Future dating experiments should thus benefit from the methodological development and validation of the Bayesian tool presented in this study.