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

This work deals with two distinct datasets, a well preserved hydrogen peroxide, H2O2, concentration data from firn cores at a high deposition location and a temperature time series, estimated from the daily records from four Antarctic stations around the Antarctic Peninsula. With them we have produced a time scale, an ice-core chronology, for the 133m deep borehole DP-07-1 from Plateau Detroit, Antarctic Peninsula. We constructed the chronology through a non-linear pairing transformation of the two series, based entirely on mathematical optimization, compensating the peroxide frequency scaling, reflecting the 5 gradual thinning of the annual firn layers with depth. We resorted to a dynamic time warping algorithm to find an optimal alignment between the two data series, allowing for the thinning of the annual firn layers with depth and the estimation of their original thicknesses at time of deposition. The core chronology spanning from Jan-1980 to Dec-2010 for the borehole reach, a time frame of a mere 30 years period, revealing a fairly stable 11 year average for the accumulation rate of 2.5m w.e./y.

This work shows that DTW is also particularly fit for compensating 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 60 content of the temperature data series as a reference in the pairing transformation through mathematical optimization, 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 PD. estimated at PD. 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 begun. The borehole temperature was fairly stable at −14.2 ± 0.1 • C at a depth of 10m, conductivity measurements on ice cores down to 20m had a modal value of 40.4µS/cm. Depths in the borehole are measured with the origin at the surface and the vertical z-axis pointing downwards. 70 The temperature time series at borehole location was estimated through an interpolation procedure on a data set of continuous temperature readings since January 1 st , 1970, at four Antarctic stations in Antarctic Peninsula. We are going to show below that the entire firn layer was accumulated in a shorter time span than the > 45y of estimated temperature time series.

