IceChrono1: a probabilistic model to compute a common and optimal chronology for several ice cores

. Polar ice cores provide exceptional archives of past environmental conditions. The dating of ice cores and the estimation of the age-scale uncertainty are essential to interpret the climate and environmental records that they contain. It is, however, a complex problem which involves different methods. Here, we present IceChrono1, a new probabilistic model integrating various sources of chronological information to produce a common and optimized chronology for several ice cores, as well as its uncertainty. IceChrono1 is based on the inversion of three quantities: the surface accumulation rate, the lock-in depth (LID) of air bubbles and the thinning function. The chronological information integrated into the model are

Published by Copernicus Publications on behalf of the European Geosciences Union.
However, prior to the interpretation of polar ice core records is the complex task of building two robust timedepth relationships, one for the tracers measured in the ice phase (e.g., water isotopes, particulates and chemical impurities) and one for those measured in the air phase (e.g., greenhouse gas concentration, isotopic composition of gases).The firn, where snow is gradually compacted into ice, constitutes the upper 50-120 m part of ice sheets.The firn is permeable and air is only locked in at its base, at a depth level called the lock-in depth (LID).As a result, the entrapped air is always younger than the surrounding ice at any depth level.Through gravitational fractionation processes, LID is closely related to the isotopic composition of δ 15 N of N 2 in air bubbles data (e.g., Buizert et al., 2015;Goujon et al., 2003;Parrenin et al., 2012a;Schwander et al., 1993).The temporal evolution of the age difference between ice and air at a given depth must therefore be estimated using firn densification modeling and air δ 15 N.This age difference is essential for clarifying the exact timing between changes in atmospheric CO 2 concentration and Antarctic surface temperature during deglaciations ( e.g., Caillon et al., 2003;Landais et al., 2013;Monnin et al., 2001;Parrenin et al., 2013;Pedro et al., 2011).However, glacial-interglacial Antarctic firn changes remain poorly understood (e.g., Capron et al., 2013).
Several strategies have been developed to build ice and gas chronologies.We briefly describe these methods, their strengths and caveats hereafter.
-Annual layer counting (e.g., Rasmussen et al., 2006;WAIS Divide Project Members, 2013;Winstrup et al., 2012).Only applicable when accumulation rates are sufficiently high to make this annual layer identification possible, this method provides accurate estimates of event durations and small uncertainties on the absolute age of the upper ice sections.However, the cumulative nature of the errors, associated with the increasing number of counted layers, leads to a decrease of the accuracy of absolute age with depth.For instance, the GICC05 (Greenland Ice Core Chronology 2005) composite timescale for Greenland ice cores (Rasmussen et al., 2006;Seierstad et al., 2014;Svensson et al., 2008) is associated with a maximum counting error of only 45 yr at ∼ 8.2 kyr B1950 (before 1950, current era).This error increases progressively with depth, reaching more than 2500 yr at ∼ 60 kyr B1950.Annual layer counting techniques cannot be applied when the annual layer thickness is too small to be resolved visually, e.g., in ice cores from central East Antarctica.
-Use of absolute age markers in ice cores.Well-dated tephra layers identified in ice cores during the last millennia provide valuable constraints (e.g., Sigl et al., 2014).Beyond that period, absolute age markers are very scarce.The links between 10 Be peaks and welldated magnetic events (Raisbeck et al., 2007) have provided an age marker for the Laschamp event (Singer et al., 2009).Promising results have recently been obtained using radiochronologic dating tools (Aciego et al., 2011;Buizert et al., 2014;Dunbar et al., 2008).
-Orbital dating in ice cores.Because there are few absolute constraints in ice cores beyond 60 kyr B1950 (limit for the layer counting in the NGRIP ice core), orbital tuning is the most effective method to provide chronological constraints on an ice core's deepest sections.In the first orbital dating exercises, tie points were determined from the tuning of water isotopic records on insolation curves (e.g., Parrenin et al., 2004), which limits further investigations of polar climate relationships with orbital forcing.More recent chronologies tried to circumvent this assumption and focused on non-climatic orbital markers.Three complementary tracers are currently used: the δ 18 O of atmospheric O 2 (δ 18 O atm ) (e.g., Bender et al., 1994;Dreyfus et al., 2007), δO 2 / N 2 (e.g., Bender, 2002;Kawamura et al., 2007;Suwa and Bender, 2008) and the total air content (e.g., Raynaud et al., 2007).While the link between δ 18 O atm and precession is explained by variations in the water cycle of the low latitudes, relationships between δO 2 / N 2 , air content and local summer insolation are understood to arise from changes in the surface snow energy budget influencing its metamorphism.Without a precise understanding of mechanisms linking these tracers to their respective orbital targets, the associated uncertainties remain large, 6 kyr for δ 18 O atm and 3-4 kyr for δO 2 / N 2 and air content (Bazin et al., 2013(Bazin et al., , 2014;;Landais et al., 2012).
-Ice core record synchronization.Inter-ice core matching exercises are undertaken to transfer absolute or orbital dating information from one ice core to another one.It generally relies on the global synchroneity of changes in atmospheric composition (CO 2 , CH 4 concentration, and δ 18 O atm ) (Bender et al., 1994;Blunier and Brook, 2001;Monnin et al., 2004), the identification of volcanic sulfate spikes within a given area (Parrenin et al., 2012b;Severi et al., 2007) or the hypothesis of synchronous regional deposition of aerosols recorded as ice impurities (Seierstad et al., 2014).In the first case, limitations are associated with the smoothing of atmospheric composition changes through firn air diffusion.In the second case, mismatches may arise through incorrect identification of events in different ice cores.
-Correlation with other well-dated climatic records.In some cases, high-resolution calcite δ 18 O records and precise U / Th dates on speleothems have been used to adjust ice core chronologies (Barker et al., 2011;Buizert et al., 2015;Parrenin et al., 2007a ) and low-latitude speleothem δ 18 O records (mostly reflecting changes in regional atmospheric water cycle).A main limitation of this method lies in the validity of this assumption.
-Modeling of the sedimentation process: snow accumulation, snow densification into ice, air bubbles trapping and ice flow (Goujon et al., 2003;Huybrechts et al., 2007;Johnsen et al., 2001).Glaciological modeling provides a chronology derived from the estimate of the annual layer thickness and, therefore, leads to more realistic event durations when the accumulation history and thinning function are well constrained.A side product of glaciological modeling is the quantification of changes in surface accumulation rates and the quantification of the initial geographical origin of ice.This additional information is necessary to convert measurements of concentrations of chemical species in ice cores into deposition fluxes and to correct ice core records from upstream origin effects (e.g., EPICA community members, 2006;Röthlisberger et al., 2008).Caveats are caused by unknown parameters of such glaciological models, such as amplitude of accumulation change between glacial and interglacial periods, the basal melting or the vertical velocity profile, which have a growing influence at depth.
A common and optimal chronology for several ice cores can be built through the combination of several of these methods in the frame of a probabilistic approach.The first attempts used absolute and orbitally tuned age markers along one ice core to constrain the unknown parameters of an ice flow model (e.g., Parrenin et al., 2001Parrenin et al., , 2004;;Petit et al., 1999).This method had however several limitations.First, the uncertainties associated with the ice flow model could not be taken into account, resulting in underestimated uncertainties.Second, the stratigraphic links between ice cores were not exploited, each ice core was dated separately, resulting in inconsistent chronologies.
A new probabilistic approach based on a Bayesian framework was subsequently introduced.The first tool, Datice, was developed by Lemieux-Dudon et al. (2010a, b).It introduced modeling errors on three canonical glaciological quantities of the dating problem: the accumulation rate, the LID of air bubbles and the thinning function (i.e., ratio between the in situ annual layer thickness on its initial surface vertical thickness).This method starts from a priori (also called "background") scenarios for the three glaciological parameters corresponding to a prior chronology for each ice core.These scenarios, deduced from a modeling of the sedimentation process, are associated with an uncertainty related to the degree of confidence in these prior scenarios.A minimization approach is then applied to find the best compromise between the prior chronological information for each ice core as well as absolute and relative age markers in the ice and in the air phases.This approach has been validated through the Datice tool and applied to build the Antarctic Ice Cores Chronology 2012 (AICC2012), producing coherent ice and air timescales for five different ice cores (Bazin et al., 2013;Veres et al., 2013): EPICA Dome C (EDC), Vostok (VK), Talos Dome (TALDICE), EPICA Dronning Maud Land (EDML) and NorthGRIP (NGRIP).Further developments of Datice were performed to incorporate additional dating constraints such as the depth intervals with known durations and correlation of errors (Bazin et al., 2014).Datice provides an excellent reference for this Bayesian approach.Still, because Datice has been developed over a long-term period with a continuous effort in calculation optimization through methodological improvement, the final code is difficult to access for a non-expert and cannot easily be used as a community tool.We thus identified the need for an open and user-friendly program with a performance similar to Datice but more easily used and implemented by different users within the ice core community.
In this paper, we present a new probabilistic model, IceChrono1, based exactly on the same approach as Datice but with improvements and simplifications in the mathematical, numerical and programming aspects.We first detail the IceChrono1 methodology highlighting the differences to Datice (Sect.2).We then perform dating experiments described in Sect. 3 using IceChrono1.We first replicate the AICC2012 experiment, and perform four additional experiments to test new functionalities of IceChrono1.The results of these experiments are discussed in Sect. 4. We summarize our main findings in the conclusions, and describe perspectives for future developments of IceChrono in Sect. 5.

