A representative density profile of the North Greenland snowpack

. Along a traverse through North Greenland in May 2015 we collected snow cores up to 2 m depth and analyzed their density and water isotopic composition. A new sampling technique and an adapted algorithm for comparing data sets from different sites and aligning stratigraphic features are presented. We ﬁnd good agreement of the density layering in the snowpack over hundreds of kilometers, which allows the construction of a representative density proﬁle. The results are supported by an empirical statistical density model, which is used to generate sets of random proﬁles and validate the applied methods. Furthermore we are able to calculate annual accumulation rates, align melt layers and observe isotopic temperatures in the area back to 2010. Distinct relations of δ 18 O with both accumulation rate and density are deduced. Inter alia the depths of the 2012 melt layers and high-resolution densities are provided for applications in remote sensing.


Introduction
In the context of global warming, the Greenland ice sheet has been identified as a so-called "tipping point" of climate change (Lenton et al., 2008).The sea level rise caused by its decay may have a severe impact on human society as well as ecological systems.Thus the difference in accumulation across the interior of the ice sheet and seasonal melting, runoff and calving at its borders, the so-called mass balance, has been in the focus of recent scientific activities in the Arctic region.The applied methods for its determination range from satellite remote sensing (e.g., Zwally et al., 2011) to regional climate modeling (e.g., Fettweis, 2007) and to large-scale climate simulations constrained by weather station data and ice core records (e.g., Hanna et al., 2011).Even though first accumulation and density measurements had already been carried out in [1952][1953][1954] (Bull, 1958) using accumulation stakes and Rammsonde measurements at a few points alongside the gravity survey of the British North Greenland Expedition, large-scale studies such as Benson (1962) are still very rare.To obtain accumulation maps of Greenland such as Bales et al. (2009), diverse data sets from ice cores, snow pits and weather stations have to be collected over several years.Recently Hawley et al. (2014) conducted a ground-penetrating radar survey alongside a traverse of about 1000 km length, supported by a few snow pits and shallow cores for bulk densities and chemical profiling.Koenig et al. (2016) used airborne snow radar to determine accumulation rates from 2009 to 2012 along flight paths of more than 10 000 km.
In summer 2012, there were 2 very warm days with temperatures above 0 • C almost all over Greenland, causing substantial melt layers (Nghiem et al., 2012).Although this was a very rare event induced by a special weather situation (Bennartz et al., 2013), the newly formed ice layers strongly influenced the physical properties of snow and firn (Nilsson et al., 2015).
We introduce a new and efficient technique for sampling the snowpack along traverses, which allows for additional lab-based measurements to gain high-resolution profiles of physical snow properties such as density.Furthermore we adapt an algorithm from speech recognition to align these spatially distributed data sets and provide further insight into Published by Copernicus Publications on behalf of the European Geosciences Union.
Table 1.Measurement sites along the traverse, see also Fig. 1.The missing liner numbers (e.g., N2E_01) result from multiple samples being taken at some locations.Nonetheless, only one profile per location was used for this study.their development with changing surrounding conditions.The method is tested with randomly generated sets of density profiles with the same statistical properties as the original measurements.As an application we present data gained along a 450 km traverse in North Greenland, deduce relations of the individual parameters (density, δ 18 O and accumulation rate) and show additional values of interest such as the depths of the 2012 melt layers.

Data acquisition and processing
In preparation for the upcoming East GReenland Ice core Project (EGRIP), the Danish Centre for Ice and Climate's dome and equipment had to be moved about 450 km from the previous drilling site, NEEM.Alongside this so-called "N2E" traverse in May 2015, several measurements of the upper part of the firn and the snow surface were taken.Amongst others, the upper 2 m of the snowpack was sampled using the "liner technique" described in detail below.Snow cores were taken approximately every 25 km at the sites shown in Fig. 1; detailed coordinates can be found in Table 1.