The H 2 O 2 Concentration Data
The H 2 O 2 concentration data has a high resolution sampling, averaging to 36 samples/year along the 98m of ice cores. It 75 has a robust seasonal signal, well preserved for the entire depth range (Potocki et al., 2016). As for other ice cores at high accumulation sites across West Antarctica, it is possible to establish a time scale 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 a considerable noise content throughout, which has to be minimized making its seasonal signal conspicuous. We produce a smooth data series C (z) by robust fitting on C (z) through a loess nonparametric 80 method (Cleveland and Grosse, 1991). The Figure 2 shows both C (z) and C (z) in micro molar (µM) concentration. It is easy to see the seasonal signal in C (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 C (z) straightforward counting of its peaks and troughs suggests the first 100m were probably accumulated in its entirety from the beginning of the 80's. Straightforward division of the total depth span by 85 the number of peaks indicates a very high deposition rate at PD, a topic we are going to address below.

Estimating a Temperature Time Series at PD
The four stations shown in Figure 1 have distinct sampling frequencies in their respective temperature records: 3, 4, 8 and 1 readings /day, respectively. We set the beginning of the estimated temperature record to January 1 st , 1970, from this date onward all four stations have continuous temperature readings. The end of the record is set on December 29 th , 2010, three years after 90 the core was drilled at PD. These limits yield a time span wide enough to safely encompass the entire deposition period of the firn horizon.
We interpolated temperature time series from the four stations shown in Figure 1 through Delaunay triangulation, weighted by the inverse of the horizontal distance between a given station to the sea-level projection of the borehole. We have considered only the maximum daily temperature reading at each station. The sea-level interpolated time series at PD, T (t), is further corrected to the height of PD at 1930m asl, using a conservative lapse rate for the decreasing of temperature with altitude of −0.55 • C/100m (Rolland, 2003).
We alleviated aliasing due to the temperature sampling by applying a two-day low-pass filter to T (t), a series with 14973 data points, far more than the 985 data points of C (z). We made the number of data points in T (t) similar to those in C (z) by decimating the former by 15×. Again we avoided aliasing and conspicuously reduced noise in T (t) by low-pass filtering 100 the decimated data series, using an eight order Chebyshev filter. We compensated 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 envelope of the not gained version of T (t). From now on T (t) will refer to the gained temperature time series.  Figure 3 shows the decimated T (t) and T (t), spanning over 41 years. The time series T (t) is quite noisy as one would have expected it to be but the T (t) proves to be a robust 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 on the passing of the years.  2 and 3 conspicuously show that the C (z) and T (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 do 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 common variable related to time.

115
Two issues to consider here: (i) C (z) and T (t) have their respective zeniths at a given summer on distinct dates, as they are distinct phenomena, and (ii) the shape of the two data series conspicuously differ from each other in terms of their spectral characteristics as easily seen comparing Figures 2 and 3. The first issue is easily 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 (t) is a function of time with a nearly constant frequency content throughout, whereas C (z) has a frequency scaling with 120 depth, a chirp behavior easily seen in Figure 2. The latter results from the gradual thinning of the firn layers due to the weight of the overburden.
The two data series C (z) and T (t) are not directly comparable, being functions of depth and time. We can make them comparable though by using a standardizing mapping procedure, . . , N . The C and T are averages and σ (C i ) and σ (T i ) are standard deviations.
The two standardized series, C and T , have the same number of data points and are zero-mean with unit standard deviation.
The standardization process minimizes eventual y−axis discrepancies between the two series, dwindling the possibility of misalignment by the DTW algorithm. The mapping (1) is invertible, allowing to go back to the original values whenever needed.

130
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 = (w 1 , w 2 , . . . , w K ) between sample and reference, where each path element w k 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 N ≤ K ≤ 2N − 1 and subject to the criteria below.
-Boundary conditions: The warp path start and end at the first and the last elements of the two sequences, w 1 = (1, 1) and 135 w K = (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 done through seeking the path W which yields the minimumdistance, where d (w k , w k+1 ) is the distance between two contiguous elements. DW should attain its minimum when the sample is corrected warped onto the reference signal (Sakoe and Chiba, 1978). We do the DTW through with an algorithm using 145 a correlation optimized 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).
The range of possible segments l is limited by an integer slack parameter, initially set to unity, s = 1. The reconstructed sample is obtained by retaining only the highest values obtained for the cumulative correlation coefficient, 150 where the summation is performed for each segment l with M points, T and C are averages, and σ T i and σ C i are standard deviations. The problem is solved by applying the COW algorithm on 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.
Begin the process of warping C onto T with the two series aligned at their respective beginnings: the borehole bottom and year span we have considered in our calculations.
Once T is warped onto C one can easily perform an inverse mapping to the original depths and time, i = 1, . . . , N → (t; z).
With that depths may be mapped onto time, directely yielding a borehole time scale, z = z (t). That is is shown in the lower panel of Figure 4 where C is plotted against deposition time in years. The conspicuous minimum on D w suggests a quantitative error estimate of 1year, significantly greater than any eventual difference between the time of occurrence of the peaks in C 170 and T within a given year.

Estimating a Borehole Time Scale and an Accumulation Rate
A simple model of an ice sheet flow considers that as an year's snowfall moves downward relative to the surface during its burial process by subsequent deposition undergoing viscoplastic deformation, becoming progressively thinner, extending laterally due to ice incompressibility. An increase in density does ensue with depth as snow slowly compacts itself into firn and from that into 175 ice. One way to simplify the process is to express all lengths in water-equivalent units (mweq), thus allowing one to disregard the compaction of snow before the complete transformation to ice. Accordingly we present depths as z mweq = z ρ(z) /ρw, where ρ (z) and ρ w are the density measured in the ice cores from PD and 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.91g/cm 3 (Cuffey and 180 Paterson, 2010). The model may have two (Herron and Langway Jr, 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.5m: ρ 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 to the borehole bottom, ρ (109m) = ρ ice , bridging the data gap with a straight line linking the imposed value to the last measured 185 density. This extrapolation will result in some inaccuracies but as at z ρ(max) = 64.5m 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 mweq for the entire borehole. We will bring this issue back below whenever appropriate.
beneath it, of a layer of thickness λ (z) since it has been deposited at the surface as an annual layer λ 0 thick and h is the total ice thickness. The model considers a steady state viscoplastic deformation with depth at the center of an ice sheet, as the annual layers are buried by subsequent deposition. From now on all length dimensions are in mweq, unless explicitly said otherwise.

195
As ice sheet is in a steady state 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 one year must be equal to the thickness of one annual layer λ (z).
As the older ice is closer to the bedrock, it is more convenient to express the vertical position of an ice particle in relation to the rock bed interface using a new vertical axis, Z = h − z. The new coordinate frame runs in the opposite direction to the 200 one we have been using so far with z > 0 pointing downwards. Assuming a steady state the distance an ice particle moves downwards in one year, or for that matter, the particle vertical velocity ν (Z), is a linear function on 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, The relation between a given depth Z to the age of ice is given by known as Nye's time scale (Nye, 1952(Nye, , 1963Cuffey and Paterson, 2010). Relation (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 bed. Notwithstanding its simplicity the Nye model still provides good estimates at shallow depths which are close to the ones from more 210 complex models, such as the Dansgaard-Johnsen model (Dansgaard and Johnsen, 1969;Cuffey and Paterson, 2010).
The warping of C onto T estimates the deposition thickness λ 0 from its observed thickness λ (z), therefore yielding a time scale z(t) spanning the entire borehole. The z(t) is shown in Figure 5 is a spline linking the center of each corrected, warped, thickness λ 0 against time. The beginning of the extrapolated densities beyond 64.5m, shows as a subtle inflection at ∼ 1990.
The annual accumulation over the period 1980-2008 as revealed by the warped thicknesses λ 0 show wider oscillations 215 towards later years, but their 11-year moving average,ā 11y ∼ = 2.5 m w.e./y appears fairly stable for the period 1980-2008.
Figure 5 also shows the 11-year annual accumulation estimated at Gomez in the south-western Antarctic Peninsula (Thomas et al., 2008), showing that a ma is more than twice as large, a very high accumulation rate.  (Thomas et al., 2008).
Apply an exponential regression on the warped data to produce estimates for the two constants (h, a /h) in relation (6). As the available data is confined to the firn layer an estimate for the total thickness h is obviously beyond our reach. Nevertheless 220 as annual accumulation rate is assumed uniform we can obtain an estimate for the the 27 years prior to the coring activity, a N = 2.82 mweq /y. The two accumulation rate estimates, a N andā 11y are compatible, providing a check on our numerical procedure.

Conclusions
We have produced a time scale, an ice-core chronology z(t), for the 133m deep borehole DP-07-1 drilled in Plateau Detroit, Those two estimates may be regarded as defining a range for a very high annual accumulation rate. That range points to an https://doi.org/10.5194/tc-2020-342 Preprint. Discussion started: 5 January 2021 c Author(s) 2021. CC BY 4.0 License. accumulation rate more than 3 times the one reported at Gomez in the south-western Antarctic Peninsula 0.8 mweq /y (Thomas et al., 2008).
By the same token we have taken advantage of the strong seasonality on H 2 O 2 data and the high accumulation rate to 240 perform the warping procedure we could have used them to estimate the mean annual accumulation rate by straightforwardly counting the peroxide annual cycles. If one counted the peaks in Figure 2 one would have estimated an accumulation of a c = 2.5 m w.e./y as it is reported elsewhere (Potocki et al., 2016). On the other hand one would not have been able to retrieve the time-corrected layer thickness accumulated at the surface at a given year λ 0 . All the results related to the accumulation at PD were obtained through a mathematical optimization process, without human intervention.