Method
The table in the Appendix B lists all symbols used to describe the IceChrono1 methodology.

Method summary
The true chronology of an ice core, i.e., the ice and air ages at any depth, is a function of three variables (also functions of the depth): the initial accumulation rate (the accumulation rate when and where the particle was at surface), the lock-in depth of the air and the thinning function (the ratio between the thickness of a layer in the ice core to the initial accumulation rate).This is what we call the forward model.These variables are unknown and to find our optimal chronology we estimate them based on -prior information about their values on each ice core -chronological observations, such as (see Fig. 1) the ice or air age at a certain depth, the time elapsed between two depths, the synchroneity between two ice or air depths within two different ice cores or the depth shift between synchronous ice and air depths within the same ice core.
All these different types of information, mathematically described as probability density functions (PDF), are assumed to be independent and are combined using a Bayesian framework to obtain posterior estimates of the three input variables (accumulation, LID and thinning) and of the resulting chronologies.Uncertainties on the prior estimates and on the observations are further assumed to be Gaussian and the forward model is linearized, which allow using the Levenberg-Marquardt (hereafter LM) algorithm (Levenberg, 1944;Marquardt, 1963) to solve this least-squares optimization problem.The philosophy of the method is similar to that of the Datice method (Lemieux-Dudon et al., 2010a, b).

