Inferring paleo-accumulation records from ice-core data by an adjoint-based method : Application to James Ross Island ’ s ice core

Abstract. Ice cores contain a record of snow precipitation that includes information about past atmospheric circulation and mass imbalance in the polar regions. We present a novel approach to reconstruct a climatic record – by both optimally dating an ice core and deriving from it a detailed accumulation history – that uses an adjoint-based method. The motivation of our work is the recent application of phase-sensitive radar which measures the vertical velocity of an ice column. The velocity is dependent on the history of subsequent snow accumulation, compaction and compression; in our inverse formulation of this problem, measured vertical velocity profiles can be utilized directly, thereby reducing the uncertainty introduced by ice-flow modelling. We first apply our method to synthetic data in order to study its capability and the effect of noise and gaps in the age–depth observations. The method is then applied to the ice core retrieved from James Ross Island, Antarctica. We show that the method is robust and that the results depend on the quality of the age–depth observations and the derived flow regime around the core site. The method facilitates the incorporation of increasing detail provided by ice-core analysis together with observed full-depth velocity in order to construct a complete climatic record of the polar regions.


Introduction
A layer of snow deposited at some position and time on the surface of an ice sheet today has a burial rate and depth that is dependent on the weight and amount of snow accumulated above it, its metamorphosis to ice, and the flow of the ice sheet.When the ice core is extracted, the chemical constituents trapped in the ice contain information on the climatic conditions during the time of deposition.
Dating the ice core is paramount in order to interpret the climatic record.Different methods are used to establish the ice-core timeline, mainly annual layer counting, localization of dated volcanic eruptions in the record, comparison of chemical species with other ice cores, and ice-flow modelling.Commonly, the ice-core chronology is obtained by a combination of methods (e.g.Hammer et al., 1978;Parrenin et al., 2007a;Lemieux-Dudon et al., 2010).Annual layer counting and the localization of known volcanic events offer accurate dating, but their use is restricted to the top of the ice core or to sparse sections of the record.
Ice-core dating relies on ice-flow models in sections where direct chronological information is unavailable or incomplete.The assumption made is that the record of snow accumulation is related to past cloud temperatures, which, in turn, are derived from ice-core isotope analysis (e.g.Watanabe et al., 2003;Parrenin et al., 2007b).However, this assumption can not be used on coastal areas of Antarctica or Greenland (Bromwich and Robasky, 1993;Kapsner et al., 1995), where the main driver of change in accumulation is atmospheric circulation and not the thermodynamic equilibrium of the atmosphere.However, recent advances in ice-core data analysis provide detailed chronologies for ice cores located in coastal areas using layer counting (e.g.Rasmussen et al., 2006;van Ommen and Morgan, 2010) or transferring timelines from reference ice cores using stratigraphic links (e.g.Ruth et al., 2007;Lemieux-Dudon et al., 2010).The aim of this paper is to use a detailed chronology such those mentioned above, de-Published by Copernicus Publications on behalf of the European Geosciences Union.

C. Martín et al.: Inferring accumulation records from ice cores
rived independently of temperature reconstructions, in order to obtain a record of palaeo-accumulation.
Previous work has estimated palaeo-accumulation and unknown ice-flow properties in order to extend the age-depth record beyond the confines of layer counting or stratigraphic links.Typically, an inverse method based on direct search is applied (e.g.Dansgaard and Johnsen, 1969) or, more recently, on Monte Carlo techniques (e.g.Parrenin et al., 2001Parrenin et al., , 2007b;;EPICA community members, 2004).Here, our approach is to apply an adjoint-based method for constrained optimization, commonly employed in computational fluid dynamics (CFD) for inverse problems that are governed by partial differential equations (e.g.Petra and Stadler, 2011).
The motivation of our work is the increasing interest in climatic records extracted from high-resolution ice cores located in coastal areas of Antarctica and Greenland (e.g. van Ommen and Morgan, 2010;Mulvaney et al., 2012).We have taken advantage of the recent development of a phase-sensitive radar that allows the direct measurement of the present vertical velocity full depth along the ice core (Corr et al., 2002;Jenkins et al., 2006); this greatly reduces the uncertainty introduced by ice-flow modelling (Kingslake et al., 2014).In addition, recent advances in multi-ice-core multi-component analysis has improved the detail in ice-core chronology (e.g.Rasmussen et al., 2006;Ruth et al., 2007;Lemieux-Dudon et al., 2010).
We first describe the method in Sect. 2 and then its application to synthetic and existing data from James Ross Island ice-core data (Mulvaney et al., 2012(Mulvaney et al., , 2014) ) in Sect.3. At the time of writing this manuscript we do not have a suitable location where both velocity derived from phase-sensitive radar and a detailed ice-core timeline is available.For simplicity, we have assumed that the flow of ice follows a Dansgaard and Johnsen (1969) analytical approach for the synthetic data (Sect.3.1) and a shallow ice approximation (e.g.Hutter, 1983) combined with the measured compaction for the James Ross Island ice-core case (Sect.3.2).