Liner technique
The sampling was done using carbon fiber tubes with sharp edges of 1 m length, 10 cm diameter and 1 mm wall thickness (called "liners").To start off, the first liner was carefully pushed and hammered into the ground until its top was parallel to the snow surface.Nonetheless in a few cases the Figure 1.The N2E traverse route with the measurement sites according to Table 1.
snow core was slightly compacted by up to 2 cm in the vertical direction, visible as a reduction of the snow level inside the tube compared to the surroundings.Subsequently a snow pit of 1 m depth was dug next to the tube and the snow was cut off at its bottom using a metal plate or small saw.The tube was removed and its openings were sealed using matching plastic bags.Then the cutting surface was cleaned and the second liner was inserted right below the first one.Finally the pit had to be deepened to 2 m to once again cut off the snow and take the second liner.Theoretically the process described can be iterated up to an arbitrary depth.However, the area of the snow pit required increases significantly with every meter of depth gained.Sampling the upper 2 m took approximately 2 h per site.

X-ray tomography
The cores were transported to the Alfred Wegener Institute, Bremerhaven, in a frozen condition.All samples were analyzed in the AWI-Ice-CT (described in detail in Freitag et al., 2013), a unique X-ray computer tomograph (CT) in a cold lab, which allows micrometer-resolution density measurements of whole 1 m core segments in 2-D and 3-D.As part of the measurement procedure a sample holder for liners was constructed, which contains several pieces of pure ice of known geometry for calibration purposes.Amongst others, the effect of the carbon fiber tube being part of the scan was corrected for using empty tube measurements.Thus, the fragile snow cores do not have to be removed from the liners.
As the required measurement time increases with resolution, we chose to do 2-D scans with a pixel size of approximately 0.128 mm.Each of these scans takes about 3 min.However, 15 min m −1 are more realistic when including sample preparation and accurate documentation.Then, the raw measurement data are automatically processed by detecting the calibration unit and directly calculating densities from the CT images.Additionally, for each liner, the mean density is determined from the mass and geometry of the snow as an independent comparison value.Figure 2 displays an example CT image with a zoomed section showing two melt layers in the snowpack aligned with the respective densities derived from 2-D analysis.

Isotope measurements
Finally, the snow was gently pushed out of the tubes and cut in samples with a vertical height of 1 cm for the 30 cm right below the surface and 2 cm otherwise.These samples were crushed and sealed in plastic bags.Finally water isotopes were measured using a Picarro L2130-i with a precision of σ = 0.1 ‰ for δ 18 O.
The snow was dated by determining and counting the maxima (summer) and minima (winter) in the seasonal δ 18 O signal.Using the density data, accumulation rates at the different sites were calculated from the snow mass for the 3-5 years worth of accumulation contained in the top 2 m of the snowpack.In the present study, we only use winter-towinter rates (separating years at the δ 18 O minima) -summerto-summer values were computed as a reference but do not show different behavior.
3 Mathematical methods