The forward model
For each ice core denoted by its index k, given three quantities (initial accumulation rate, air LID and incompressible ice thinning function), the model computes at any depth z k the age for the ice matrix χ k (in years relative to the drilling date) and the age for the air/hydrates contained in the ice ψ k (which, for simplicity, are assumed unique; we do not consider the age distributions due to gas diffusion and progressive lock in with depth) using the following formulas: where z k is the depth in the ice core (0 is the surface), D k is the relative density of the material with respect to pure ice (treated as a known time series), a k is the initial accumulation rate (expressed in meters in ice equivalent per year), τ k is the incompressible thinning function, d k is the depth (the depth shift between synchronous events in the ice and air phases, taken as dependent on the air depth, by convention), l k is the LID (taken as dependent on the air depth, by convention), D firn k is the average relative density of the firn when depth z k was at LID (treated as a known time series, dependent on the air depth, by convention) and z ie k is the ice equivalent depth: Note that in our convention inherited from glaciology, accumulation a k is expressed in ice equivalent per year and the thinning function τ k is expressed with respect to this iceequivalent accumulation rate, hence the appearance of D k in Eq. ( 1).We used this convention because most of the ice flow models, which give prior estimates on the thinning function (see below), only consider pure ice; i.e., they assume the snow is instantaneously densified into ice when it falls at the surface of the ice sheet.
The first equation integrates along the depth axis the number of annual layers per meter from the surface.The second equation means that the air age at depth z k is equal to the ice age at depth z kd k .This is the definition of depth.The third equation means that if one virtually unthins a depth interval between an ice depth and the synchronous air depth, one gets the unthinned lock-in depth in ice equivalent (ULI-DIE; see Fig. 2).The right-hand side of the third equation is an approximation, since we assume that the incompressible thinning function in the firn is constant through time.But because the incompressible thinning function is usually very close to 1 in the firn (the ice flow is almost negligible in the first ∼ 100 m of the ice sheet), this assumption is almost verified.No other approximation is introduced in this set of equations.However, there are additional approximations in the sedimentation models that provide the values of a k , l k and τ k and we will discuss the errors linked to these approximations below.
Let us assume we have prior estimates of the accumulation (denoted a b k ), of the LID (denoted l b k ) and of the ice thinning function (denoted τ b k ).The prior information can come from a combination of models and data.Typically, ice flow models can give an estimate of the thinning function, ice δD and δ 18 O are often used to deduce accumulation rates (e.g., Parrenin et al., 2007b), and firn densification models or δ 15 N measurements in air bubbles are used to deduce the LID and the average firn relative density.We now define logarithmic correction functions: so that the forward model effectively takes as inputs the three logarithmic correction functions instead of the three raw glaciological quantities.This change of variable allows us to transform Jeffreys variables into Cartesian variables (Tarantola, 2005) so as to express our problem into a leastsquares problem and will allow us to reduce the number of variables to be inverted (see below).
The three prior quantities a b k , l b k and τ b k as well as the relative density D k are discretized onto the same fine-depth grid: z k,i (i being the index of the node), called the age equation grid, which may not be regular.This allows us to compute prior ice and air age scales.Then, we discretize the correction functions c a k , c l k and c τ k onto coarser grids, which also may not be regular and which depend on the age for a and l and on the depth for τ : ) i as the correction vectors.These grids are coarser than the age grid to limit the number of variables to be inverted.The correction functions are then transferred to the age equation grid using a linear interpolation (for a and l, this step requires using the prior age scale).This allows us to compute the three glaciological quantities a k , l k and τ k onto the age equation grid.We note X as the state vector, that is, the vector containing all the correction vectors on the coarse grids for all the ice core sites.
Equation (1) for the ice age is solved first.Then, Eq. (3) for depth is solved, which allows us to deduce the air age from Eq. (2).To solve Eq. ( 3), we first integrate D k /τ k from the surface down to every depth in the age equation grid; i.e., we have a correspondence table between real depths and unthinned-ice-equivalent (UIE) depths.Then, for every air real depth in the age equation grid, we obtain the air UIE depth from the correspondence table.We then subtract the second member of Eq. ( 3) from this air UIE depth to get the ice UIE depth.Finally, we use the correspondence table to obtain the ice real depth and the depth.When we need to compute the ice age, air age or depth at depths which are not nodes of the age equation grid (for example when comparing the model with observations, see below), we use a linear interpolation.

The cost function
In probabilistic terms, one combines different sources of information which are assumed to be independent (the prior and the observations) and one looks for the most probable scenario, i.e., the most probable X (Tarantola, 2005).In mathematical terms, it corresponds to multiplying the PDFs of the prior and of the observations, the result of this multiplication being called the likelihood function L. Here, we assumed the PDFs to be independent multivariate Gaussian distributions.The posterior likelihood function L can therefore be written as where J , the cost function, is a sum of least-squares terms, each corresponding to an independent multivariate Gaussian PDF.Maximizing the likelihood function therefore corresponds to minimizing the cost function.In this case, we assume the information on each ice core and on each ice core pair to be independent.The cost function can therefore be written as a sum of terms: where J k is the term related to ice core number k and J k,m is the term related to the ice core pair (k, m).J k measures the misfit of the model with respect to the prior and the observations for ice core number k.It is written as the sum of independent terms: where J a k , J l k and J τ k are the prior terms for respectively a, l and τ , J ih k is linked to ice-dated horizons, J ah k is linked to air-dated horizons, J ii k is linked to ice intervals with known duration, J ai k is linked to air intervals with known duration and J d k is linked to depth observations.J k,m measures the misfit of the model with respect to the stratigraphic links in between ice cores k and m.It is the sum of four independent terms: where J ii k,m is linked to ice-ice stratigraphic links, J aa k,m is linked to air-air stratigraphic links, J ia k,m is linked to ice-air We first describe the prior terms: where P a k , P l k , P τ k are the correlation matrices of the prior and where R a k , R l k , R τ k are the residual vectors: with (σ a k,i ) i , (σ l k,i ) i , and (σ τ k,i ) i the standard deviation vectors for respectively the prior accumulation, LID, and thinning function.These three terms J a k , J l k and J τ k bring the "glaciological constraints" of the problem given by the sedimentation models.For example, they ensure that the optimum values for a, l and τ will be close to the prior values and also that their rates of change with respect to depth will be close to the rates of change of the prior values.That means the model giving the prior scenario for the accumulation, the LID, and the thinning function should have a quantified error.The correlation matrices define the correlation of uncertainties of the correction functions at different depth levels or ages.In practice, the IceChrono code is flexible, and any standard deviation vector and correlation matrix can be prescribed.
The other terms of J k are simply comparisons to observations, assumed independent: where R is the residual vector, (m i ) i are the model realizations of the observations, (o i ) i are the observations, (σ i ) i are their standard deviations and P is their correlation matrix.These terms are given in more details in the Appendix.
Again, the IceChrono model is flexible, and whatever standard deviation vectors and correlation matrices can be prescribed.
For J ih k and J ah k , one prescribes depth, age and standard deviation of ice-and air-dated horizons and the correlation matrix.For J ii k and J ai k , one prescribes the top depth, bottom depth, duration and standard deviation of ice and air intervals with known durations and the correlation matrix.For J d k , one prescribes air depth, depth and standard deviation of depth markers and the correlation matrix.For J ii k,m , J aa k,m , J ia k,m , J ai k,m , one prescribes the depth in the first core, the depth in the second core and a standard deviation in years of ice-ice, air-air, ice-air or air-ice stratigraphic links and the correlation matrix.
The determination of the standard deviation vectors and of the correlation matrices of the prior (resp.the observations) can be a difficult problem which requires an in-depth analysis of the methods used to determine the prior (resp.the observations).