Description of the problem: forward model
The combined result of snow accumulation, compaction and compression is that a particle of ice deposited in the surface of a summit of ice thickness H is transported with a velocity (u, v, w) within the ice.The age of a particle of ice A, represents the time elapsed for a particle of ice to get transported from the surface to a given position (x, y, z), where z represents the vertical.
The age of a particle of ice can be expressed with the pure advection equation (e.g.Hindmarsh, 2001), because the increase in the age of a particle of ice is identical to the increase in time, where D represents the material derivative.In this paper, we consider that the distribution of the age varies in time, from a given initial state A 0 to the present (t = T ), but we neglect the effect of horizontal advection and the local variation of thickness through time ( ∂H ∂t = 0).The age of ice can then be expressed as We assume that the spatial variation of the velocity only depends on rheology and not on time-dependent variables such as surface accumulation or ice thickness (e.g.Parrenin et al., 2007a).This assumption is supported by the perturbation analysis of Wilchinsky and Chugunov (1997) for isothermal ice.Considering the conservation of mass, the vertical velocity is where a is the surface accumulation rate, m the basal melt rate, and η is the shape function.For the sake of simplicity, we refer to the surface accumulation rate and basal melt rate as accumulation and melting, as there is no possible ambiguity in this paper.

Inverse model
The inverse problem can be formulated as finding the records of accumulation a(t) and melting m(t), as well as the initial distribution of age-depth A 0 (z), that minimize the difference between observations and model results: where γ R a , γ R m and γ R A 0 denote the weight of the Tikhonov regularization for accumulation, melting and initial agedepth distribution.A d is the measured age of ice, σ d its estimated uncertainty and P d is the point measurement operator for discrete data (Petra and Stadler, 2011), where N d is the number of point measurements and δ(z−z j ) is the Dirac delta function.a(t), m(t) and A 0 (z) are functions of time or space that we want to find; A depends on them and can be obtained by solving Eq. (1).J is a functional that depends on a(t), m(t), A 0 (z) and A(z, t).
We follow a Lagrangian approach (e.g.Tröltzsch, 2010) that includes a Lagrangian multiplier function λ in order to find the solution of the minimization problem (Eq. 3) constrained by the forward model (Eq.1).
L(a, m, A 0 ; A, λ) = J (a, m, A 0 ; A) The optimal solution is obtained by setting the variations (functional derivatives) of the Lagrangian L with respect to all variables to zero.The Lagrangian in Eq. ( 5) confirms that setting the functional derivative of the Lagrangian with respect to the Lagrangian multiplier to zero (δ λ L = 0) is equivalent to the forward problem, also known as state equation (Eq.1).The details of the derivation of the optimal system of equations are described in the Appendix.
Setting the variations of L with respect to the age (δ A L = 0) leads to the adjoint problem for λ, which in the strong form can be written as Finally, the so-called control equations are obtained by calculating the variations of L with respect to accumulation, melting and initial distributions of age and setting them to zero (g a = δ a L = 0, g m = δ m L = 0 and g A = δ A 0 L = 0):