Automatic alignment of stratigraphic features
In order to efficiently analyze the data sets generated along the traverse, we investigated several ways to automatically detect coherent signals at the different sites.A well-known matching method is maximizing the cross-correlation.However, determining a constant shift in depth between two profiles is not suitable for our case as the accumulation rate, and thus the vertical spacing of layers, is subject to change going eastwards.Under the assumption of constant accumulation over time and no significant compaction in the top 2 m, one would expect a shift which linearly increases with depth and has a slope equal to the ratio of accumulation rates.Then again, local environmental conditions such as wind speed and direction influence the mass accumulated by a certain deposition event (Fisher et al., 1985).Therefore we aimed to align snow of the same origin and its properties with continuously changing shifts, a problem that has already been worked on at a lower vertical resolution for alpine snow (e.g., Hagenmuller and Pilloix, 2016).
The dynamic time warping (DTW) method, which was introduced in speech recognition in the 1970s (Itakura, 1975), provides an efficient algorithm for that purpose.It has already been applied in numerous fields, e.g., for the tracking of ice floes in synthetic aperture radar images (McConnell et al., 1991).For a detailed review of DTW, see Senin (2008).
The basic idea is to discretize the two data sets to be compared with the same step size l (resulting in two vectors S and T of length n and m) and then consecutively assign the values of one to another, whereby each value can be matched with multiple values of the other data set.To find the best fit, one calculates a matrix D, where D[i, j ] indicates the error of the best path that leads to the ith element of the first data set being connected to the j th element of the second one.
The original algorithm starts by calculating the matrix in the upper left corner, fixing the first elements of both data sets to be linked with each other.Then it proceeds through the matrix by taking the path with the minimal error leading to the respective cell and adding the local error, i.e., for i = 0 and j = 0 Finally, on arrival at cell D[n, m] it backtraces the path of minimal errors to D[0, 0], obtaining the best fit of the complete data sets in the given norm • .
For our application -matching measurements of the upper 2 m of the snowpack -we do not aim to fit complete data sets, but rather allow for different offsets at the top and bottom.The former may be caused by variations of the snow surface due to current conditions, the latter by different accumulation rates, leading to data at the bottom of the liners not having any physical relation apart from being the deepest snow analyzed at the given location.To accomplish that, we expand the idea of Sakurai et al. (2007), introducing maximal surface and bottom index offsets s and b.Then we initialize D by www.the-cryosphere.net/10/1991/2016/The Cryosphere, 10, 1991-2002, 2016 and before proceeding through the matrix.Finally instead of backtracing simply from D[n, m], we end our fitting path at and search a trace back to any of the initialized elements.
Thereby we find the best matching of subsets of S and T with a maximal shift of s • l at the top and b • l at the bottom.
In between, we verify that a linearly increasing maximal shift is not exceeded.
The simple way we proceed through the matrix so far, often referred to as "stepping pattern", is unrealistic for our case as a single value of one data set could be fit to arbitrarily many values of the other.Along the traverse we find the maximal ratio of the respective accumulation rates between two sites to be a little smaller than 2. Therefore, we apply constrained stepping as presented by Sakoe and Chiba (1978) such that each value of one data set can be fit to at most two values of the other.This is obtained by Figure 3 illustrates the different patterns for proceeding through the matrix.Here, usage of cell [i, j ] refers to S[i] being assigned to T [j ].In the aftermath, the backtracing has to occur according to the implemented stepping.
Finally, we do not only want to fit one type of data (e.g., densities) but combine the available information in the profiles to gain a robust picture of the developing stratigraphy along the traverse.In a first step, we match the δ 18 O signal, which shows a clear seasonal behavior but almost no small-scale variations, as the high-frequency component is lost by diffusion.Then, we use the obtained depth assignment of the two different sites to resample the measured densities to a common depth scale.In a second step, we apply the algorithm to these densities at a much higher resolution to fine-tune our depth alignment according to small-scale stratigraphic features.As a norm we use the Euclidean distance divided by the path length (i.e., the root mean square error), which means that we have to keep track of the path lengths in a second matrix.Table 2 summarizes the final set of parameters.The maximum allowed offsets for the coarse fitting have been chosen according to the measured height of variations in the snow surface (e.g., dunes) and the maximum ratio of estimated accumulation rates.In the second step we allow for fine-tuning up to the maximum remaining shift, which was manually identified by aligning the vertical centers of the 2012 melt layers.This method does not only allow us to compare data from two sites, but also to obtain a moving depth alignment by fitting the profiles to the first data set one by one.The result, a continuous image of the snow layering, can be compared with other indicators such as the melt layer positions.In addition, being able to align densities and stratigraphic features all along the traverse enables us to provide a representative density profile for the region.For its construction, we first use the continuous layering to transform all density curves to the first depth scale (NEEM) and average them.This, however, is not yet a representative density profile as all profiles now replicate the layering at NEEM; e.g., a layer that is very thin there but thicker at most sites would be considered thin.
To overcome this, we calculate the mean shifts applied to the values that were aligned and thus averaged.On average, i.e., for constant accumulation rates, we would expect these shifts to go linear with depth for the layering to be representative.Thus we calculate a linear least squares regression and correct the depth accordingly.
Nonetheless, the depth scale still represents the accumulation rate at NEEM.To transfer the average profile to any location X in the sampling area of known accumulation (not necessarily one of the N2E sites), we need to calculate a linear rescaling factor f X for the depth d X that fulfills We expect f X to be determined by the accumulation rate, or rather its ratio to the one at NEEM.