Optimization
We now try to find the X input vector that minimizes the cost function defined in Eq. ( 9).This least-squares optimization problem is solved using a standard LM algorithm (Levenberg, 1944;Marquardt, 1963).The LM algorithm iteratively converges toward a minimum of the cost function and stops when a convergence criteria is met.To converge, the LM algorithm uses the Jacobian of the observation operator (the operator which output all the residuals from the X vector), which is here calculated using a finite difference approach.The LM algorithm gives as a result the optimized values of the model input variables vector X opt (containing the correction functions defined on the coarse grids) as well as an approximate estimation of its a posteriori error covariance matrix C X : the algorithm approximates the model by its linear tangent around the solution.
We note G(X) as the forward model vector containing all the variables for which we want to compute the optimal values and the uncertainties (in particular the ice and air ages on the age equation grids).The optimal model is therefore G(X opt ).From the model Jacobian G at X opt , which is calculated using a numerical finite difference approach, and from C X , one can compute an approximated value of the error covariance matrix for the model as The diagonal elements of C G give the variance of the model's variables.In practice, the matrix C G is never completely formed to limit the memory size of the program.Only diagonal bloc elements corresponding to a particular variable (e.g., ice age) on a particular ice core are calculated.
To ensure that the LM algorithm converges toward a global minimum of the cost function, we initialize it with different initial conditions X 0 all taken randomly according to the prior density of probability.We then check that the LM algorithm always converges toward the same solution.

Programming aspects
IceChrono1 is entirely coded in the python v2 programming language.IceChrono1 both resolves the optimization problem and provides figures to display the results.The core of Geosci.Model Dev., 8, 1473-1492, 2015 www.geosci-model-dev.net/8/1473/2015/ the code is entirely separated from the experiment setup directory which also contains the results of the run and which is composed of general parameter files, a directory for each ice core (which contains the parameters and observations for the given ice core) and a directory for each ice core pair (which contains the observations for the given ice core pair).
Only a basic understanding of python and an understanding of the structure of the experiment setup directory are needed to run IceChrono1.A detailed tutorial on how to use the IceChrono1 software is also available within the code.The core of the code is about 1000 lines long (including white lines and comments) and is built using an objectoriented paradigm.In such an object oriented language, apart for the classical type of variables (integer, real, characters, etc.), one can define his or her own classes of objects, containing variables and functions.In IceChrono1, a class exists for the ice core objects.It contains the variables related to this ice core: the age equation grid, the correction function grids, the prior scenarios and their associated standard deviations and correlation matrices, the relative density profile, the correction functions, the observations and their associated standard deviations and correlation matrices and the resulting calculated variables (accumulation, LID and thinning, ice and air ages, depth, ice and air layer thickness, etc.).It also contains functions performing the following tasks: the initialization of the ice core (i.e., reading of the parameters, priors and observations), the calculation of the age model, the calculation of the residuals, the calculation of the forward model Jacobian, the calculation of the standard deviations, the construction of the figures (for ice age, air age, accumulation, LID, thinning, ice layer thickness and depth) and the saving of the results.A class also exists for the ice core pair objects.It contains all the stratigraphic links and their associated standard deviation and correlation matrices relative to this ice core pair.It also contains functions that perform the following tasks: the initialization of the ice core pair, the calculation of the residuals, and the construction of the figures (for ice-ice links, air-air links, ice-air links and air-ice links).The main program is kept as simple and straightforward as possible.
We used the LM algorithm as implemented in the leastsq function from the scipy.optimizelibrary, which also provides an automatic convergence criteria.It does not try to minimize directly the cost function but rather a residual vector, the components of this residual vector being supposed independent from each other and with a unit standard deviation.Inside each term of the cost function We allow defining a correlation matrix P so that the residuals can actually be correlated.We thus used a Cholesky decomposition of P, and a change of variable, to transform the residual vector into a vector composed of independent variables with unit standard deviation.The associated term of the cost function can now be written as that is, the residuals are now independent and with a unit standard deviation.

Dating experiments
In this section, we describe the setups and the results of various dating experiments.
In this AICC2012-VHR (AICC2012 very high resolution) experiment, IceChrono1 is tested on the exact same setup, using the same observations and the same definition of the prior information as Datice was used.Only one aspect is modified: in AICC2012, the prior correlation matrices for the thinning function and LID are supposed to have an across-diagonal Gaussian decrease in amplitude.This Gaussian shape leads to a too high correlation for neighboring points.As a consequence, these matrices are numerically very difficult to invert (with a poor conditional state).We therefore opt for correlation matrices with an across-diagonal linear decrease in amplitude.We use regular grids with a time step of 1 kyr for the accumulation and LID correction functions and with 1001 nodes for the thinning correction functions.For the age equation grid, we use a step of 0.55 m for EDC, 1 m for VK, 1 m for TALDICE, 1 m for EDML and 1 m for NGRIP.
The IceChrono1 AICC2012-VHR experiment results are provided as a Supplement.This also includes all the figures produced by the IceChrono1 model.

