Articles | Volume 15, issue 7
The Cryosphere, 15, 3495–3505, 2021
The Cryosphere, 15, 3495–3505, 2021

Research article 27 Jul 2021

Research article | 27 Jul 2021

Reconstruction of annual accumulation rate on firn, synchronising H2O2 concentration data with an estimated temperature record

Reconstruction of annual accumulation rate on firn, synchronising H2O2 concentration data with an estimated temperature record
Jandyr M. Travassos1, Saulo S. Martins2, Mariusz Potocki4,5, and Jefferson C. Simões3,4 Jandyr M. Travassos et al.
  • 1Graduate Program in Geophysics, Universidade Federal do Pará (UFPA), Rua Augusto Corrêa n, 01, Belém, Pará, Brazil
  • 2Graduate Program in Geology, Universidade Federal Rural do Rio de Janeiro (UFRRJ), 465 km 7, Seropédica, Brazil
  • 3Centro Polar e Climático, Universidade Federal do Rio Grande do Sul (UFRGS), Porto Alegre, Brazil
  • 4Climate Change Institute, University of Maine, Orono, ME 04469, USA
  • 5School of Earth and Climate Sciences, University of Maine, Orono, ME 04469, USA

Correspondence: Saulo S. Martins (


This work deals with reconstructing firn layer thicknesses at the deposition time from the firn's observed thickness in ice cores, thus reconstructing the annual accumulation, yielding a timescale and an ice-core chronology. We employed a dynamic time warping algorithm to find an optimal, non-linear alignment between an H2O2 concentration data series from 98 m worth of ice cores of a borehole on the central ice divide of the Detroit Plateau, the Antarctic Peninsula, and an estimated local temperature time series. The viability and the physical reliability of the procedure are rooted in the robustness of the seasonal marker H2O2 in a high-accumulation context, which brought the entire borehole to within the operational life span of four Antarctic stations around the Antarctic Peninsula. The process was heavily based on numerical optimisation, producing a mathematically sound match between the two series to estimate the annual layering efficiently on the entire data section at once, being disposition-free. The results herein confirm a high annual accumulation rate of aN=2.8 m w.e./yr, which is of the same order of magnitude as and highly correlated with that of the Bruce Plateau and twice as large as that of the Gomez Plateau, 300 and 1200km further south, respectively.

1 Introduction

Ice cores provide a continuous record of climatic and environmental data series based on ice’s physical and chemical properties, reflecting past atmospheric composition and climatic variability, (e.g. Masson-Delmotte et al.2006). Snow deposited on the ice surface is gradually compressed into firn and ice, having the ability to preserve a very reliable climate record with a low risk of missing years, provided that the accumulation rate is sufficiently high. A vital issue in the palaeoclimatic reconstruction is dating the stratigraphic sequence through different techniques, including 1-D to 3-D flow models (Nye1952; Dansgaard and Johnsen1969; Gillet-Chaulet et al.2012; Passalacqua et al.2016); the counting cycles of seasonally varying quantities; reference horizons, most commonly layers of high concentrations of sulfuric acid related to volcanic events (Vinther et al.2006); and layer identification through peaks of radioactive isotopes (Vinther et al.2006; Cuffey and Paterson2010). Often those techniques are combined, e.g. incorporating stable water isotope δ18Oatm into an ice flow model (Capron et al.2010).

Hydrogen peroxide H2O2 is produced by photochemical reactions in natural waters exposed to solar irradiation, surficial and atmospheric. It is the most stable of the reactive oxygen species created in the atmosphere through a chemical reaction requiring ultraviolet light. A kinetic model has explained that 76.7 % of the variation in H2O2 concentrations is due to solar irradiance and temperature variation only (Sigg and Neftel1988). In particular that production in Antarctica has a pronounced regular seasonality resulting from cycles of complete darkness in midwinter to 24 h daylight in midsummer. This gives a phenomenological basis for quasi-sinusoidal variability in H2O2 atmospheric concentration with maxima occurring during the sunlit summer (Steig et al.2005; Frey et al.2006). The H2O2 is an exceptionally robust marker for ice cores at high-accumulation sites in Antarctica where post-depositional losses are minimised, resulting in excellent preservation of the records, with summer-to-winter ratios of over 5 (Sigg and Neftel1988; Hutterli et al.2003; Frey et al.2006).

The H2O2 concentration data come from ice cores extracted from borehole DP-07-1 drilled in December 2007 at the ice divide of the Detroit Plateau (DP), at 640507′′ S, 593842′′ W, 1930 m above sea level (m a.s.l.). DP-07-1 reveals well-resolved seasonal cycles of H2O2 concentration data in the context of a very high deposition rate (Potocki et al.2016). We take advantage of the observed strong seasonality in the H2O2 record to estimate a core timescale spanning the entire firn horizon. That is done by synchronising the concentration data to an estimated temperature time series at the borehole location.

The maxima of H2O2 production and surficial atmospheric temperature occur during the sunlit months of the austral summer, allowing us to seek a correlation between their respective maxima. They do not necessarily coincide, but they both occur during summertime; the time difference between them is a fraction of a year. A temperature record at the borehole location on the DP may be estimated by interpolating the historical temperature recordings from six Antarctic stations not too far from the DP: Bellingshausen, Esperanza, Faraday–Vernadsky, Marambio, O'Higgins and Rothera, which have almost continuous meteorological observations from the late 1950s. We have discarded Bellingshausen and O'Higgins, the first for being heavily biased by maritime conditions. The second is a relatively short record with a sizable gap in it, leaving us with four stations forming the vertices of a polygon having the DP within its perimeter. Only Marambio lies on the eastern side, which may imply some unknown bias towards the western temperature regime of the Antarctic Peninsula. Figure 1 shows the locations of the Antarctic stations on an outline of the northern Antarctica Peninsula.

Figure 1The four Antarctic Stations, Esperanza (ES), Marambio (MA), Faraday–Vernadsky (FV) and Rothera (RO), and the borehole at the DP, on the northern Antarctica Peninsula. The white arrow in the lower right corner inset shows the location of the DP on the peninsula. Both maps were modified from a pan-sharpened scene of the Landsat Image Mosaic of Antarctica (LIMA) by USGS (, last access: February 2020).

The synchronisation of the concentration data to a temperature series is warranted here due to the local accumulation rate, high enough to bring the entire firn horizon deposition period within reach of the four stations' operational span. Both data series independently follow the same seasonal variation and the passing of the years, albeit in their particular manners; the H2O2 concentration displays a frequency scaling with depth, a result from the accumulated vertical strain, whereas the temperature series has a uniform frequency behaviour. The frequency scaling reflects the gradual thinning of the annual firn layers, which manifests itself as a frequency chirp in the H2O2 concentration series.

We have allowed for the frequency scaling of the peroxide concentration series concerning the uniform frequency temperature content by resorting to dynamic time warping (DTW). DTW is a fast and efficient algorithm for finding an optimal alignment between two sequences through a non-linear warping of one onto the other along the time–depth axis (Rabiner et al.1978; Sakoe and Chiba1978). We have worked with standardised versions of the peroxide and temperature series, using the distance between them as a measure for their resemblance (Rabiner et al.1978; Sakoe and Chiba1978). Once this is optimally found, the peroxide series becomes warped onto the temperature series, allowing for the observed frequency scaling.

Notwithstanding DTW being associated with speech recognition (Rabiner et al.1978; Sakoe and Chiba1978; Gilbert et al.2010), it has proved to be useful in several other applications. These encompass handwriting recognition (Jayadevan et al.2009), image and shape matching (Wang and Gasser1997; Latecki et al.2007), analysis and classification of the land cover of remotely sensed images (Verbesselt et al.2010; Xue et al.2014), gene expression and protein structure (Criel and Tsiporkova2005; Legrand et al.2008), and even brain activity (Chaovalitwongse and Pardalos2008). Speech recognition has been used to detect layers in deep Greenland ice cores, using a hidden Markov model (Winstrup et al.2010).

This work shows that DTW is also particularly fit for compensating for the peroxide frequency scaling with depth, realigning it to a temperature data time series and, at the same time, quantifying their dissimilarities. We have used the constant spectral content of the temperature data series as a reference in the pairing transformation through mathematical optimisation, thus yielding an estimate of a relation of depth to time without human intervention. Moreover, the procedure has also confirmed a very high deposition rate for the entire firn horizon at the DP.

2 The data sets

We deal with two independent data sets, a H2O2 concentration from the 133 m deep borehole and a temperature time series estimated at the DP. We have also collected a record of the stable water isotope deuterium, which was not used due to its poor seasonal variability (Potocki et al.2016). The borehole yielded intact ice cores down to z=109.3±0.5m, from where brittle ice began. The borehole temperature was fairly stable at -14.2±0.1C at a depth of 10 m. Depths in the borehole are measured with the origin at the surface and the vertical z axis pointing downwards.

The temperature time series at the borehole location was estimated through an interpolation procedure on a data set of continuous temperature readings since 1 January 1970, at four Antarctic stations on the Antarctic Peninsula. We will show below that the entire firn layer was accumulated in a shorter period than the >45 years of the estimated temperature time series.

2.1 The H2O2 concentration data

The H2O2 concentration data were retrieved from the first 98 m of ice cores with high-resolution sampling, with an average of 36 samples/year. It is a robust seasonal signal, well preserved for the entire depth range of ice cores (Potocki et al.2016). As for other ice cores at high-accumulation sites across West Antarctica, it is possible to establish a timescale for the core through straightforward counting of the annual cycles (Sigg and Neftel1988; Frey et al.2006).

The H2O2 concentration record, C(z), has considerable noise content throughout, which has to be minimised, making its seasonal signal conspicuous. We produce a smooth data series 𝒞(z) by robust fitting on C(z) through a loess nonparametric method (Cleveland and Grosse1991). Fig. 2 shows both C(z) and 𝒞(z) in micromolar (µM) concentration. It is easy to see the seasonal signal in 𝒞(z) as well as the effect of the accumulated vertical strain with depth on the annual firn layers. The latter manifests itself as a gradual thinning of the annual firn layers.

Figure 2The grey dots are the raw data C(z) and the solid line is their smoothed version 𝒞(z), both expressed in micromolar (µM) concentration. For the sake of visualisation, we have omitted just two data points with concentrations of C(z)>20µM at depths approximately equal to 5 m.


Notwithstanding some residual noise left on 𝒞(z), straightforward counting of its peaks and troughs suggests the first 98 m was probably accumulated in its entirety from the beginning of the 1980s. Direct division of the total depth span by the number of peaks indicates a very high deposition rate at the DP, which we will address below.

2.2 Estimating a temperature time series at the Detroit Plateau

The four stations shown in Fig. 1 have distinct sampling frequencies of temperature, varying from 1 to 8 readings/day. We set the beginning of the estimated temperature record to 1 January 1970; from this date onward, all four stations have continuous temperature readings. The end of the record is set to 29 December 2010, 3 years after the core was drilled at the DP. These limits yield a period wide enough to encompass the entire deposition period of the firn horizon safely.

We interpolated the daily temperature time series from the four stations shown in Fig. 1 through Delaunay triangulation, having the DP borehole sea-level projection inside the convex hull formed by the station set. That is a linear interpolation weighted by the inverse of the horizontal distance between a given station and the borehole projection. It is noteworthy that all stations but Marambio are located on the occidental part of the peninsula; nevertheless this station shares the most significant weights with Esperanza. Some bias towards the western climate regime somewhat compensated for by Marambio is thus expected, which is a fact we have to live with anyway.

Only the maximum daily temperature reading at each station was used in the interpolation process. The sea-level-interpolated time series at the DP, T(t), is further corrected to the height of the DP at 1930 m a.s.l., with a lapse rate in temperature with altitude of −0.55C/100 m (Rolland2003). Even taking care to obtain the best temperature estimates from the interpolation process at the DP, the accuracy of a particular temperature estimate is not crucial to our results. We use only the location in time of given summertime peak temperatures for synchronisation.

We alleviated aliasing due to the temperature sampling by applying a 2 d low-pass filter to T(t), a series with 14 973 data points, far more than the 985 data points of C(z). We made the number of data points in 𝒯(t) similar to those in C(z) by decimating the former by 15×. Again we avoided aliasing and conspicuously reduced noise in 𝒯(t) by low-pass filtering the decimated data series, using an eighth-order Chebyshev filter. We compensated for the amplitude losses incurred throughout the conditioning process by a constant multiplicative gain, bringing the amplitudes of the filtered temperature time series somewhat back to the original levels of the unfiltered T(t). The multiplicative factor is estimated in successive time windows as the quotient of the envelope of the original T(t) by an envelope of the non-gained version of 𝒯(t). From now on 𝒯(t) will refer to the accumulated temperature time series.

Figure 3 shows the decimated T(t) and 𝒯(t), spanning 41 years. The time series T(t) is quite noisy as one would have expected it to be, but the 𝒯(t) series proves to be a powerful depiction of the annual summer–winter cycles. It has a smaller amplitude than that of T(t), which is hardly an issue here as we are not looking for individual temperature figures but rather a reliable counter for the passing of the years.

Figure 3The grey dots are the interpolated and decimated temperature time series T(t). The solid curve is its smoothed and gained counterpart 𝒯(t).


3 Results

3.1 Warping H2O2 concentration data onto the temperature series

Figures 2 and 3 conspicuously show that the 𝒞(z) and 𝒯(t) data series record the passing of the years through their annual cycles of peaks and troughs, summer to winter, respectively. Nevertheless the two data series record the annual cycles in distinct manners, the former against depth and the latter against time, their similar shapes suggesting we could employ a simple mapping procedure from depth to year of deposition to a standard variable related to time.

There are two issues to consider here: (i) 𝒞(z) and 𝒯(t) have their respective zeniths in a given summer on different dates, as they are distinct phenomena, and (ii) the shapes of the two data series conspicuously differ from each other in terms of their spectral characteristics as is quickly seen comparing Figs. 2 and 3. The first issue is efficiently dealt with as peaks will differ from each other within a fraction of a given summertime, a noise source one just needs to be aware of. The second point is more involved as 𝒯(t) is a function of time with a nearly constant frequency content throughout, whereas 𝒞(z) has a frequency scaling with depth, a chirp behaviour easily seen in Fig. 2. The latter results from the gradual thinning of the firn layers due to the weight of the overburden.

The two data series 𝒞(z) and 𝒯(t) are not directly comparable, being functions of depth and time. We can make them comparable though by using a standardising mapping procedure:

(1) C i C i ^ = 1 σ C i C i - C T i T i ^ = 1 σ T i T i - T ,

where 𝒞i≡𝒞(z) and 𝒯i≡𝒯(t), with i=1,,N. C and T are averages, and σ(𝒞i) and σ(𝒯i) are standard deviations. The two standardised series, C^ and T^, have the same number of data points and have zero mean with unit standard deviation. The standardisation process minimises eventual y-axis discrepancies between the two series, decreasing the possibility of misalignment by the DTW algorithm. The mapping (system labelled as Eq. 1) is invertible, allowing a return to the original values whenever needed.

We warp the series C^, call it the sample, onto the reference series, T^, allowing for layer thinning with depth in the sample. In applying DTW we construct a warp path W=w1,w2,,wK between sample and reference, where each path element wk is linked to the two series indexes i,i, for the N elements in C^ and T^, respectively. The warp path length W is bounded to NK2N-1 and subject to the criteria below.

  • Boundary conditions. The warp path starts and ends at the first and the last elements of the two sequences, w1=1,1 and wK=N,N, all elements considered.

  • Continuity. The warping procedure preserves the ordering of the two aligned sequences:

  • Monotonicity. The elements of W are monotonically spaced in the independent variable, thus preventing big jumps:


The process of warping the sample onto the reference series is carried out by seeking the path W, which yields the minimum distance:

(2) D W = 1 2 N min k = 1 K d w k , w k + 1 ,

where dwk,wk+1 is the distance between two contiguous elements. DW should attain its minimum when the sample is correctly warped onto the reference signal (Sakoe and Chiba1978). We do the DTW through an algorithm using a correlation optimised warping, or COW, which aligns the sample onto the reference by piecewise linear stretching and compression of the warping segments with variable lengths l (Nielsen et al.1998; Pravdova et al.2002; Tomasi et al.2004). An integer slack parameter limits the range of possible segments l, initially set to unity, 𝔰=1. The reconstructed sample is obtained by retaining only the highest values obtained for the cumulative correlation coefficient:

(3) ξ T ^ , C ^ = l T i ^ - T ^ C i ^ - C ^ M - 1 σ T i ^ σ C i ^ ,

where the summation is performed for each segment l with M points, T^ and C^ are averages, and σTi^ and σCi^ are standard deviations. The problem is solved by applying the COW algorithm to all N/l segments through dynamic programming (Nielsen et al.1998; Pravdova et al.2002; Tomasi et al.2004). A complete description of the DTW and COW algorithms is well beyond the scope of this work; the reader is kindly referred to the literature cited herein.

The analysis proceeds as follows: begin the process of warping C^ onto T^ with the two series aligned at their respective beginnings, the borehole bottom and 1 January 1970, respectively. Warp C^ and retain the value of the total distance DW. Discard the year 1970 on T^, which now begins on 1 January 1971, and repeat the warping procedure with the entire C^ record; retain the new value for the total distance DW. Continue moving forward to the beginning of the T^ record in 1-year steps, storing the values of Dw estimated at each iteration. Continue this process of advancing the beginning of the T^ in 1-year steps, monitoring the evolution of the estimated values of Dw.

We observed a decreasing trend in the estimates of Dw retained at each round of warping described above, which reached a conspicuous minimum with the beginning of the T^ series aligned on 1 January 1980. Further 1-year steps on the starting date of the temperature ensured an increasing trend with a faster pace. We stopped the 1-year-step warping process on the increasing branch of Dw 5 years after reaching its minimum. Figure 4 shows both the original and the warped versions of series C^, with the borehole bottom, aligned with T^ on 1 January 1980. The figure also shows the behaviour of Dw for the entire year span we have considered in our calculations.

Figure 4Panel (a) shows the unwrapped C^ and T^ series in standardised ordinates, with the index i=0 corresponding to the mouth of the core. Panel (b) shows the two warped series with their abscissas i mapped back to time, expressed in years beginning on 1 January 1980: 𝒯, 𝒞(t). In both panels, C^ is shown as a dotted curve and T^ is shown as a solid curve, with ordinates in arbitrary units. Panel (c) shows the behaviour of distance Dw for the year we have performed the wrappings.


Once T^ is warped onto C^, one can easily perform an inverse mapping to the original depths and time, with i=1,,Nt;z. With that, depths may be mapped onto time, directly yielding a borehole timescale, z=z(t). That is shown in the lower panel of Fig. 4 where C^ is plotted against deposition time in years. The conspicuous minimum on Dw suggests a quantitative error estimate of 1 year, significantly greater than any eventual difference between the time of occurrence of the peaks in C^ and T^ within a given year.

3.2 Estimating a borehole timescale and accumulation rate

A simple model of an ice sheet flow considers that as a year's snowfall moves downwards relatively to the surface during its burial process by subsequent deposition undergoing viscoplastic deformation, it becomes progressively thinner, extending laterally due to ice incompressibility. An increase in density ensues with depth as the snow slowly compacts itself into firn and from that into ice. One way to simplify the process is to express all lengths in water-equivalent units (m w.e.), thus allowing one to disregard the compaction of snow before the complete transformation to ice. Accordingly we present depths as zmw.e.=zρzρw, where ρ(z) and ρw are the ice cores' density estimates and the density of pure water, respectively.

Figure 5The solid line shows the annual accumulation rate estimates at the DP and the dash-dotted line gives their 11-year moving average. Use the right ordinate for the annual accumulation rate and the left ordinate for borehole depths.


We use the measured ρ(z) from the ice cores to estimate an empirical model of firn densification, which assumes the density change with depth is proportional to the deviation relative to the density of pure glacier ice ρice=0.91 g/cm3 (Cuffey and Paterson2010). The model may have two (Herron and Langway 1980) or even three (Ligtenberg et al.2011) distinct firn densification stages, spanning from the surface to the zone of pore close-off. The adopted model has one densification stage from the surface down to the last available density estimate at zρ(max)=64.5 m: ρz=0.339 z0.1853, with R2=0.97 (Travassos et al.2018). The density measurements beyond zρ(max) were accidentally lost; so we impose the density of glacier ice onto the core bottom, ρ(109 m)=ρice, bridging the data gap with a straight line linking the imposed value to the last-measured density. This extrapolation will result in some inaccuracies, but as at zρmax=64.5 m the power law has already reached its slowest increase rate with depth, it may be reasonable to assume they are relatively small. On the other hand, that allows for the transformation of length dimensions to metres water equivalent (m w.e.) for the entire borehole. We will refer to this issue below whenever appropriate.

In the simplest model for an ice sheet flow, the total vertical strain of any layer is equal to the total vertical strain of the ice beneath it:

(4) λ z λ 0 = 1 - z h ,

where λ(z) is a layer of thickness at a depth z, which had accumulated as an annual layer of thickness λ0, and h is the total ice thickness, from the surface to the bedrock. The model considers a steady-state viscoplastic deformation with depth at the centre of an ice sheet, as the annual layers are buried by subsequent deposition. From now on, all length dimensions are in metres water equivalent (m w.e.), unless explicitly said otherwise. As the ice sheet is steady, we assume that accumulation and vertical thinning are constant in time and that a layer thickness does not vary horizontally. If those assumptions hold, the distance an ice particle moves downwards during 1 year must be equal to the thickness of one annual layer λ0.

As the older ice is closer to the bedrock, it is more convenient to express the vertical position of an ice particle concerning the rock bed interface using a new vertical axis, Z=h-z. The new coordinate frame runs in the opposite direction to the one we have been using so far, with z>0 pointing downwards. Assuming a steady state, the distance an ice particle moves downwards in 1 year, or, for that matter, the vertical particle velocity ν(Z), is a linear function of Z, and therefore, the thinning rate dνdZ is constant. The velocity at the surface equals the accumulation rate νh=-a, and at the bed ν(0)=0, the velocity being negative in the new reference frame as it points downwards:

(5) ν Z = - a Z h .

The relation between a given depth Z to the age of the ice is provided by

(6) t = h Z ν - 1 d Z Z = h exp - a h t ,

known as Nye's timescale (Nye1952, 1963; Cuffey and Paterson2010). The relation depicted in Eq. (6) provides the simplest model for describing how a layer of thickness λ0 deposited at the surface thins to λ(Z) when it is at a distance Z from the bedrock. Notwithstanding its simplicity, the Nye model still provides good estimates at shallow depths, close to the ones from more complex models, such as the Dansgaard–Johnsen model (Dansgaard and Johnsen1969; Cuffey and Paterson2010).

The warping of C^ onto T^ estimates the deposition thickness λ0 from its observed thickness λ(z), therefore reconstructing the accumulation as well as yielding a timescale z(t) spanning the entire borehole. The accumulation over the period 1980–2008, as revealed by the warped thicknesses λ0, shows wider oscillations towards later years. An 11-year moving average of accumulation shows a fairly stable regime for the period 1980–2008 of a11y2.5mw.e./yr. The small relative increase in accumulation from 1980–1990 to 1990–2008 seen in Fig. 5 is affected by the estimated densities deeper than 64.5 m used to transform depths. Moreover, the statistical significance of an 11-year moving average within a 28-year period is limited; we use it to compare with literature results, where the solar cycle period is often used.

We apply an exponential regression to the warped data to produce estimates for the two constants h,ah in the relation of Eq. (6). As the available data are confined to the firn layer, an estimate for the total thickness h is obviously beyond our reach.

Nevertheless, as the annual accumulation rate is assumed uniform, we can obtain an estimate for the 27 years before the coring activity of aN=2.82 m w.e./yr. Peak counting in Fig. 2 yields an estimated accumulation of ac=2.5mw.e./yr, equal to a figure reported elsewhere (Potocki et al.2016). The two accumulation rate estimates, aN and a11y, differ by ∼10 %, demonstrating reasonably compatibility considering the assumptions leading to the relation of Eq. (6) and providing a weak check on our numerical procedure.

Goodwin et al. (2016)Thompson et al. (1994)Thomas et al. (2008)

Table 1Location of third-party ice cores sites on the Antarctic Peninsula with their distances to the DP ice core. zmax is the maximum depth, and the period ΔT, in years, is shown in parentheses.

Download Print Version | Download XLSX

Figure 6Annual snow accumulation in ice cores from the DP (solid line), Bruce Plateau (dashed line), Gomez Plateau (dash-dotted line) and Dyer Plateau (dotted line) for the period 1980–2010.


It is worthwhile ending this section by comparing our estimated annual accumulation variability with data from the three ice cores listed in Table 1, all south of the DP in the Antarctic Peninsula. Figure 6 shows that the accumulation rates at the DP and Bruce Plateau are compatible throughout, an indication that both sites may have been subject to similar high-accumulation regimes, twice as large as that of the Gomez Plateau. Figure 6 also suggests annual snow accumulation for the period 1980–2010, giving a stable accumulation for all four ice cores. Nevertheless, the period spanned by our data is too short for probing multi-decadal trends; it has been reported that the Antarctic Peninsula has been experiencing an increased rate since 1900 (Thomas et al.2017). In particular, the Bruce Plateau ice core suggests an increase in snow accumulation during the late twentieth century, increasing at a rate of 0.19mmw.e./yr since the 1950s (Goodwin et al.2016).

4 Conclusions

Stratigraphic dating of ice cores is rooted in the use of reference horizons and annually resolved data to count annual layers to establish a core chronology. The latter uses outward data, e.g. volcanic events, to measure annual layers. This work has resorted to an independent data set, recorded temperature series, as a time reference to reconstruct a given layer thickness λ0 at deposition time from its observed thickness λ(z), thus reconstructing the annual accumulation and thereby a timescale or an ice-core chronology z(t).

The adopted non-linear numerical algorithm warped the H2O2 concentration data from borehole DP-07-1 onto an estimated local temperature record by aligning their respective summertime peaks, an interannual process with a ≃0.5-year time accuracy. Both the viability and the physical reliability of the procedure were rooted in the robustness of H2O2 as a seasonal marker associated with the observed high accumulation rate, which brought the entire borehole to within the operational life span of the Antarctic stations.

The considerable noise content in both series was alleviated through a nonparametric loess filter, which produced clean, smoothed versions of the data series albeit still retaining their complexity, as seen in Figs. 2 and 3. Any time difference between the summertime temperature and peroxide concentration peaks falls necessarily within the interannual process’ time accuracy of ≃0.5 years. The whole process was based on numerical optimisation, producing a mathematically sound match between the two series.

The secular variation in accumulation has revealed a high annual accumulation rate of aN=2.8 m w.e./yr, with the large variability seen in Fig. 5. The high accumulation rate observed at the DP is of the same order as the one reported for the Bruce Plateau, and they are highly correlated throughout the observational period considered here. The DP regime occurs 1 year earlier than at Bruce in a couple of time sections in Fig. 6, a small but detectable discrepancy, probably related to the distinct dating approaches. The conspicuous correlation of the DP and Bruce is an indication that the northern tip of the Antarctic Peninsula has been under a high-snow-accumulation regime, twice as large as that of the Gomez Plateau further south. The short period reported here is incapable of revealing multi-decadal trends; nevertheless, it is reasonable to suggest the DP may have been experiencing a similar increase in snow accumulation in the late twentieth century, similar to the one reported at the Bruce Plateau.

The limited time window of the period of our data reveals relatively stable behaviour throughout the 27 years before coring, with an 11-year moving average of the accumulation of a11y2.5mw.e./yr. A regularity in snow deposition preserved a reliable climate record, minimising post-depositional losses in the concentration of H2O2. We should expect a relatively short temporal range for firn layer ice cores in the northern Antarctic Peninsula by the same token, turning that region into a valuable climate record ranging through 3 decades before coring. The top DP layer should be now, almost 15 years after drilling, halfway through the firn layer if assuming deposition rate stability.

Mathematical procedures for annual layer counting are notoriously more laborious than manual counting; nevertheless, the latter has no other intrinsic quality but its easiness; quality or effectiveness cannot be technically guaranteed. As is the case of the present work, the former approach is indisputably rigorous, able to efficiently estimate the annual layering on the entire data section and disposition-free. The layer counting applied to our data produced annual accumulation figures that differ from those presented here by up to 40 %, with 17 % on average. All that considered, the choice ultimately remains with the investigator weighing in on an acceptable level of chronological inaccuracy in their work.

Comparison of algorithm results with simple layer counting performed on the smoothed versions of our data set suggests inaccuracies are non-uniform and within ±1 year. Notwithstanding that the algorithm is potentially usefully applied to other data sets where manual counting is more challenging than in the present case, it is not case-specific, and it is not restricted to the dyad peroxide temperature either; it can deal with other kinds of annually laminated data, not necessarily of a related origin, even among different wells. We are convinced there may be many different situations where there is the need to synchronise particular data sets where the procedure shown here may prove helpful.

Code availability

We have modified somewhat a COW code from (Tomasi et al.2004). The original code is available at

Data availability

The interpolated temperature daily series at the sea-level projection of the borehole DP-07-1 we have used is in the file “daily_temperature.asc” published at (Martins and Travassos2021). The peroxide concentration data can be requested from Mariusz Potocki (

Author contributions

JMT worked with the synchronisation of H2O2, the temperature data series and its accrued accumulation rate and wrote the manuscript. SSM estimated the temperature data series and contributed with additional data processing. MP processed the original H2O2 data series from the ice cores. JCS worked with all aspects of acquiring the ice cores in the field and contributed to several glaciological aspects of this work. All authors reviewed and agreed on the final manuscript.

Competing interests

The authors declare that they have no conflict of interest.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work was fully supported by the Brazilian Antarctic Program (PROANTAR) through the CNPq and CAPES. The present work is part of the ice-core programme Climate of Antarctica and South America (CASA) in association with the Climate Change Institute, University of Maine. The authors also acknowledge the Chilean Antarctic Institute (INACH), the Chilean Air Force (FACh), the Brazilian Air Force and the Brazilian Navy.

Financial support

This research has been supported by the INCT da Criosfera (CAPES project (grant no. 88887.136384/2017-00)) and the PROANTAR (CNPq project (grant no. 442755/2018-0)).

Review statement

This paper was edited by Michiel van den Broeke and reviewed by two anonymous referees.


Capron, E., Landais, A., Lemieux-Dudon, B., Schilt, A., Masson-Delmotte, V., Buiron, D., Chappellaz, J., Dahl-Jensen, D., Johnsen, S., Leuenberger, M., et al.: Synchronising EDML and NorthGRIP ice cores using δ18O of atmospheric oxygen (δ18O atm) and CH4 measurements over MIS5 (80–123 kyr), Quaternary Sci. Rev., 29, 222–234, 2010. a

Chaovalitwongse, W. and Pardalos, P.: On the time series support vector machine using dynamic time warping kernel for brain activity classification, Cybern. Syst. Anal., 44, 125–138, 2008. a

Cleveland, W. S. and Grosse, E.: Computational methods for local regression, Stat. Comput., 1, 47–62, 1991. a

Criel, J. and Tsiporkova, E.: Gene Time Expression Warper: a tool for alignment, template matching and visualization of gene expression time series, Bioinformatics, 22, 251–252, 2005. a

Cuffey, K. and Paterson, W.: The Physics of Glaciers, 4th edn., Academic Press, Elsevier,​​​​​​​ Burlington, MA 01803, USA, ISBN 9780123694614, 2010. a, b, c, d

Dansgaard, W. and Johnsen, S.: A flow model and a time scale for the ice core from Camp Century, Greenland, J. Glaciol., 8, 215–223, 1969. a, b

Frey, M. M., Bales, R. C., and McConnell, J. R.: Climate sensitivity of the century-scale hydrogen peroxide (H2O2) record preserved in 23 ice cores from West Antarctica, J. Geophys. Res.-Atmos., 111, D21301,,​​​​​​​ 2006. a, b, c

Gilbert, J. M., Rybchenko, S. I., Hofe, R., Ell, S. R., Fagan, M. J., Moore, R. K., and Green, P.: Isolated word recognition of silent speech using magnetic implants and sensors, Med. Eng. Phys., 32, 1189–1197, 2010. a

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576,, 2012. a

Goodwin, B. P., Mosley-Thompson, E., Wilson, A. B., Porter, S. E., and Sierra-Hernandez, M. R.: Accumulation variability in the Antarctic Peninsula: The role of large-scale atmospheric oscillations and their interactions, J. Climate, 29, 2579–2596, 2016. a, b

Herron, M. M. and Langway Jr., C. C.: Firn densification: an empirical model, J. Glaciol., 25, 373–385, 1980. a

Hutterli, M. A., McConnell, J. R., Bales, R. C., and Stewart, R. W.: Sensitivity of hydrogen peroxide (H2O2) and formaldehyde (HCHO) preservation in snow to changing environmental conditions: Implications for ice core records, J. Geophys. Res.-Atmos., 108, 4023,, 2003. a

Jayadevan, R., Kolhe, S., and Patil, P.​​​​​​​: Dynamic time warping based static hand printed signature verification, Journal of Pattern Recognition Research, 1, 52–65, available at: (last access: 17 July 2021), 2009. a

Latecki, L. J., Megalooikonomou, V., Wang, Q., and Yu, D.: An elastic partial shape matching technique, Pattern Recogn., 40, 3069–3080, 2007. a

Legrand, B., Chang, C., Ong, S., Neo, S.-Y., and Palanisamy, N.: Chromosome classification using dynamic time warping, Pattern Recogn. Lett., 29, 215–222, 2008. a

Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semi-empirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819,, 2011. a

Martins, S. and Travassos, J.: daily_temperature.asc, figshare [data set],, 2021. a

Masson-Delmotte, V., Dreyfus, G., Braconnot, P., Johnsen, S., Jouzel, J., Kageyama, M., Landais, A., Loutre, M.-F., Nouet, J., Parrenin, F., Raynaud, D., Stenni, B., and Tuenter, E.: Past temperature reconstructions from deep ice cores: relevance for future climate change, Clim. Past, 2, 145–165,, 2006. a

Nielsen, N.-P. V., Carstensen, J. M., and Smedsgaard, J.: Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping, J. Chromatogr. A, 805, 17–35, 1998. a, b

Nye, J.: The mechanics of glacier flow, J. Glaciol., 2, 82–93, 1952. a, b

Nye, J.: Correction factor for accumulation measured by the thickness of the annual layers in an ice sheet, J. Glaciol., 4, 785–788, 1963. a

Passalacqua, O., Gagliardini, O., Parrenin, F., Todd, J., Gillet-Chaulet, F., and Ritz, C.: Performance and applicability of a 2.5-D ice-flow model in the vicinity of a dome, Geosci. Model Dev., 9, 2301–2313,, 2016. a

Potocki, M., Mayewski, P. A., Kurbatov, A. V., Simoes, J. C., Dixon, D. A., Goodwin, I., Carleton, A. M., Handley, M. J., Jaña, R., and Korotkikh, E. V.: Recent increase in Antarctic Peninsula ice core uranium concentrations, Atmos. Environ., 140, 381–385,,​​​​​​​ 2016. a, b, c, d

Pravdova, V., Walczak, B., and Massart, D.: A comparison of two algorithms for warping of analytical signals, Anal. Chim. Acta, 456, 77–92, 2002. a, b

Rabiner, L., Rosenberg, A., and Levinson, S.: Considerations in dynamic time warping algorithms for discrete word recognition, IEEE T. Acoust. Speech, 26, 575–582, 1978. a, b, c

Rolland, C.: Spatial and seasonal variations of air temperature lapse rates in Alpine regions, J. Climate, 16, 1032–1046, 2003. a

Sakoe, H. and Chiba, S.: Dynamic programming algorithm optimization for spoken word recognition, IEEE T. Acoust. Speech, 26, 43–49, 1978. a, b, c, d

Sigg, A. and Neftel, A.: Seasonal variations in hydrogen peroxide in polar ice cores, Ann. Glaciol, 10, 157–162, 1988. a, b, c

Steig, E. J., Mayewski, P. A., Dixon, D. A., Kaspari, S. D., Frey, M. M., Schneider, D. P., Arcone, S. A., Hamilton, G. S., Spikes, V., Albert, M., Meese, D., Gow, A. J., Shuman, C. A., White, J. W. C., Sneed, S., Flaherty, J., and Wumkes, M.​​​​​​​: High-resolution ice cores from US ITASE (West Antarctica): Development and validation of chronologies and determination of precision and accuracy, Ann. Glaciol., 41, 77–84,, 2005. a

Thomas, E. R., Marshall, G. J., and McConnell, J. R.: A doubling in snow accumulation in the western Antarctic Peninsula since 1850, Geophys. Res. Lett., 35, L01706,,​​​​​​​ 2008. a

Thomas, E. R., van Wessem, J. M., Roberts, J., Isaksson, E., Schlosser, E., Fudge, T. J., Vallelonga, P., Medley, B., Lenaerts, J., Bertler, N., van den Broeke, M. R., Dixon, D. A., Frezzotti, M., Stenni, B., Curran, M., and Ekaykin, A. A.: Regional Antarctic snow accumulation over the past 1000 years, Clim. Past, 13, 1491–1513,, 2017. a

Thompson, L., Peel, D., Mosley-Thompson, E., Mulvaney, R., Dal, J., Lin, P., Davis, M., and Raymond, C.: Climate since AD 1510 on Dyer Plateau, Antarctic Peninsula: Evidence for recent climate change, Ann. Glaciol., 20, 420–426, 1994. a

Tomasi, G., Van Den Berg, F., and Andersson, C.: Correlation optimized warping and dynamic time warping as preprocessing methods for chromatographic data, J. Chemometr., 18, 231–241, 2004 (data available at:, last access: 12 July 2021). a, b, c

Travassos, J. M., Martins, S. S., Simões, J. C., and Mansur, W. J.: Radar diffraction horizons in snow and firn due to a surficial vertical transfer of mass, Brazilian Journal of Geophysics, 36, 507–518, 2018. a

Verbesselt, J., Hyndman, R., Newnham, G., and Culvenor, D.: Detecting trend and seasonal changes in satellite image time series, Remote Sens. Environ., 114, 106–115, 2010. a

Vinther, B. M., Clausen, H. B., Johnsen, S. J., Rasmussen, S. O., Andersen, K. K., Buchardt, S. L., Dahl-Jensen, D., Seierstad, I. K., Siggaard-Andersen, M.-L., Steffensen, J. P., Svensson, A., Olsen, J., and Heinemeier, J.​​​​​​​: A synchronized dating of three Greenland ice cores throughout the Holocene, J. Geophys. Res.-Atmos., 111, D13102,,​​​​​​​ 2006. a, b

Wang, K. and Gasser, T.: Alignment of curves by dynamic time warping, Ann. Stat., 25, 1251–1276,, 1997.  a

Winstrup, M., Svensson, A. M., Rasmussen, S. O., Winther, O., Steig, E. J., and Axelrod, A. E.: An automated approach for annual layer counting in ice cores, Clim. Past, 8, 1881–1895,, 2012. a

Xue, Z., Du, P., and Feng, L.: Phenology-driven land cover classification and trend analysis based on long-term remote sensing image series, IEEE J. Sel. Top. Appl., 7, 1142–1156, 2014. a

Short summary
This paper gives a timescale estimation and the yearly accumulation rate from ice cores encompassing the entire firn layer at the Detroit Plateau, the Antarctic Peninsula, through a non-linear pairing transformation of high-resolution H2O2 concentration data to a local temperature time series. An 11-year moving average of the yearly ice accumulation rate may suggest an increase in the span of 30 years, with an average of 2.5–2.8 m w.e./year.