Significance testing and surrogate density profiles
Any alignment method will increase the covariance between records even if they are not related (Haam and Huybers, 2010).Therefore, to test the statistical significance of our density alignment, we generate sets of surrogate density profiles with similar statistical properties independently for each site and process them the same way as the original data.
Alongside the artificial density profiles, the real δ 18 O signals are used for the coarse fitting step.
The complexity of the density signal consisting of slow variations, sharp property changes as well as strong melt layer and wind-crust-related density spikes, inhibits the use of simple surrogate construction methods such as autoregressive processes.Instead we propose the following algorithm.
For each site, as a base curve, we identify the δ 18 O component of the density signal by linear regression, using the same step size l low as for the coarse (δ 18 O-based) fitting step.This can be done because we rely on δ 18 O to follow a seasonal cycle -otherwise water isotope dating would be impossible.Let ρ base be the base density from δ 18 O, r low the autocorrelation and σ base the standard deviation of the fluctuations of the measured density (averaged to resolution l low ) around ρ base for lag l low .We start generating an artificial low-resolution density profile ρ low by where Here ν ∼ N (0, σ base ) implies that the ν i are distributed normally with mean zero and standard deviation σ base .In the following, U(0, 1) will represent a continuous uniform distribution for the interval [0, 1].The inclusion of higher autocorrelation lengths is straightforward.r low has to be replaced by the autocorrelation matrix, which is multiplied by a vector of the preceding ε i .Second, on the fine scale (step size l high ), we have a look at the differences between the interpolated low-resolution density and the high-resolution density values from the measurements.As we find the distribution to be trimodal, we split the differences into three components -low amplitude variations within snow of similar properties (henceforth denoted "noise" even though they might partly have physical origin), fast and moderate amplitude changes in the density due to layering or wind crusts ("shocks") and rapid high amplitude changes at melt layers ("melt").Again, we compute the autocorrelation factor r high for lag l high .Nonetheless, this time, the standard deviations σ noise , σ shocks and σ melt and the means µ shocks and µ melt have to be calculated separately.Furthermore we need to estimate the probabilities P shocks and P melt of beginning a shock or a melt layer at a specific position.For this purpose, we determine the number of melt layers N melt , the number of shocks N shocks and the average distance to the previous shock d avg .
In addition, we denote the total number of data points by N and the distance to the last shock at a given position i by d i .
Finally, the basic model to generate a random density profile ρ high is for i = 0 or P > P melt + P shocks N (µ shocks , σ shocks ) for i = 0 and P melt < P ≤ P melt + P shocks N (µ melt , σ melt ) for i = 0 and P ≤ P melt ( 12) where The same approach as before can be used to expand to higher autocorrelation lengths.However, we use the model in the presented form as it already provides realistic density surrogates.