Sensitivity test to resolution
In this experiment AICC2012-V2HR (AICC2012 very very high resolution), we test the impact of the resolution of the correction functions on the results of the optimization.We use exactly the same setup as for the AICC2012-VHR experiment, except for the resolution of the correction functions which is twice as high: we use regular grids with a time step of 0.5 kyr for the accumulation and LID correction functions and with 2001 nodes for the thinning correction functions.
Figure 3 compares the AICC2012-V2HR and AICC2012-VHR experiments for the ice age.The differences are minor and do not exceed 60 yr.
Figure 4 compares the resulting chronologies for the EDC ice core using the IceChrono1 (experiment AICC2012-V2HR) and Datice codes, as well as their uncertainties.Note that we did not use the official AICC2012 uncertainties for the ice and air ages (Bazin et al., 2013;Veres et al., 2013), but rather the raw uncertainties as calculated by the Datice code.The difference in both chronologies is less than 200 yr for the last 600 kyr and less than 500 yr all along the record.The uncertainties also show comparable shapes.We remark a few noticeable differences, though.There is a high frequency noise in the air age difference, especially for the most recent periods.This is due to the fact that the δ 15 N data which are used to infer the LID contain high-frequency variations and that Datice uses a higher resolution than IceChrono for these periods (1.65 m for EDC, 2 m for VK, 3 m for EDML, 1 m for TALDICE and 1 m for NGRIP).For the last 60 kyr, the uncertainties are almost always smaller in Datice than in IceChrono for the ice and almost always smaller in IceChrono than in Datice for the air ages.For example, between 15 and 60 kyr B1950, the uncertainty on the air age often decreases below 1000 yr in IceChrono1 but never in Datice.Also, the uncertainties of Datice are smaller in the ice phase than in the air phase.

Test experiment for ice intervals with known durations
We now perform a simplified dating experiment to test the inclusion of ice intervals with known durations.This is a new type of observation introduced in IceChrono, not available in the initial version of Datice (Lemieux-Dudon et al., 2010a, b) and not used in the AICC2012 experiment.It has only recently been developed in Datice (Bazin et al., 2014).This experiment is called "test-ice-intervals-with-known-durations".The experiment only contains the NGRIP ice core.The resolution chosen for the correction functions is 200 yr for the accumulation and LID correction functions and 501 nodes for the thinning correction function.
The prior accumulation is deduced directly from the deuterium content of the ice using an exponential relationship: where a 0 = 0.256 ice m yr −1 and β = 0.0193.The accumulation covariance matrix is defined by a uniform variance value of 20 % and with a correlation matrix with an acrossdiagonal linear decrease duration of 4000 yr.The prior thinning is deduced using the so-called pseudo-steady assumption (Parrenin et al., 2006): IceChrono-Datice IceChrono confidence interval with µ = 11 % being the ratio between melting and accumulation assumed constant through time and with ω being the flux-shape function (Parrenin and Hindmarsh, 2007), defined using a Lliboutry analytical model (Lliboutry, 1979), with p = −0.61.The covariance matrix of the thinning function is defined as a variance, which is assumed linearly related to the ice-equivalent depth z ie : where H ie is the ice-equivalent ice thickness and k = 0.5.The value of a 0 , β, µ and p were adjusted to obtain a good fit with the layer-counted GICC05 (Svensson et al., 2008) age markers sampled every 5 kyr along the core.
Here we include only one type of observations, ice intervals with known durations, which result from the layercounted GICC05 chronology (Svensson et al., 2008, and references therein) covering the last 60 kyr.This information is sampled every 1 kyr.The standard deviation of each observation on every 1 kyr interval is taken as half the variation of the maximum counting error (MCE) on this interval.All observations are assumed independent, i.e., the error on one 1 kyr interval is not correlated with the error on another 1 kyr interval.
The NGRIP ice age results of IceChrono1 on this optimization experiment are plotted in Fig. 6.As expected, the timescale is tightly constrained by the GICC05 ice intervals with known durations down to 60 kyr.The uncertainty of the NGRIP prior (blue) and posterior (black) chronologies and the posterior confidence interval (grey area) when using independent 1 kyr long intervals with known durations from the GICC05 layercounted timescale (green rectangles; Svensson et al., 2008).The 1σ uncertainty of the posterior timescale is also shown on the left (pink).This figure is automatically produced by IceChrono1.
optimized timescale is smaller than half the MCE: for example, at 60 kyr B1950, the uncertainty of the optimized age scale is 230 yr while half the MCE is 1300 yr.This is because the MCE is cumulative while the uncertainties of the 1 kyr long intervals with known durations are assumed to be independent and thus tend to cancel out.Indeed, the squared standard deviation of the sum of independent Gaussian variables is the sum of the squared standard deviation of the independent Gaussian variables.Beyond 60 kyr, when there are no ice interval constraints anymore, the uncertainty of the optimized age scale then increases quickly to reach ∼ 7 kyr at the bottom of the core.

Test experiment for correlated errors of observations
We now perform a simplified dating experiment to test the inclusion of correlated uncertainties.This is a new functionality in IceChrono with respect to the first version of Datice (Lemieux-Dudon et al., 2010a, b) and which has not been used in the AICC2012 experiment.It has only recently been developed in Datice (Bazin et al., 2014).This experiment is called "test-ice-intervals-with-known-durationswith-correlation".Again, this experiment contains only the NGRIP ice core.The resolution chosen for the correction functions is 200 yr for the accumulation and LID correction functions and 501 nodes for the thinning correction function.The prior information is exactly the same as the previous experiment.Only one type of observations is included: ice intervals with known durations.Again, this information comes from the layer-counted timescale GICC05 (Svensson et al., 2008, and references therein) and is sampled every 1 kyr.The standard deviation of each observation on every 1 kyr interval is taken as half the variation of the MCE on this interval.Contrary to the previous experiment, the observations are not assumed independent: the correlation matrix is assumed to contain ones on its diagonal and 0.5 outside of it.
The NGRIP ice age results of IceChrono1 on this optimization experiment are plotted in Fig. 7.As expected, the timescale is well constrained by the GICC05 intervals with known durations of observations down to 60 kyr and explodes beyond this age.Contrary to the previous experiment, the uncertainties of the intervals with known durations do not cancel out since they are correlated.Indeed, the standard deviation of the sum of completely dependent Gaussian variable is the sum of the standard deviations of the completely dependent Gaussian variables.As a consequence, the a posteriori uncertainty of the optimized timescale increases to ∼ 900 yr at 60 kyr, comparable to a MCE equal to ∼ 2600 yr at this age.