Numerical solution
The system of equations composed the state equation (Eq.1), the adjoint equation (Eq.6) and the gradients (Eq.7) represents a complex non-linear problem.In order to solve the system, we use a steepest descent method (e.g.Petra and Stadler, 2011).
First we assume initial solutions for accumulation, melting and initial distribution of age (a k , m k and A k 0 ).We then solve the state equation (Eq. 1) and subsequently, using the inferred solution of age A k , solve the adjoint equation backward in time to obtain the solution of the Lagrangian multiplier λ k .
In order to update the solution for accumulation, melting, and initial distributions of age and temperature, we increase them in the direction of steepest descent of the Lagrangian.
where β a , β m are β A are parameters.
The values of β a , β m and β A are calculated each iteration using a backtracking line search algorithm (e.g.Dennis and Schnabel, 1996, Sect. 6.3.2).The method assumes sufficiently large initial values for β a , β m and β A , which are successively halved until the corresponding decrease in the objective function J (Eq. 3) fulfil the Armijo (1966) where x represents the state variables a, m or A 0 ; g x represents the respective gradients of the objective function g a , g m or g A ; β x represents the parameters β a , β m or β A ; and the constant c is fixed to c = 10 −2 .The values of γ R a , γ R m and γ R A 0 , are chosen by inspection of the L-curve (e.g.Hansen, 2000).The Lcurve is a plot of the norm of the regularization term ( ) versus the difference in age-depth between observation and model results ( for different values of the regularization parameter.The L-curve has two distinctive sections: for sufficiently small values of the regularization parameter, the differences between data and model dominate the objective function, and a gradual increase in the smoothing of the solution implies a gradual increase in the differences between data and model; conversely, for sufficiently large values of the regularization parameter, the regularization term dominates the objective function, and a gradual increase in the smoothing of the solution translates into sharper differences between data and model.The curve is typically L-shaped, hence its name.The regularization parameters are chosen so that their corresponding terms in the objective function sit on the L-curve between these two sections, that is, the corner of the L. Regarding the numerical techniques, we solve the state equation (Eq. 1) and the adjoint equation (Eq.6) by using an operator-combined semi-Lagrangian Crank-Nicholson algorithm (SLCN; Spiegelman and Katz, 2006) in an equispaced grid of elevation {z i }.
The second-order two-level semi-Lagrangian scheme (e.g.Staniforth and Côté, 1991) solves the advection terms as an

C. Martín et al.: Inferring accumulation records from ice cores
ordinary differential equation along the trajectory that connects from some take-off elevation z * at time t = (j ) t to a regular grid point z i at time t = (j + 1) t in the state equation, which is solved forward in time, and at time t = (j − 1) t in the adjoint equation, which is solved backward in time ( t → − t, w → −w).The Crank-Nicolson algorithm is used to discretize the diffusive terms with a central difference.
The discretized state equation (Eq. 1) and the adjoint equation (Eq.6 ), using the SLCN method, are respectively where δ denotes the gradient discrete operator and * denotes the value at the take-off position z * and time t = (j ) t.

Recovering accumulation from synthetic age-depth observations
In order to study the capability of our method to recover palaeo-accumulation from age-depth observations, we have generated synthetic observations by adding random noise to a reference solution.
The reference accumulation solution varies over 10 kyr as where τ 1 = 10 kyr and τ 2 = 1 kyr.Figure 1 shows the reference accumulation and age-depth.
In this section, we assume that there is no melting, and that the shape function (Eq.2) follows Dansgaard and Johnsen (1969), where z k is a parameter that we fix to z k = 0.3 H .In addition, we assume that the initial age distribution (A(t = 0)) is exact and that it is a steady-state solution of the velocity distribution at time t = 0.
We first check our model by applying it to the agedepth derived from the reference accumulation.The algorithm (Sec.2.3) converges in 100-200 iterations, depending on the initial guess and the stopping criteria.After several numerical tests we used a grid with 1 m of vertical resolution and 50-year time steps.In this simulation we assume perfectly sampled data, i.e. the age-depth distribution is available at the numerical grid points.The results are not reported here as the recovered accumulation and age-depth converges to the exact reference accumulation.
Then, we add a 1 % random noise to the reference agedepth solution and apply our method to recover accumulation and age-depth.That is, to each sample of the reference age-depth value A i we add an independent value sampled uniformly in the interval [−0.01A i , 0.01A i ].The results are shown in Fig. 1.The recovered age-depth is nearly identical to the reference age-depth (under 0.1 % difference), even in the deepest area of the record.On the other hand, the random noise introduces differences between recovered and reference accumulation: the difference increases with time from the present up to a value of about 10 %.
Next, we investigate the effect of increasing noise.We compare the recovered accumulation from synthetic data obtained by adding 1, 5 and 10 % random noise to the reference age-depth.The results are shown in Fig. 2. The solution deteriorates when random noise increases, but the recovered accumulation only loses information for small wavelengths.
In the final numerical experiment (Fig. 3), we study the effect of not having data in an area.In the simulation we assumed that there are no data in 400 m< z < 600 m and that there is 1 % random noise in the synthetic age-depth.We find that the effect of not having any data over a given depth range of the borehole on the inversion is a period of time (area highlighted in Fig. 3) over which the retrieved accumulation is mainly constrained by the initially assumed profile and the regularization applied.

Recovering palaeo-accumulation from James Ross
Island ice core In this section, we apply the method to the ice core extracted from James Ross Island, a large island off the south-east side and near the north-eastern extremity of the Antarctic Peninsula.The data used in this section and their climatological implications are discussed in Mulvaney et al. (2012) and Mulvaney et al. (2014).
Figure 4 shows the timeline for the James Ross Island ice core, which is provided by fixed markers derived from local and global volcanic events with unequal estimated uncertainty (details in Mulvaney et al., 2012, Table S1 in the Supplement).
We assume that the shape function (Eq.2) can be expressed as where ρ is the measured density at a given depth, ρ I = 910 kg m −3 is the reference density of ice, and η I is the shape function for ice.We further assume that η I follows a shallow ice approximation (e.g.Hutter, 1983), where s is the normalized depth (s = (H − z)/H ).In order to measure the sensitivity to the uncertainty in ice flow in the recovered age-depth distribution and accumulation record, we assume that the rheological index n varies in the range 2 ≤ n ≤ 5. We assume that there is no melting at the base of the dome, as a temperature of −8.5 • C was recorded at the bottom of the ice core.Also, because the time interval we are interested in is only 3 kyr, we avoid the problem of guessing A 0 by expanding the simulation time until 10 kyr before present.That way, our solution is not sensitive to the value of A 0 .Based on different numerical tests we use a grid with 0.3 m resolution and 10-year time step.
The recovered age-depth (Fig. 4) follows the confidence intervals of the observed measurements, and the uncertainty introduced by the ice-flow parameter n is smaller than 1 %, indiscernible in the figure.However, the recovered palaeoaccumulation (Fig. 5) roughly follows the stable water isotope (δD), which is a proxy for temperature, but the result is not clear due to the small number of age markers (12 markers in 3 kyr) and the uncertainty introduced by the ice-flow parameter n.The uncertainty increases with time from the present and it becomes about 20 % of the estimated value after 3 kyr.
Finally, to estimate the sensitivity of the method to agedepth uncertainty, we use a basic Monte Carlo method (e.g.Campolongo et al., 2000).The error propagation is obtained by calculating the standard deviation, at a given depth, from a set of 100 results for both accumulation and age-depth.For each one of the 100 simulations, we resample the measured age-depth observations assuming a normal distribution with standard deviation equal to one-third of the uncertainty estimation (so that 99 % of the samples lie within the range of observations).
The results are shown in Fig. 6.The estimation of uncertainty in observed age-depth is variable, from 1 % at certain depths to about 10 % in the deepest sections of the ice core.Similar to the uncertainty propagation introduced by the ice flow (included in the figure), uncertainty in accumulation increases with time from the present but reaches a maximum of about 15 % after 3000 years.

Discussion
The quality of the information recovered from an inverse problem is always related to the quality of the observed measurements and to the nature of the assumptions made in the forward model.In this paper, we have analysed the effect on palaeo-accumulation (recovered information) of random noise, gaps and uncertainty in the age-depth distribution (observed measurements).We also evaluate the propagation of uncertainty introduced by assumptions in the forward model, which follows the compaction and compression of ice with time (Fig. 5).That uncertainty in ice flow is simplified in this study as the uncertainty in the rheological index n.
The experiments with synthetic data (Sect.3.1) show that our method is robust: in perfectly sampled data (Fig. 2), it recovers the reference accumulation in the presence of random noise and, even when the increase in noise deteriorates the accumulation, it recovers time-averaged values as the details are first lost in the high temporal frequencies.In addition, when there are gaps in the data (Fig. 3), the retrieved profile is constrained by the regularization.We find that the smoothing in the accumulation transforms, after solving the forward Figure 5. Recovered accumulation (blue) from James Ross ice core observed age-depth and stable w ter isotope δD (black).Red symbols represent the uncertainty in age-depth at the points where absol markers where found (Mulvaney et al., 2012).The area in blue represents uncertainty in recovered ac mulation due to unknown ice-flow parameter n (2 ≤ n ≤ 5 ).The stable water isotope δD has not be used for the inversion and is plotted only for reference.

26
Figure 5. Recovered accumulation (blue) from James Ross ice core observed age-depth and stable water isotope δD (black).Red symbols represent the uncertainty in age-depth at the points where absolute markers were found (Mulvaney et al., 2012).The area in blue represents uncertainty in recovered accumulation due to unknown ice-flow parameter n (2 ≤ n ≤ 5 ).The stable water isotope δD has not been used for the inversion and is plotted only for reference.model, into a smooth age-depth (Fig. 3).In contrast, a typical polynomial interpolation of the age-depth leads to spurious accumulation records due to the non-linearity introduced by the age equation (Eq.1; Fudge et al., 2014).
The sensitivity analysis to the James Ross ice-core data (Fig. 6) shows that between 1 and 10 % of uncertainty in the observations results in similar uncertainty in palaeoaccumulation.On the other hand, the propagation of uncertainty from the unknown parameter in the ice-flow model has a stronger effect on the accumulation (about 20 %) than in the age-depth (under 1 %).
It is not the objective of this study to evaluate the use of modelling in ice-core dating, but this later result shows that, in our case study, there is a relatively small sensitivity (un- der 1 %) of age-depth to uncertainty in ice-flow modelling.Age-depth depends on accumulation, and in the ice cores where accumulation can be extracted from temperature reconstructions and age-depth is constrained by fixed markers, our test case supports the use of modelling in ice-core dating (e.g.Parrenin et al., 2007a;Ruth et al., 2007;Lemieux-Dudon et al., 2010).
We only consider simplified ice-flow shape functions (derived from Dansgaard and Johnsen (1969) in Sect.3.1 and shallow ice approximation in Sect.3.2) as they are easier to analyse and because our aim is to incorporate measured full-depth vertical velocities as they become available in the future.Our main assumption is that the spatial variability of the velocity depends on compaction and rheology, and the time dependence is proportional to the combined result of accumulation and surface thinning.This assumption is partially supported by the perturbation analysis of Wilchinsky and Chugunov (1997), which assumes isothermal ice; however we do not consider the effect of time evolution on firnification (e.g.Arthern and Wingham, 1998), temperature along the core (e.g.Nereson and Waddington, 2002), ice fabric (e.g.Martín and Gudmundsson, 2012), dome migration (e.g.Martín et al., 2009), or horizontal advection or thinning (e.g.Huybrechts et al., 2007).
More complete ice-flow models can be incorporated into the inverse problem.The method consists of adding the weak form of the system of equations for the velocity (ice-flow model) to the Lagrangian L in Eq. ( 5) (e.g.Petra and Stadler, 2011, Sect. 3).The challenge is that the data extracted from the ice core, which are constraining the climatic reconstruction, are one-dimensional but the ice-flow problem is threedimensional with a poorly constrained rheology.The main issue when inferring palaeo-accumulation is that the ice along the core could have been deposited at the surface at a different position to that of the actual ice-core surface position, and ice domes near the coast are areas where the gradients in surface accumulation can be high (Lenaerts et al., 2014).
We expect that the analysis of full-depth velocity data sets near ice rises (Kingslake et al., 2014) provides information on spatial variations of ice flow on ice domes and helps to reconstruct the palaeo-accumulation records.

Conclusions
The present study describes a novel application of an adjointbased method to derive palaeo-accumulation records from ice-core chronology.The method is designed for ice cores in coastal areas of Antarctica and Greenland because it does not rely on temperature reconstructions.
Our approach is to follow a Lagrangian method for constrained optimization, commonly employed in CFD inverse problems governed by partial differential equations.We apply our method to synthetic and existing data from the James Ross Island ice core.We show that the method is robust, as it retrieves palaeo-accumulation from data with random noise and it reacts well to gaps in the data.We also show that the quality of the palaeo-accumulation depends on the quality of the age-depth observations, uncertainty and sampling, as well as the uncertainty in ice flow.
Furthermore, in our case study on the James Ross Island ice core, where full-depth velocity measurements are not available, the largest uncertainty in the inferred palaeoaccumulation is a consequence of the unknown ice-flow parameters.It increases with time until a maximum of about 20 %.In addition, an uncertainty in the chronology between 1 and 10 % propagates into a similar uncertainty in palaeoaccumulation.The Lagrangian L (Eq. 5) has been defined so that its minimization is equivalent to the minimization of J (Eq. 3), the difference between observations and model results, and the forward model (Eq.1).In this section, we calculate the functional derivative of L with respect to the Lagrangian multiplier λ, to check that the minimization of L is constrained to verify the forward model (Eq.1); with respect to the age, to obtain the adjoint problem (Eq.6); and with respect to the control variable a, to calculate its gradient g a (Eq.7).The gradients for melting (g m ) and initial age distribution (g A ) in Eq. ( 7) can be obtained in a similar way.By definition, the functional derivative of L with respect to the Lagrangian multiplier λ in a direction λ is

Figure 1 .Figure 1 .Figure 2 .
Figure 1.Recovered accumulation (red, left panel) and age-depth distribution (red, right panel) from age-depth synthetic data (green, right panel).The synthetic age-depth distribution has been obtained by adding 1% uniform noise to the reference solution (blue right-panel) derived from the reference accumulation record (blue, left panel).Inset in the right panel zooms into the age-depth profile in an area that is highlighted both in the accumulation record and in the age-depth distribution.

Figure 3 .Figure 3 . 25 Figure 4 .
Figure 3. Recovered accumulation (red, left panel) and age-depth distribution (red, right panel) from age-depth synthetic data (green, right panel).The synthetic age-depth distribution has been obtained by adding 1% uniform noise to the reference solution (blue right-panel) derived from the reference accumulation record (blue, left panel), and considering a gap in the data (400m < z < 600m).Inset in the right panel zooms into the age-depth results where there is a gap in the age-depth data.

Figure 6 .Figure 6 .
Figure6.Sensitivity analysis to uncertainty in chronology (red) and shape function η (blue).The uncertainty in chronology is assumed to follow a normal distribution.The uncertainty in the shape function η is represented by the range 2 ≤ n ≤ 5 for the rheological index in the Glen flow law.


δ λ L(a, m, A 0 ; A, λ)( λ) = d d L(a, m, A 0 ; A, λ + λ) =0 ,where is a variable and λ is an arbitrary function.The result of setting to zero the functional derivative of L with respect to the Lagrangian multiplier isδ λ L(a, m, A 0 ; A, λ)( 0) [A(z, 0) − A 0 ] dz = 0,which is the weak form of the forward model (Eq.1).That is, because λ is an arbitrary function, the equation is true when the age distribution satisfies the forward model.Similarly, by setting to zero the functional derivative of L with respect to the age, we obtainδ A L(a, m, A 0 ; A, λ)( Ā) = t Ā + w ∂ z Ā) dz dt + H 0 λ(z, 0) Ā(z, 0) dz = 0.We then integrate by parts λ∂ t Ā with respect to time and λw ∂ z Ā with respect to z, 0) Ā(z, 0) dz = 0.After setting Ā(H, t) to zero, as the age-depth distribution is null at the surface, we obtainH 0 Ā λ + [A(z, T ) − A d ] t)λ(0, t)w(0)dt = 0,which is the weak form of the adjoint problem (Eq.6).Finally, we calculate the gradient of L with respect to the accumulation as δ a L(a, m, A 0 ; A, λ)( ā) = γ R integration by parts to the first term of the right-hand side, can be written as δ a L(a, m, A 0 ; A, λ)( dt + ā(T )∂ t a(T ) − ā(0)∂ t a(0) = 0, which is the weak form of the gradient defined in Eq. (7).