Profile alignment
As an example of the matching process, we present a fit of data from N2E_11 to the first site (NEEM) in Fig. 4. The distance between the two locations is about 240 km, i.e., a little more than half of the total traverse length.First the δ 18 O profiles are matched, yielding an approximately linearly increasing coarse shift.In the second step the densities are finetuned, which results in small shifts fluctuating around zero and never reaching the allowed maximum of 0.1 m.To provide an overview of the changing snow structure, we fitted all combinations of profiles from two sites and plotted the matrix of the root mean square errors (RMSEs) in Fig. 5.A remarkable change in the pattern of the fitting errors occurs between the fourth and fifth site along the traverse.
Figure 6 shows the continuous depth alignment obtained by fitting all liners along the traverse to the first site (NEEM).There were no notable differences when another location (e.g., EGRIP) was chosen as the reference or the fitting was done consecutively.For comparison, the melt layer positions detected during the CT measurements (cf.have been included.In addition, selected density profiles are displayed.Using the previously calculated depth alignment, density records were stacked to obtain a representative density profile (Fig. 7).The gray area indicates a 1 standard deviation error band.Comparing the necessary rescaling factors (known from the construction of the stacked profile) to the ratio of accumulation rates, we apply linear least squares to find where ȧX denotes the mean annual accumulation rate at site X.The coefficient of determination is R 2 = 0.82.3).
At the base resolution of 0.1 cm we find a mean shared variance of R 2 = 0.56 between the average and the individual density profiles.It can be increased by smoothing and obtains a maximum of R 2 = 0.71 when using a 4.3 cm moving average.In comparison, for 1000 randomly generated density data sets (e.g., Fig. 8), the respective stacked profiles share an average of R 2 = 0.44 with their components at base resolution.The maximum is R 2 = 0.61.We determine a p value (probability of finding such high R 2 by chance) of 0.015 for Table 3. Melt layers, the water isotopic season of origin for the surrounding snow and mean annual accumulation rates for each site.The given depths indicate the vertical center of the respective melt layer.The upper two melt layers are always located in snow from summer 2012.For the lower ones, the season of origin for the surrounding snow is given, where S indicates summer and W winter.The accumulation rates are annual mean values for all available years at the particular location.

Site
Depth  the measured profiles within the distribution; i.e., the high shared variance of the measured profiles is statistically significant.

Raw densities, isotope extrema and accumulation rates
For all sites we find at least two melt layers in the snow isotopically dating back to the summer of 2012.In addition, some liners show melt layers which are surrounded by snow dating to winter 2011/2012 or summer 2011.For an overview of all melt layers, see Table 3 or Fig. 6.From the raw density profiles, we obtain Fig. 9, which shows the average densities of the top meter and decimeter, which do not contain any prominent melt layers.The density in the top meter tends to decrease from the maximum of 332 kg m −3 at NEEM down to a minimum of 297 kg m −3 roughly 150 km from EGRIP before slightly increasing again.For 15 out of 18 sites the surface density is higher; nonetheless both parameters evolve similarly along the traverse.Table 3 displays the mean annual accumulation rates along the traverse.Starting with a maximum of 225 kg m −2 a −1 at NEEM, the values steadily decrease down to the minimum of 115 kg m −2 a −1 about 100 km from EGRIP before increasing again to 140 kg m −2 a −1 at EGRIP.Comparing average values for the different years, there is neither a trend nor considerable variations in the accumulation rate (cf.Table 4).However, we observe much higher differences between successive years within the same core (average change 34.67 kg m −2 a −1 ), where we mainly see alternating behavior of high-and low-accumulation years.
Of the 5 years contained in our data, 2012 had the isotopically warmest summer for 83 % of the sites.At the three remaining locations (N2E_11, N2E_16 and EGRIP), the highest δ 18 O values occur in 2014.For the winters, 2014/2015 was isotopically coldest in 51 % of the cases, 2011/2012 in 19 % and 2010/2011 in 30 %. Regarding annual δ 18 O averages of all available sites (Table 4), we also find the highest δ 18 O values for 2012.

Linking accumulation, δ 18 O and density
Comparing the annual average δ 18 O values with the accumulation rates we obtain Fig. 10.Positive linear relations were fit to the data of 2012, 2013 and 2014 respectively, showing that within 1 year higher temperatures coincide with higher accumulation.The coefficient of determination is highest for 2012, while we have larger spreads for the other 2 years, in particular 2013.
To relate the density with the seasonal, low-frequency δ 18 O signal at NEEM, we applied a 10 cm running mean to the stacked high-resolution density profile in Fig. 11.On average, snow with a high δ 18 O value (considered summer snow) has a low density and the other way around.The only exception is the summer of 2012, where we find high-density values in summer, too.

