Research article 27 Jul 2021
Research article  27 Jul 2021
Reconstruction of annual accumulation rate on firn, synchronising H_{2}O_{2} concentration data with an estimated temperature record
 ^{1}Graduate Program in Geophysics, Universidade Federal do Pará (UFPA), Rua Augusto Corrêa n, 01, Belém, Pará, Brazil
 ^{2}Graduate Program in Geology, Universidade Federal Rural do Rio de Janeiro (UFRRJ), 465 km 7, Seropédica, Brazil
 ^{3}Centro Polar e Climático, Universidade Federal do Rio Grande do Sul (UFRGS), Porto Alegre, Brazil
 ^{4}Climate Change Institute, University of Maine, Orono, ME 04469, USA
 ^{5}School of Earth and Climate Sciences, University of Maine, Orono, ME 04469, USA
 ^{1}Graduate Program in Geophysics, Universidade Federal do Pará (UFPA), Rua Augusto Corrêa n, 01, Belém, Pará, Brazil
 ^{2}Graduate Program in Geology, Universidade Federal Rural do Rio de Janeiro (UFRRJ), 465 km 7, Seropédica, Brazil
 ^{3}Centro Polar e Climático, Universidade Federal do Rio Grande do Sul (UFRGS), Porto Alegre, Brazil
 ^{4}Climate Change Institute, University of Maine, Orono, ME 04469, USA
 ^{5}School of Earth and Climate Sciences, University of Maine, Orono, ME 04469, USA