Test experiment for air intervals with known durations
We now perform a simplified dating experiment to test the inclusion of air intervals with known durations, which is a new type of observation in IceChrono with respect to the first version of Datice (Lemieux-Dudon et al., 2010a, b), not used in AICC2012.This experiment is called test-air-intervals-withknown-durations.The experiment contains only the EDC ice core.The resolution chosen for the correction functions is 200 yr for the accumulation and LID correction functions and 501 nodes for the thinning correction function.We perform the dating experiment only down to a depth of 1400 m (age of ∼ 100 kyr B1950).The prior information is the same as for the AICC2012-VHR experiment.We include observations as air intervals with known durations.This information comes from the layer-counted timescale GICC05 of NGRIP (Svensson et al., 2008, and references therein) and is transferred onto EDC through a NGRIP-δ 18 O ice / EDC-CH 4 synchronization.Each air-dated interval is defined by the onsets of DO events (including the Younger Dryas/Preboreal -YD/PB -transition) down to DO12 (∼ 46.8 kyr B1950).The top and bottom depths of the air intervals with known durations are directly taken from fast CH 4 transitions (Loulergue et al., 2007, Table 2).The durations of the air intervals with known durations are directly taken from the GICC05 timescale (Svensson et al., 2008, Table 1).The standard deviation of each interval observation is taken as half the MCE on this interval.The observations are assumed uncorrelated: the correlation matrix is assumed to be equal to the identity matrix.We also add the GICC05 age of the CH 4 YD/PB transition as an airdated horizon.This is equivalent to using an air interval from the LID to the CH 4 YD/PB transition with known duration, chronologies and the posterior confidence interval (grey area) when using correlated 1 kyr long intervals with known durations from the GICC05 layer-counted timescale (green rectangles; Svensson et al., 2008).The 1σ uncertainty of the posterior timescale is also shown on the left (pink).This figure is automatically produced by IceChrono1.EDC prior (blue) and posterior (black) chronologies and the posterior confidence interval (grey area) when using independent CH4 intervals with known durations from the GICC05 layer-counted timescale (green rectangles; Svensson et al., 2008).The 1σ uncertainty of the posterior timescale is also shown on the left (pink).This figure is automatically produced by IceChrono1.
but it is numerically more stable since the LID is not fixed during the iterative optimization procedure.
The EDC air age results of IceChrono on this optimization experiment are plotted in Fig. 8.As expected, the timescale is well constrained by the GICC05 intervals with known du-rations of observations down to DO12 and explodes beyond this age.The uncertainty on the absolute age also tends to increase inside intervals with known durations.This is expected since the total duration of an interval is constrained, but inside this interval only the prior information constrains the variation of the age.This feature is more visible here than in the two previous experiments because of the uncertainty on the LID which has a direct impact on the uncertainty of the air age from Eq. (2).By comparison, the uncertainty of the accumulation rate only has an impact on the uncertainty of the derivative of the ice age from Eq. ( 1).

Test experiment for mixed ice-air stratigraphic links
We now perform a simplified dating experiment to test the inclusion of mixed ice-air stratigraphic links, which is a new type of observation in IceChrono1 with respect to the first version of Datice (Lemieux-Dudon et al., 2010a, b), not used in the AICC2012.This experiment is called test-mixed-strati.
The experiment contains only the EDC and NGRIP ice cores.
The resolution chosen for the correction functions is 200 yr for the accumulation and LID correction functions and 501 nodes for the thinning correction function.We perform the dating experiment only down to a depth of 1400 m for the EDC ice core (age of ∼ 100 kyr B1950).The prior information for both the NGRIP and EDC ice cores are the same as for the test-ice-dated-intervals and test-air-dated-intervals experiments respectively.Only one type of observations is included: mixed ice-air stratigraphic links.This information comes from the synchronization of the EDC CH 4 record with the NGRIP δ 18 O record.Each stratigraphic link is the onset of one DO event (including the Younger Dryas/Preboreal transition) down to DO12 (∼ 46.8 kyr B1950).The EDC CH 4 depths are taken from Loulergue et al. (2007, Table 2).The NGRIP δ 18 O depths are taken from Svensson et al. (2008, Table 1).The standard deviation of each observation is taken as 100 yr for simplicity.The observations are assumed uncorrelated: the correlation matrix is assumed to be equal to the identity matrix.
The NGRIP air age vs. EDC ice age results of IceChrono1 on this optimization experiment are plotted in Fig. 9.While the prior scenario is in poor agreement with stratigraphic observations the optimized chronology agrees with the observations within their confidence intervals.

Robustness to a change of resolution
It is important to assess whether the formulation of IceChrono1 is robust to a change of resolution: when the resolution increases, the simulations should converge toward a meaningful result.IceChrono1 uses two different types of grids to optimize the ice cores age scales: the age equation grids and the correction function grids.
The age equation grids are used to solve Eqs. ( 1), ( 2) and (3).Equation ( 2) is the value of the ice age function at a given depth, so it is clearly robust to a change of resolution.Equations (1) and (3) are integrals and are therefore also robust to a change of resolution.
Concerning the correction functions grids, we made two test experiments with different resolutions: AICC2012-VHR and AICC2012-V2HR.The fact that the AICC2012-VHR and AICC2012-V2HR experiments agree well indicates that the formulation of the optimization problem in IceChrono1 is robust to a change of resolution of the correction functions.

IceChrono-Datice comparison on the AICC2012 experiment
The fact that IceChrono1 and Datice agree fairly well on this AICC2012 experiment indicates that both codes, which have been developed independently using different programming languages and different numerical schemes, are correct.But one has to keep in mind that both codes are based on the same main principles, so this is not a confirmation of these principles.Concerning the differences in the posterior confidence intervals observed for the last 60 kyr, we note that Datice uses an approximated version of Eq. (3) and that some limitations regarding the calculation of the uncertainties in Datice are known (see SOM, p. 280-289, of Bazin et al., 2014).These limitations are not present in IceChrono1 and all the uncertainties are calculated using the same formula (Eq.20).Some limitations have been corrected in a more recent version of Datice (Bazin et al., 2014).During the last glacial period, there are many CH 4 Antarctica-NGRIP stratigraphic links with uncertainties of a few centuries.With the NGRIP chronology being tightly constrained to GICC05 within 50 yr, it is expected that the air age uncertainty at EDC will sometimes decrease below 1000 yr during this time period.The posterior uncertainty calculated by IceChrono1 is therefore consistent with the chronological information used, in contradiction to that calculated by Datice.