New methodology
The liner technique allows us to retrieve non-disturbed snow samples from the field and thereby conduct lab-based analysis (such as high-resolution density measurements) to gain  further insight in the development of physical snow properties over large distances.This is a major improvement compared to previous methods, e.g., for measuring snow density, which was so far mainly done by weighing a known volume of snow where we have a trade-off of accuracy (bulk density) and resolution (density cutters).Both horizontal resolution and vertical depth can be adjusted to fit the needs of the respective study.
Figure 4 illustrates that we are able to align δ 18 O and density data down to small stratigraphic features very well over a distance of over 200 km.Along the traverse, one observes a clear change in the RMSE (cf.Fig. 5) and thereby the snow structure at the fourth site, indicated by significantly different fitting errors.This coincides with the location where the ice divide was left eastwards and thereby the traverse entered a different accumulation regime in agreement with the drainage systems given by Zwally et al. (2011).
Furthermore the continuous depth alignment agrees very well with the melt layer positions detected during the CT measurements (Fig. 6).Stratigraphic features are still well aligned over the complete traverse distance of almost 450 km.We obtain a clear picture of the layering of the snowpack along the traverse.In comparison to radar measurements, which are limited to centimeter vertical resolution but can resolve annual layers down to 12 m (Hawley et al., 2006), we can give a much more precise picture and observe smallscale structures like wind crusts.In exchange we are limited to shallower depths -the maximum we plan to access in the near future is 6 m in a trench at the EGRIP drilling site.
For rescaling the stacked profile to any location in the area with known annual accumulation, we obtain a linear relation of the depth factor with the ratio of accumulation rates.This is plausible because, on average, we find linearly increasing shifts for the matching.Furthermore we do not observe significant densification in the upper 2 m of the snowpack, and therefore the depth of snow from the same deposition event is primarily determined by the accumulation rate.In addition, the relation has a high coefficient of determination for the applied linear least squares.
As the stratigraphy does not seem to change remarkably along the traverse apart from the effect of the decreasing accumulation rate, we consider the profile in Fig. 7 to be representative of the whole traverse region, potentially even most of North Greenland.For the given error band, there is an overlap of uncertainty in the depth alignment (x direction) with the uncertainty in density (y direction).The former is mainly caused by the variability of the snow mass accumulated from a single deposition event.Regarding the latter, the average density of the snowpack greatly varies as can be seen in Fig. 9. Thus, for the second meter, even though it is contained in the uncertainty band, we do not expect a straight line, but rather an alternation of high-and low-density layers similar to the upper meter.
A statistical test using surrogate density profiles shows that the high shared variance of the measured profiles is statistically significant (p = 0.015), even though the actual difference in numbers is quite small.This underlines that the density alignment provides additional information as we tried to use the most realistic surrogates (original δ 18 O signal, seasonal cycle, three component stratigraphy model).Furthermore, a coefficient of determination of R 2 = 0.56 between the stacked and the individual profiles shows how much of the layering does reappear.Smoothing increases R 2 as it steadily transforms the profile to the low-resolution density curve that shows seasonal behavior (see Fig. 11) while smaller local variations vanish.