Correspondence: Saulo S. Martins (ssmsaulo@gmail.com)
Hide author detailsCorrespondence: Saulo S. Martins (ssmsaulo@gmail.com)
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 icecore chronology. We employed a dynamic time warping algorithm to find an optimal, nonlinear alignment between an H_{2}O_{2} 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 H_{2}O_{2} in a highaccumulation 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 dispositionfree. The results herein confirm a high annual accumulation rate of a_{N}=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 1200 km further south, respectively.
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. MassonDelmotte 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 1D to 3D flow models (Nye, 1952; Dansgaard and Johnsen, 1969; GilletChaulet 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 Paterson, 2010). Often those techniques are combined, e.g. incorporating stable water isotope δ^{18}O_{atm} into an ice flow model (Capron et al., 2010).
Hydrogen peroxide H_{2}O_{2} 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 H_{2}O_{2} concentrations is due to solar irradiance and temperature variation only (Sigg and Neftel, 1988). 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 quasisinusoidal variability in H_{2}O_{2} atmospheric concentration with maxima occurring during the sunlit summer (Steig et al., 2005; Frey et al., 2006). The H_{2}O_{2} is an exceptionally robust marker for ice cores at highaccumulation sites in Antarctica where postdepositional losses are minimised, resulting in excellent preservation of the records, with summertowinter ratios of over 5 (Sigg and Neftel, 1988; Hutterli et al., 2003; Frey et al., 2006).
The H_{2}O_{2} concentration data come from ice cores extracted from borehole DP071 drilled in December 2007 at the ice divide of the Detroit Plateau (DP), at $\mathrm{64}{}^{\circ}{\mathrm{05}}^{\prime}{\mathrm{07}}^{\prime \prime}$ S, $\mathrm{59}{}^{\circ}{\mathrm{38}}^{\prime}{\mathrm{42}}^{\prime \prime}$ W, 1930 m above sea level (m a.s.l.). DP071 reveals wellresolved seasonal cycles of H_{2}O_{2} 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 H_{2}O_{2} 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 H_{2}O_{2} 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.
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 H_{2}O_{2} 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 H_{2}O_{2} 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 nonlinear warping of one onto the other along the time–depth axis (Rabiner et al., 1978; Sakoe and Chiba, 1978). 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 Chiba, 1978). 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 Chiba, 1978; 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 Gasser, 1997; 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 Tsiporkova, 2005; Legrand et al., 2008), and even brain activity (Chaovalitwongse and Pardalos, 2008). 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.
We deal with two independent data sets, a H_{2}O_{2} 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=\mathrm{109.3}\pm \mathrm{0.5}$ m, from where brittle ice began. The borehole temperature was fairly stable at $\mathrm{14.2}\pm \mathrm{0.1}$ ^{∘}C 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 H_{2}O_{2} concentration data
The H_{2}O_{2} concentration data were retrieved from the first 98 m of ice cores with highresolution 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 highaccumulation sites across West Antarctica, it is possible to establish a timescale for the core through straightforward counting of the annual cycles (Sigg and Neftel, 1988; Frey et al., 2006).
The H_{2}O_{2} 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 Grosse, 1991). 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.
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 sealevel 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 sealevelinterpolated 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.55 ^{∘}C/100 m (Rolland, 2003). 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 lowpass 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 lowpass filtering the decimated data series, using an eighthorder 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 nongained 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.
3.1 Warping H_{2}O_{2} 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:
where 𝒞_{i}≡𝒞(z) and 𝒯_{i}≡𝒯(t), with $i=\mathrm{1},\mathrm{\dots},N$. $\stackrel{\mathrm{\u203e}}{\mathcal{C}}$ and $\stackrel{\mathrm{\u203e}}{\mathcal{T}}$ are averages, and σ(𝒞_{i}) and σ(𝒯_{i}) are standard deviations. The two standardised series, $\widehat{\mathcal{C}}$ and $\widehat{\mathcal{T}}$, have the same number of data points and have zero mean with unit standard deviation. The standardisation process minimises eventual yaxis 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 $\widehat{\mathcal{C}}$, call it the sample, onto the reference series, $\widehat{\mathcal{T}}$, allowing for layer thinning with depth in the sample. In applying DTW we construct a warp path $W=\left({w}_{\mathrm{1}},{w}_{\mathrm{2}},\mathrm{\dots},{w}_{K}\right)$ between sample and reference, where each path element w_{k} is linked to the two series indexes $\left(i,{i}^{\prime}\right)$, for the N elements in $\widehat{\mathcal{C}}$ and $\widehat{T}$, respectively. The warp path length W is bounded to $N\le K\le \mathrm{2}N\mathrm{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, ${w}_{\mathrm{1}}=\left(\mathrm{1},\mathrm{1}\right)$ and ${w}_{K}=\left(N,N\right)$, all elements considered.

Continuity. The warping procedure preserves the ordering of the two aligned sequences:
$$\begin{array}{rl}& {w}_{k}\left(i,{i}^{\prime}\right)\to {w}_{k+\mathrm{1}}\left(\widehat{i},\widehat{{i}^{\prime}}\right)\Rightarrow i\le \widehat{i}\le \left(i+\mathrm{1}\right)\\ & \mathrm{and}\phantom{\rule{0.33em}{0ex}}{i}^{\prime}\le \widehat{{i}^{\prime}}\le \left({i}^{\prime}+\mathrm{1}\right).\end{array}$$ 
Monotonicity. The elements of W are monotonically spaced in the independent variable, thus preventing big jumps:
$$\begin{array}{rl}& {w}_{k}\left(i,{i}^{\prime}\right)\to {w}_{k+\mathrm{1}}\left(\widehat{i},\widehat{{i}^{\prime}}\right)\Rightarrow \left(i\widehat{i}\right)\ge \mathrm{0}\\ & \mathrm{and}\phantom{\rule{0.25em}{0ex}}\left({i}^{\prime}\widehat{{i}^{\prime}}\right)\ge \mathrm{0}.\end{array}$$
The process of warping the sample onto the reference series is carried out by seeking the path W, which yields the minimum distance:
where $\mathrm{d}\left({w}_{k},{w}_{k+\mathrm{1}}\right)$ is the distance between two contiguous elements. D_{W} should attain its minimum when the sample is correctly warped onto the reference signal (Sakoe and Chiba, 1978). 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:
where the summation is performed for each segment l with M points, $\stackrel{\mathrm{\u203e}}{\widehat{\mathcal{T}}}$ and $\stackrel{\mathrm{\u203e}}{\widehat{\mathcal{C}}}$ are averages, and $\mathit{\sigma}\left(\widehat{{\mathcal{T}}_{{i}^{\prime}}}\right)$ and $\mathit{\sigma}\left(\widehat{{\mathcal{C}}_{i}}\right)$ 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 $\widehat{\mathcal{C}}$ onto $\widehat{\mathcal{T}}$ with the two series aligned at their respective beginnings, the borehole bottom and 1 January 1970, respectively. Warp $\widehat{C}$ and retain the value of the total distance D_{W}. Discard the year 1970 on $\widehat{\mathcal{T}}$, which now begins on 1 January 1971, and repeat the warping procedure with the entire $\widehat{C}$ record; retain the new value for the total distance D_{W}. Continue moving forward to the beginning of the $\widehat{T}$ record in 1year steps, storing the values of D_{w} estimated at each iteration. Continue this process of advancing the beginning of the $\widehat{\mathcal{T}}$ in 1year steps, monitoring the evolution of the estimated values of D_{w}.
We observed a decreasing trend in the estimates of D_{w} retained at each round of warping described above, which reached a conspicuous minimum with the beginning of the $\widehat{\mathcal{T}}$ series aligned on 1 January 1980. Further 1year steps on the starting date of the temperature ensured an increasing trend with a faster pace. We stopped the 1yearstep warping process on the increasing branch of D_{w} 5 years after reaching its minimum. Figure 4 shows both the original and the warped versions of series $\widehat{C}$, with the borehole bottom, aligned with $\widehat{T}$ on 1 January 1980. The figure also shows the behaviour of D_{w} for the entire year span we have considered in our calculations.
Once $\widehat{T}$ is warped onto $\widehat{C}$, one can easily perform an inverse mapping to the original depths and time, with $i=\mathrm{1},\mathrm{\dots},N\mapsto \left(t;z\right)$. 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 $\widehat{C}$ is plotted against deposition time in years. The conspicuous minimum on D_{w} suggests a quantitative error estimate of ≲1 year, significantly greater than any eventual difference between the time of occurrence of the peaks in $\widehat{C}$ and $\widehat{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 waterequivalent units (m w.e.), thus allowing one to disregard the compaction of snow before the complete transformation to ice. Accordingly we present depths as ${z}_{\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.}=\frac{z\phantom{\rule{0.125em}{0ex}}\mathit{\rho}\left(z\right)}{{\mathit{\rho}}_{\mathrm{w}}}$, where ρ(z) and ρ_{w} are the ice cores' density estimates and the density of pure water, respectively.
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/cm^{3} (Cuffey and Paterson, 2010). 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 closeoff. 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 z^{0.1853}, with R^{2}=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 lastmeasured density. This extrapolation will result in some inaccuracies, but as at ${z}_{\mathit{\rho}\left(max\right)}=\mathrm{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:
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 steadystate 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=hz$. 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 $\frac{\mathrm{d}\mathit{\nu}}{\mathrm{d}Z}$ is constant. The velocity at the surface equals the accumulation rate $\mathit{\nu}\left(h\right)=a$, and at the bed ν(0)=0, the velocity being negative in the new reference frame as it points downwards:
The relation between a given depth Z to the age of the ice is provided by
known as Nye's timescale (Nye, 1952, 1963; Cuffey and Paterson, 2010). 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 Johnsen, 1969; Cuffey and Paterson, 2010).
The warping of $\widehat{C}$ onto $\widehat{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 11year moving average of accumulation shows a fairly stable regime for the period 1980–2008 of ${\stackrel{\mathrm{\u203e}}{a}}_{\mathrm{11}\mathrm{y}}\cong \mathrm{2.5}$ $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}./\mathrm{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 11year moving average within a 28year 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 $\left(h,\frac{a}{h}\right)$ 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 a_{N}=2.82 m w.e./yr. Peak counting in Fig. 2 yields an estimated accumulation of a_{c}=2.5 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}./\mathrm{yr}$, equal to a figure reported elsewhere (Potocki et al., 2016). The two accumulation rate estimates, a_{N} and ${\stackrel{\mathrm{\u203e}}{a}}_{\mathrm{11}y}$, 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)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 highaccumulation 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 multidecadal 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.19 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}./\mathrm{yr}$ since the 1950s (Goodwin et al., 2016).
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 icecore chronology z(t).
The adopted nonlinear numerical algorithm warped the H_{2}O_{2} concentration data from borehole DP071 onto an estimated local temperature record by aligning their respective summertime peaks, an interannual process with a ≃0.5year time accuracy. Both the viability and the physical reliability of the procedure were rooted in the robustness of H_{2}O_{2} 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 a_{N}=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 highsnowaccumulation regime, twice as large as that of the Gomez Plateau further south. The short period reported here is incapable of revealing multidecadal 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 11year moving average of the accumulation of ${\stackrel{\mathrm{\u203e}}{a}}_{\mathrm{11}\mathrm{y}}\cong \mathrm{2.5}$ $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}./\mathrm{yr}$. A regularity in snow deposition preserved a reliable climate record, minimising postdepositional losses in the concentration of H_{2}O_{2}. 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 dispositionfree. 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 nonuniform 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 casespecific, 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.
We have modified somewhat a COW code from (Tomasi et al., 2004). The original code is available at http://models.kvl.dk/dtw_cow.
The interpolated temperature daily series at the sealevel projection of the borehole DP071 we have used is in the file “daily_temperature.asc” published at https://doi.org/10.6084/m9.figshare.14946177.v1 (Martins and Travassos, 2021). The peroxide concentration data can be requested from Mariusz Potocki (mariusz.potocki@maine.edu).
JMT worked with the synchronisation of H_{2}O_{2}, 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 H_{2}O_{2} 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.
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 icecore 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.
This research has been supported by the INCT da Criosfera (CAPES project (grant no. 88887.136384/201700)) and the PROANTAR (CNPq project (grant no. 442755/20180)).
This paper was edited by Michiel van den Broeke and reviewed by two anonymous referees.
Capron, E., Landais, A., LemieuxDudon, B., Schilt, A., MassonDelmotte, V., Buiron, D., Chappellaz, J., DahlJensen, D., Johnsen, S., Leuenberger, M., et al.: Synchronising EDML and NorthGRIP ice cores using δ^{18}O of atmospheric oxygen (δ^{18}O atm) and CH_{4} 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 centuryscale hydrogen peroxide (H_{2}O_{2}) record preserved in 23 ice cores from West Antarctica, J. Geophys. Res.Atmos., 111, D21301, https://doi.org/10.1029/2005JD006816, 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
GilletChaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sealevel rise from a newgeneration icesheet model, The Cryosphere, 6, 1561–1576, https://doi.org/10.5194/tc615612012, 2012. a
Goodwin, B. P., MosleyThompson, E., Wilson, A. B., Porter, S. E., and SierraHernandez, M. R.: Accumulation variability in the Antarctic Peninsula: The role of largescale 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 (H_{2}O_{2}) and formaldehyde (HCHO) preservation in snow to changing environmental conditions: Implications for ice core records, J. Geophys. Res.Atmos., 108, 4023, https://doi.org/10.1029/2002JD002528, 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: https://api.semanticscholar.org/CorpusID:9609636 (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 semiempirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819, https://doi.org/10.5194/tc58092011, 2011. a
Martins, S. and Travassos, J.: daily_temperature.asc, figshare [data set], https://doi.org/10.6084/m9.figshare.14946177.v1, 2021. a
MassonDelmotte, 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, https://doi.org/10.5194/cp21452006, 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., GilletChaulet, F., and Ritz, C.: Performance and applicability of a 2.5D iceflow model in the vicinity of a dome, Geosci. Model Dev., 9, 2301–2313, https://doi.org/10.5194/gmd923012016, 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, https://doi.org/10.1016/j.atmosenv.2016.06.010, 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.: Highresolution ice cores from US ITASE (West Antarctica): Development and validation of chronologies and determination of precision and accuracy, Ann. Glaciol., 41, 77–84, https://doi.org/10.3189/172756405781813311, 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, https://doi.org/10.1029/2007GL032529, 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, https://doi.org/10.5194/cp1314912017, 2017. a
Thompson, L., Peel, D., MosleyThompson, 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: http://models.kvl.dk/dtw_cow, 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., DahlJensen, D., Seierstad, I. K., SiggaardAndersen, 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, https://doi.org/10.1029/2005JD006921, 2006. a, b
Wang, K. and Gasser, T.: Alignment of curves by dynamic time warping, Ann. Stat., 25, 1251–1276, https://doi.org/10.1214/aos/1069362747, 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, https://doi.org/10.5194/cp818812012, 2012. a
Xue, Z., Du, P., and Feng, L.: Phenologydriven land cover classification and trend analysis based on longterm remote sensing image series, IEEE J. Sel. Top. Appl., 7, 1142–1156, 2014. a