IceChrono new functionalities
IceChrono1 has four new functionalities with respect to the initial Datice model (Lemieux-Dudon et al., 2010a, b): -observations as ice intervals with known durations -observations with correlated uncertainties -observations as air intervals with known durations -observations as mixed ice-air stratigraphic links.
We performed several tests which indicate that these functionalities work correctly.Note that the first two of these functionalities have also been implemented in a recent version of Datice (Bazin et al., 2014).The observations of intervals with known durations will be useful for future dating efforts to account for the layer counting information.This information is generally correlated since counting errors on two different intervals using the same method could be biased in the same way.Mixed ice-air stratigraphic observations will also be useful when accounting for the Antarctica-CH 4 /Greenland-δ 18 O synchronization links.

IceChrono-Datice codes comparison
Although IceChrono1 follows an approach similar to that of Datice, there are mathematical, numerical and programming differences with the initial version of Datice (Lemieux-Dudon et al., 2010a, b) that we will detail below.
1. Datice does not solve Eq. ( 3) but an approximation of it: that is, τ is assumed constant between the synchronous ice and air depths.
2. There is only one depth grid per drilling in Datice while the correction functions in IceChrono1 are discretized on different, coarser grids than the age equation grid.This allows reducing significantly the number of variables to be inverted and therefore to decrease the computation time.These coarser grids were not necessary in Datice since the analytical gradients already reduce the computation time.
Geosci 4. In IceChrono1, we allow for mixed ice-air and air-ice stratigraphic links in between ice cores.This is new with respect to Datice.A concrete example of the use of mixed ice-air stratigraphic links could be the synchronization of Dansgaard-Oechger events recorded in the methane records in Antarctica and in the ice isotope records in Greenland, as illustrated in the testmixed-strati experiment.This is especially useful when methane records for some Greenland ice cores are not yet available at sufficient resolution.
5. IceChrono1 allows for ice and air intervals with known durations.The functionality of ice intervals with known durations has also recently been implemented in Datice (Bazin et al., 2014) to take into account the information from layer counting.An example of the use of air intervals with known durations could be the dating of one Antarctic ice core using its methane record synchronized to the NGRIP ice isotopic record dated by layer counting.
6. IceChrono1 allows for correlated errors in observations for every kind of observation, while the original version of Datice does not allow for correlated errors for observations (Bazin et al., 2013;Lemieux-Dudon et al., 2010a, b;Veres et al., 2013).Note however that a new version of Datice implements this new functionality (Bazin et al., 2014).
7. In IceChrono1, the Jacobians of the model and of the observation operator are computed numerically by a finite difference approach while they are computed analytically in Datice.The Jacobian of the observation operator is needed by the minimizer to find the optimal solution X opt and its uncertainty C X .When the optimal solution is found, the Jacobian of the model allows evaluating the uncertainty of the model C G through Eq. ( 20).In Datice, analytical expressions of the Jacobians with respect to X have been derived and these expressions are used to numerically compute the Jacobians for a particular X.In IceChrono, each component of the X vector is alternatively perturbed and the forward model is run to evaluate how the model G(X) and the observation operator are modified.In other words, the Jacobians are evaluated by a finite difference approach.While a numerical computation of the Jacobian leads to a slower computation time, it leads to a more flexible use of the model since, if one modifies the formulation of the cost function or of the model, one does not need to derive again analytical expressions for the Jacobians, which is a complex task.
8. IceChrono1 is coded in a simple, flexible and straightforward way using the object-oriented python language.
It is very simple in IceChrono1 to modify the parameters of the problem, e.g., the age equation grid and the correction vectors grids.By comparison, IceChrono1 is about 1000 lines long (including the construction of the figures) while Datice is about 30 000 lines long of fortran code (without any figure construction).The reduction of the code length in IceChrono is in particular due to (1) the use of the python programming language, which was not advanced enough at the time when the Datice code was developed; (2) the use of the python leastsq function, which automatically calculates the gradient of the observation operator and the posterior C X variance-covariance matrix; and (3) the use of a numerical gradient of the model G'.Datice also implements posterior diagnostics of the data assimilation system which make its code length larger (Bazin et al., 2014).The simplicity of IceChrono will make it easy to add new functionalities to the code.

Current limitations of IceChrono and possible perspectives
IceChrono1 is already a useful tool to define a common and optimized chronology for several ice cores all together.However, it has several limitations that we will discuss below.
1.All uncertainties are assumed Gaussian.While the Gaussian probability density functions are the most often encountered in science, they are not always appropriate.For example, the volcanic synchronization of ice cores (e.g., Parrenin et al., 2012b) is ambiguous: a volcanic peak in one ice core can sometimes be synchronized to several other volcanic peaks in the other ice core.IceChrono1 therefore assumes that the features' recognition has been performed without ambiguity and that the only uncertainty remaining is the exact match of the two features within their durations.
2. The forward dating model is assumed to be "almost linear" in the plausible vicinity of the solution.Further developments would be necessary to diagnose if this assumption is justified.
3. IceChrono1 is not appropriate for optimizing the chronology of many ice cores at a very high resolution: the computation time would be too high, due to the numerical finite differences approach to evaluate the Jacobian matrices.In practice this is not a problem as a high resolution is only necessary for recent periods.If, in the future, the need for a more efficient dating model appears, we could develop an analytical gradient for the forward model, as it is done in the Datice software.dating methods are the same).We might allow in the future to take into account the correlation of observations between ice cores and ice core pairs.5. IceChrono1 requires the need for external sedimentation models to infer prior scenarios and their uncertainties.In practice, these sedimentation models also need to be optimized using age observations.In a next step, we will incorporate sedimentation models directly into IceChrono.The uncertainties of these sedimentation models could be inferred automatically by comparing them to the age observations.
6.The stratigraphic observations between the different ice cores need to be derived externally and imported as stratigraphic observations into IceChrono1.This step also requires some prior knowledge about the sedimentation process.Therefore, it would be best to incorporate it directly into the IceChrono software.Automatic methods for ice core synchronization would eliminate this step, which is most of the time done visually, in a subjective, fastidious and undocumented way.
7. The layer counting dating of ice cores needs to be done externally and imported as intervals with known durations of observations into IceChrono1.Again, this step also requires prior knowledge about the sedimentation process (e.g., the typical annual layer thickness).Therefore, it would be best to incorporate it directly into the IceChrono software.An automatic method for layer counting has already been proposed (Winstrup et al., 2012).
8. The definition of realistic prior correlation matrices is a difficult issue which will be dealt with in details in future studies.