Temporal and regional variability of snow properties
The vast majority of melt layers are found in snow dating back to the very warm summer of 2012 (Nghiem et al., 2012).Moreover, above most of the melt layers within older snow, we find clear signs of percolation (cf.Fig. 2).Therefore we assume that 2012 was the only year in the period 2010-2015 with significant melt occurring in the observed area.From Fig. 9 we can infer that on the one hand, the average density of the snow in the top 2 m at a certain location can already be deduced from the surface density.On the other hand, the surface snow in May is among the denser ones within the year, thereby rather representing a spring or even winter signal than a summer one (compare Fig. 11).Furthermore we are able to visually identify many layers of homogeneous density, often clearly separated by wind crusts, which seem to contain snow from single deposition events.For the accumulation rate (see Table 3), the 1964-2005 average of 220 kg m −2 a −1 determined from the NEEM ice core (Steen-Larsen et al., 2011) agrees very well with the 225 kg m −2 a −1 that we obtain from the corresponding snow liner.In addition, both accumulation maps from field measurements (Bales et al., 2009) and regional climate models (Fettweis, 2007) show the same behavior towards the east.While Table 4 shows no significant interannual changes in the average accumulation rate for the study area, we observe www.the-cryosphere.net/10/1991/2016/ The Cryosphere, 10, 1991Cryosphere, 10, -2002Cryosphere, 10, , 2016 high fluctuations in the local annual values, a feature consistent with the strong influence of stratigraphic noise in single profiles (Münch et al., 2016).These can be explained by the accumulation of every year compensating previous local variations in the snow surface before new structures are introduced by wind-induced drift and dunes.Nonetheless, they also might partly originate from the uncertainty of separating the years only according to the δ 18 O extrema.
In the majority of cases we find the highest isotopic summer temperatures and average δ 18 O values for 2012, underlining the exceptional warmth of this year.The values for 2014 indicate that it was still warmer than the other contained years, in particular 2010, which was formerly regarded as very warm (Harper et al., 2012).The picture for the winters is less clear.Indeed, we assume that the isotopic signal of the fresh snow from winter 2014/2015 might still change.

Relations of density, δ 18 O and accumulation rate
We find a positive linear relationship of annual mean δ 18 O and accumulation rate (Fig. 10) with similar slopes for 2012 and 2014.This relation might partly originate from the changing surrounding conditions (e.g., elevation) along the traverse.The offset between the years could potentially be caused by the very high temperatures and the consequential surface melting in 2012 as we find the relation for 2013 to be a lot closer to 2014 than 2012.The dependence of the offset on the annual mean temperature (which is quite similar along the traverse) could explain why previous attempts to link both parameters by averaging data from several years (e.g., Weißbach et al., 2016) show results that are less clear.
We observe a clear anticorrelation of low-resolution density and δ 18 O in Fig. 11.This agrees with the widely accepted conceptual model of Shimizu (1964), which states that snow has lower densities in summer and higher ones in winter.The main causes given are the increased packing due to stronger winds in winter and the larger size of precipitation particles in summer.For the summer of 2012, the high average densities are caused by the prominent melt layers, superimposing the original signal of the snow.

Summary and conclusions
We introduced the liner technique, which allows the very efficient retrieval of high-quality samples from the upper meters of the snowpack.To support this new sampling technique, we adapted a robust fitting algorithm from acoustic signal processing for the diverse data sets produced by such studies.This enables us to identify characteristic changes in the snowpack according to surrounding conditions as well as to generate continuous depth alignment using features from all available records.
To demonstrate their feasibility we applied the described methods to the upper 2 m of snow along a traverse in North Greenland.We obtain a record up to May 2015 of the depths of the 2012 melt layers and submillimeter-resolution densities.By combining these with δ 18 O measurements, which indicate temperature, we are able to reconstruct accurate accumulation rates for the years 2010-2014 along a distance of about 400 km.
We combine isotope and density data as inputs for the matching algorithm.Thereby we are able to identify the different accumulation regimes along the traverse and resolve the continuous stratigraphy of the snow over the whole distance.This allows us to create a representative density profile for the study area, whose quality is proven by comparison with randomly generated data based on a statistical density model.The profile is available at a resolution of 0.1 cm and only has to be rescaled according to accumulation rate.Thus it is ready to act as a benchmark for snowpack models or be applied for the conversion of volume to mass and the detection of strong density gradients as potential reflectors in remote sensing (compare e.g., Hurkmans et al., 2014).
The success of fitting density and isotope profiles over hundreds of kilometers shows that even though there is a local component in the snow stratigraphy (e.g., layer thickness, average density), the general pattern is dominated by non-local processes in North Greenland.We assume that an important factor for that is the origin of weather and precipitation as air masses dominantly move in from the west to the east (Chen et al., 1997).
We observe large interannual accumulation variations locally but almost none on average, which can be explained by the smoothing of the surface by accumulation before new surface structures are caused by dunes and drift.The exceptionally warm summer of 2012 is clearly visible in the water isotope data; additionally 2014 shows the second highest summer values of δ 18 O within the study period.
Relating the various snow properties, we find a distinct anticorrelation of smoothened density and δ 18 O in accordance with previous literature.Furthermore we deduce a positive linear relation between δ 18 O and accumulation rate, whose slope seems to be constant for the period considered, while the offset varies between the years, and thus might be temperature-dependent.This, however, poses the question as to whether models commonly used in the dating of deep ice cores (e.g., Parrenin et al., 2007, for the EPICA Dome C ice core) do correctly reconstruct accumulation rates from the δ 18 O values, especially for times with significantly differing annual mean temperatures such as glacials.
Future work should include the automatic recognition of wind crusts and layering from CT images and the application of the described methods on different scales for both Antarctica and Greenland to gain further insight into the variability of physical properties in the snowpack.