Conclusions
We have developed and made accessible a new open-source probabilistic model to produce a common and optimized age scale for several ice cores, taking into account modeling and observation information.The code is similar in scope to Datice but has mathematical, numerical and programming differences: IceChrono1 is simpler to use, more flexible to develop and more powerful than Datice, although it might be slower to run depending on the chosen resolution.When compared to an AICC2012-like experiment, IceChrono1 and Datice generally produce similar results, which is a confirmation of both codes but which is not a confirmation of their principles which are identical.There are some differences in the evaluation of AICC2012 uncertainties for the last glacial period, and IceChrono appears to be more consistent with the chronological information which have been used.We also tested four new features of IceChrono1 with respect to Datice: the use of ice intervals with known duration, the use of correlated observations, the use of air intervals with known durations, and the use mixed ice-air stratigraphic links.Although primarily built for ice cores, IceChrono1 can also be used to date other paleoclimatic archives like marine cores, lake cores, speleothems, etc.The flexibility of IceChrono now opens interesting perspectives: the allowance of inter-cores and inter-core-pair correlation, the coupling with sedimentation models, the coupling with automatic synchronization methods, and the coupling with automatic layer counting methods.These developments will be made available in future versions of IceChrono.

Geosci
where z ih k,i is the depth of the ith ice-dated horizon in the kth ice core, χ obs k,i is its age, σ ih k,i is its standard deviation and where P ih k is the correlation matrix.
where z ah k,i is the depth of the ith air-dated horizon in the k th ice core, ψ obs k,i is its age, σ ah k,i is its standard deviation and where P ah is the correlation matrix.
where z ii,t k,i and z ii,b k,i are the top and bottom depths of the ith ice-dated interval in the kth ice core, χ obs k,i is its duration, σ ii k,i is its standard deviation and where P ii k is the correlation matrix.
where z ai,t k,i and z ai,b k,i are the top and bottom depths of the ith air-dated interval in the kth ice core, ψ obs k,i is its duration, σ ii k,i is its standard deviation and where P ai k is the correlation matrix.
where z d k,i is the depth of the ith depth observation in the kth ice core, d obs k,i is its observed value, σ d k,i is its standard deviation and where P d k is the correlation matrix.
where z ii,1 k,m,i and z ii,2 k,m,i are the depths in the kth and mth ice cores of the ith ice-ice stratigraphic link in the (k, m) pair of ice cores, σ ii k,m,i is its standard deviation and where P ii k,m is the correlation matrix.where z aa,1 k,m,i and z aa,2 k,m,i are the depths in the kth and mth ice cores of the ith air-air stratigraphic link in the (k, m) pair of ice cores, σ aa k,m,i is its standard deviation and where P aa k,m is the correlation matrix.
where z ia,1 k,m,i and z ia,2 k,m,i are the depths in the kth and mth ice cores of the ith ice-air stratigraphic link in the (k, m) pair of ice cores, σ ia k,m,i is its standard deviation and where P ia k,m is the correlation matrix.

Figure 2 .
Figure 2. Scheme of different representations of a firn column.(Left) When the firn column is at surface, its height is equal to the lock-in depth.(Middle) If one virtually converts the firn column into cumulated ice-equivalent accumulations, one gets the unthinned lock-in depth in ice equivalent (ULIDIE).(Right) When the firn column is buried down into the ice sheet and encounters vertical thinning, its height decreases to depth.

Figure 3 .
Figure 3. Sensitivity of the AICC2012-like experiment to a change of resolution of the correction functions.The black line is the age difference between the AICC2012-V2HR and AICC2012-VHR experiments (see text).

Figure 4 .
Figure 4. Comparison of IceChrono and Datice chronologies of the EDC ice core in the AICC2012 experiment.IceChrono is the black plain line and its standard deviation is the black dashed line.The AICC2012 standard deviation is represented by the grey area.(top) Ice chronologies.(bottom) Air chronologies.

Figure 5 .
Figure 5.Comparison of IceChrono and Datice chronologies of the EDC ice core in the AICC2012 experiment, for the last 60 kyr.IceChrono is the black plain line and its standard deviation is the black dashed line.The AICC2012 standard deviation is represented by the grey area.(top) Ice chronologies.(bottom) Air chronologies.

Figure 8 .
Figure8.Experiment test-air-intervals-with-known-durations.EDC prior (blue) and posterior (black) chronologies and the posterior confidence interval (grey area) when using independent CH4 intervals with known durations from the GICC05 layer-counted timescale (green rectangles;Svensson et al., 2008).The 1σ uncertainty of the posterior timescale is also shown on the left (pink).This figure is automatically produced by IceChrono1.

4.
The age observations on different cores or on different ice core pairs are not always independent (because the www.geosci-model-dev.net/8 ,i and z ai,2 k,m,i are the depths in the kth and mth ice cores of the ith air-ice stratigraphic link in the (k, m) pair of ice cores, σ ai k,m,i is its standard deviation and where P ai k,m is the correlation matrix.Discretized depth z k for the solving the age equations t a k,i Discretized time for the accumulation correction function t l k,i Discretized time for the lock-in depth correction function z Thinning prior standard deviation vector X opt Optimized model state vector C X A posteriori covariance matrix of the state vector G' Model Jacobian C G A posteriori covariance matrix of the model
www.geosci-model-dev.net/8/1473/2015/Geosci.Model Dev., 8, 1473-1492, 2015 Uncertainties are assumed Gaussian on depth observations in IceChrono1 while they are assumed lognormal in Datice.In practice, if the uncertainty on the depth observation is small in front of the value of the depth observation, this should make little difference.

. Model Dev., 8, 1473-1492, 2015 www
Below, we describe in detail the terms of the cost function which are linked to observations.