Figure 2 .
Figure 2. Example 2-D CT image of a 1 m liner (1-2 m depth) and a zoomed section showing two melt layers aligned with the respective densities.In the left image a distinct density layering (e.g., blue triangle), several melt layers (e.g., blue circle) and wind crusts (e.g., blue square) are visible.Above the lower zoomed melt layer a clear percolation pattern (blue arrow) can be seen on the right-hand side of the snow core.

Figure 3 .
Figure 3. (a) Basic and (b) constrained stepping patterns for the DTW algorithm.Usage of cell [i, j ] indicates that the ith element of the first and the j th element of the second data set were matched.The basic pattern allows for a single value to be assigned to arbitrarily many values of the other data set, while for the constrained stepping, each value can only be identified with one or two others.

Figure 4 .Figure 5 .
Figure 4. Alignment of the data from NEEM and N2E_11.(a) First, the raw δ 18 O data from N2E_11 (orange) are fit to those of NEEM (blue) resulting in the red curve.(b) Then, the calculated (coarse) shifts are applied to the raw N2E_11 density data to obtain the red curve as an input for a second alignment with the raw NEEM density profile (blue).We end up with the pink curve as a final result.(c) The applied coarse (black) and fine (gray) shifts.

Figure 6 .
Figure6.Continuous depth alignment, example density profiles and melt layers.A color map was applied uniformly at the first site (NEEM) and then transformed the same way as the depths were assigned.Thus snow within the same color band was matched during the fitting process.Linear interpolation was used between the sampled sites.In black, measured density profiles for the labeled locations are shown at the same scale, centered around their respective mean values.The white lines and points indicate the melt layer positions detected from the CT scans (cf.Table3).

Figure 7 .
Figure7.Representative density profile for the traverse region.The gray area indicates a 1 standard deviation error band in both x and y directions as there are uncertainties in the depth alignment as well as the averaged densities of all sites.Here, the depth scale was adjusted to the NEEM accumulation rate and has to be rescaled according to the accumulation rate for different sites.

Figure 8 .
Figure 8.The measured density profile and three surrogates for the first site (NEEM).The artificial profiles are based on the seasonal δ 18 O component of the density and have the same statistical properties as the original curve.Each profile is displayed at the same scale and has been centered around its mean.

Figure 9 .
Figure 9. Average densities along the traverse through North Greenland (May 2015) in the top 1 and 0.1 m derived from CT data.

Figure 10 .Figure 11 .
Figure 10.δ 18 O signal vs. accumulation rate for the years 2012-2014.The lines were obtained by linear least squares fitting with coefficients of determination of R 2 = 0.52 for 2012, R 2 = 0.27 for 2013 and R 2 = 0.37 for 2014.The data points for 2013 show the largest spread and were omitted for clarity.

Table 2 .
Fitting parameters for our adaption of the DTW algorithm.

Table 4 .
Mean deviations of the given year from the average local annual (winter-to-winter) accumulation rate and δ 18 O.For each year, data from all available sites were used.