Unlocking annual firn layer water equivalents from ground-penetrating radar data on an Alpine glacier

The spatial representation of accumulation measurements is a major limitation for current glacier mass balance monitoring approaches. Here, we present a method for estimating annual accumulation rates on a temperate Alpine glacier based on the interpretation of internal reflection horizons (IRHs) in helicopter-borne ground-penetrating radar (GPR) data. For each individual GPR measurement, the signal travel time is combined with a simple model for firn densification and refreezing of meltwater. The model is calibrated at locations where GPR repeat measurements are available in two subsequent years and the densification can be tracked over time. Two 10.5 m long firn cores provide a reference for the density and chronology of firn layers. Thereby, IRHs correspond to density maxima, but not exclusively to former summer glacier surfaces. Along GPR profile sections from across the accumulation area we obtain the water equivalent (w.e.) of several annual firn layers. Because deeper IRHs could be tracked over shorter distances, the total length of analysed profile sections varies from 7.3 km for the uppermost accumulation layer (2011) to 0.1 km for the deepest (i.e. oldest) layer (2006). According to model results, refreezing accounts for 10 % of the density increase over time and depth, and for 2 % of the water equivalent. The strongest limitation to our method is the dependence on layer chronology assumptions. We show that GPR can be used not only to complement existing mass balance monitoring programmes on temperate glaciers but also to retrospectively extend newly initiated time series.


Introduction
Mountain glaciers are not only known to be excellent indicators of climatic change (Haeberli and Beniston, 1998;Kaser et al., 2006;IPCC, 2013) but also important to the hydrology and ecology of Alpine regions (Beniston, 2003).In order to better understand their response and adaption, glacier changes are monitored around the globe (Barry, 2006;Zemp et al., 2009).Conventional mass balance measurements are prone to large uncertainties when compared to geodetically derived volume change.The reconciliation of the two approaches is the subject of ongoing research (Huss et al., 2009;Fischer, 2011;Zemp et al., 2013).Ablation can be measured accurately at many locations across a glacier.For many study sites, the main drawback is a lack of accumulation measurements that typically involve the time-consuming excavation of snow pits (Østrem and Brugman, 1991).If measurements of accumulation are under-represented, the spatial variability is often not resolved correctly in modern monitoring approaches (e.g., Huss et al., 2014).To compensate for this drawback, an adequate complement to conventional accumulation measurements is necessary, which should preferably be non-destructive and efficient.
In glaciology, ground-penetrating radar (GPR) is used for a variety of purposes (Plewes and Hubbard, 2001).The low electrical conductivity of snow and ice facilitates the deep penetration of the signal, and the use of narrow-bandwidth, high-resolution systems (Ulriksen, 1982).Because the signal is reflected from boundaries with a contrast in dielectric permittivity, not only a change in the related material properties can be detected, such as from ice to bedrock, water content, or impurities, but also changes in density are resolved (Plewes and Hubbard, 2001;Woodward and Burke, 2007).Thus, GPR applications range from ice thickness measurements (e.g.Robin et al., 1969;Bauder et al., 2003), to the mapping of snow accumulation distribution (Machguth et al., 2006;Sold et al., 2013;Helfricht et al., 2014), to the observation of internal layers within the snow cover (e.g., Heilig et al., 2010).Thereby, airborne GPR surveys can cover large areas in short periods.On ice sheets with rather homogeneous accumulation patterns, internal reflection horizons (IRHs) within firn and ice can be tracked over long-ranging GPR profiles (Arcone et al., 2004;Huybrechts et al., 2009;Miège et al., 2013;Hawley et al., 2014).Past accumulation rates are typically estimated from a set of pronounced IRHs that do not correspond to annual layers.Based on a given depth-age relation and the density, an average mass balance is obtained for the period that is covered by each pair of IRHs.Only a few studies have observed and analysed the water equivalent (w.e.) of annually spaced IRHs in radar data (e.g., Kohler et al., 1997;Hawley et al., 2006;Kruetzmann et al., 2011;Van Pelt et al., 2014), but not on temperate mountain glaciers in mid-latitudes to low latitudes.
On temperate glaciers, or glaciers that are constrained by topography and have complex flow fields, the continuous tracking of IRHs is often difficult (Karlsson et al., 2009).Therefore, studies have focused on the mapping of facies zones and other general glaciological characteristics (Pälli et al., 2003;Langley et al., 2008;Eisen et al., 2009;Dunse et al., 2009).On the other hand, GPR not only could accompany ongoing monitoring programmes to complement the sparse accumulation measurements with high spatial resolution but also could be used to retrospectively extend newly initiated time series into the past.
To extract past annual accumulation rates from the GPR signal, the IRHs must correspond to previous summer surfaces.In polar and sub-polar regions this is confirmed by several studies (e.g., Wadham et al., 2006;Van Pelt et al., 2014).For mid-latitude glaciers with a complex firn stratigraphy it can be difficult to establish this link (Kohler et al., 1997), although the large number of melt-refreeze cycles suggests the generation of a high-density or ice layer at the snow surface during summer.Thus, this precondition must be verified by independent layer dating information.
Furthermore, an estimate of the propagation velocity is required to convert the GPR travel time to depth.Because the velocity depends on the density of the material, the latter needs to be measured or estimated (Plewes and Hubbard, 2001).Thus, density introduces an uncertainty, as it also does for conventional accumulation measurements that are based on the determination of layer thickness.In order to obtain the density and water equivalent of firn layers, meltwater that percolates from the surface into the firn must be taken into account.Aside from its effect on the propagation velocity (Schmid et al., 2014), it changes the density and water equivalent of firn layers when it refreezes.
Studies that analyse IRHs in firn along GPR profiles are typically complemented by cores to provide density, layer dating, and potentially, dielectric characteristics (Pälli et al., 2002;Arcone et al., 2004;Eisen et al., 2006).However, firn cores are expensive in terms of cost and effort and provide point measurements only.Physical models can be used to estimate firn density profiles along GPR profiles.They require various input data such as temperature, precipitation, wind, and terrain to calculate the accumulation rate (Ligtenberg et al., 2011).
In this study we obtain IRH travel times from GPR data measured in May 2012 by 400 MHz helicopter-borne GPR on Findelengletscher, a temperate Alpine valley glacier in Switzerland.The analysis is limited, however, to the upper part of the glacier where sufficient accumulation allows for distinguishing multiple IRHs.Due to IRH discontinuities and variations in the GPR signal penetration depth, the analysed subset consists of disjoint profile sections.The chronology of the IRHs is directly retrieved from two shallow firn cores.Comparison with density and impurity profiles shows that not all IRHs represent former summer surfaces.The measured travel times of annual IRHs serve as input to a simple transient model for the density of firn layers and the refreezing of meltwater.The model is calibrated with an observed change in the IRH travel times at multiple locations.They are obtained from a repeat survey conducted in April 2013.The model is run for all analysed GPR measurement locations separately and provides estimates for the density, water equivalent, and meltwater refreezing in annual firn layers.Thus, on the basis of the layer chronology from firn cores, we estimate annual accumulation rates for Findelengletscher over a multiyear period from airborne GPR measurements.

Study site and field data
Findelengletscher is a 13 km 2 temperate Alpine valley glacier in the vicinity of Zermatt, Switzerland, ranging from 2600 to 3900 m a.s.l.It is situated directly below the main Alpine ridge with a north-westerly exposure (Fig. 1).A monitoring programme for annual mass balance was started in 2004 (Machguth, 2008) with a network of 13 ablation stakes and two snow pits for density measurements in the accumulation area (Fig. 1).Since 2009, the data record is complemented by measurements of winter mass balance using conventional snow probings and snow pits.The seasonal mass balance measurements are extrapolated to the entire glacier by constraining a distributed accumulation and temperatureindex melt model operating on a daily resolution (Huss et al., 2009).The model is calibrated for each year individually using the point measurements of ablation and accumulation and takes into account spatial variations of potential solar radiation and snow accumulation.Meteorological data are available from automatic weather stations at Zermatt (1638 m a.s.l.; 6 km west of the glacier terminus) and Stockhorn (3421 m a.s.l.; 2.5 km south of the glacier terminus).A comprehensive description of the model is given in Huss et al. (2009).Since the onset of measurements, the average annual mass balance of Findelengletscher has been negative which has also been accompanied by a retreat of the glacier snout (Glaciological Reports, 2011).
Helicopter-borne constant-offset 400 MHz GPR measurements were taken at Findelengletscher on 10 May 2012 and 14 April 2013.Thus, the uppermost firn layer is covered by the recent winter accumulation layer.The surveys were conducted along regular 500 m × 500 m grid lines covering the entire glacier and had a total profile length of 79 km.With a flying speed of approx.10 m s −1 , measurements were taken from 5 to 20 m above the surface at time intervals of 0.02 s, corresponding to a trace spacing of approximately 0.2 m.Together with the position obtained from a differential global positioning system (DGPS), a time window of 250 ns was recorded for each trace at a waveform sampling interval of 0.24 ns.
To obtain ground reference data, two shallow firn cores were recovered by mechanical drilling in the upper accumulation area of Findelengletscher on the same date as the GPR measurements in May 2012 (Fig. 1).The choice for the two sampling sites on a smooth ridge at 3495 and 3530 m a.s.l. was based on the availability of GPR profiles of preceding years at these locations (Sold et al., 2013) and an expected slow glacier flow.The recovered cores were about 10.5 m long and 58 mm in diameter.

Firn sampling, analysis, and dating
The two 10.5 m firn cores recovered during the campaign in spring 2012 using the electromechanical drill FELICS (Fast, Electromechanical Lightweight Ice Coring System) small (Ginot et al., 2002) were stored in polyethylene tubes in sections of 0.5-0.7 m length at −20 • C.After a visual inspection for melt features or dust layers, the outer part of the core was removed and the sections were cut with a resolution of approx.0.06 m using a bandsaw with stainlesssteel blades, Teflon-covered tabletops, and saw guides to derive decontaminated samples with dimensions of approx.60 mm × 17 mm × 17 mm They were melted in individual pre-cleaned 50 mL plastic tubes under inert gas (N 2 ) before the analysis.The density of each sample was determined gravimetrically from its measured volume.The concentration of the major soluble inorganic cations and anions including ammonium (NH + 4 ) and sulfate (SO 2− 4 ) was determined by ion chromatography (Metrohm 850 Professional IC combined with an 872 Extension Module and an 858 Professional Sample Processor autosampler).Additionally, the water stable isotope ratio (δD) was measured using Wavelength-Scanned Cavity Ring Down Spectrometry (Picarro L2130-i).
Dating of the two firn cores was performed by annual layer counting, using the pronounced seasonality of ammonium concentrations and water stable isotope ratios δD, both peaking in summer (Dansgaard, 1964;Eichler et al., 2000) (Fig. 4).A distinct dust layer from the Saharan dust fall on 28 May 2008 (Hostettler and Bader, 2008) found at 10.2 m (core 1) and 8.5 m (core 2) depth was used to corroborate the dating.Core 1 and 2 cover the time periods of winter 2007/08 to spring 2012 and summer 2007 to spring 2012, respectively.

GPR data processing
A sequence of processing steps was applied to the GPR data in order to reduce noise and to improve the general visibility of reflectors.The selection of individual processing steps and their parameter settings depends on the field conditions and survey intention; thus, no universal procedure is available (Annan, 1993;Fisher et al., 1996;Ulriksen, 1982).However, we applied a scheme that was previously used when mapping the snow accumulation distribution on Findelen- gletscher (Sold et al., 2013) but with some adjustments to the parameter sets.The scheme consists of a spatial interpolation of traces to a 0.3 m spacing to correct for variations in the helicopter's velocity, a frequency bandpass filter and background removal to reduce noise, and a gain function to account for signal attenuation with depth.Additionally, the profiles were flattened to a horizontal snow surface that is represented by the first and most pronounced reflector.Migration was found to not improve the data quality for the given survey design, due to the moderate surface slopes at the analysed GPR measurement locations (95 % less than 15 • ).The processing was performed with the software Reflexw (Sandmeier Scientific Software) and sufficiently improved data quality without the risk of over-processing and generation of artificial reflectors (Fig. 2).The analysis of firn layers is limited to locations where accumulation occurs.To extract an annual firn layer, at least two IRHs must be present, i.e. representing two subsequent summer surfaces.Furthermore, it must be ensured that between the upper-and lowermost detected IRHs all other IRHs of summer surfaces are present and can be distinguished.Thus, areas with little or no accumulation, as well as profile sections where the interpretation was ambiguous, were excluded from the data set.The reduced set of processed GPR data consists of disjunct profile sections in the accumulation area of Findelengletscher (7.3 km total length for the uppermost firn layer in May 2012) (Fig. 1).The penetration depth of the GPR signal depends on the height of the helicopter above ground and the physical properties of the subsurface.Therefore, the IRHs of deeper, i.e. older, layers appear in shorter sections, leading to differences in the spatial extent of extracted layers.We found reflections down to about 20 m below the snow surface and a maximum of eight IRHs.They were digitised manually in the processing software due to noise and lateral variability in layer thickness which impede a reliable automatic tracking of IRHs.For the evaluation of the GPR repeat measurements in May 2012 and April 2013, we aimed at locating the correspondent of each reflector at its new depth in the second year.A matching procedure was applied on the travel-time domain at 12 locations where a sequence of IRHs was found in both GPR data sets.(Fig. 1).Thereby, the uppermost IRH in the second year was omitted, representing the added first year's accumulation.A free scaling of ±20 % of the depth of the uppermost reflector in the first year was applied, which experiences the strongest compaction over the year.Then, the one-to-one mapping of reflectors was found with the minimum distance in the time domain.

Modelling firn density
An estimate for the dielectric permittivity of the firn is required to obtain the radar-wave velocity, converting the twoway travel time of the GPR signal from time to depth domain.Density strongly affects the dielectric properties of snow and firn and several empirical formulas exist to estimate the permittivity (e.g.Kovacs et al., 1993;Frolov and Macheret, 1999).In addition, the density is required to ultimately determine the water equivalent of firn layers.
We used the approach by Reeh (2008) to model firn layer compaction.It is based on the empirical steady-state model of Herron and Langway (1980), which assumes density changes over time t to be linearly related to the pressure of overlying snow or firn: where ρ f is the density of the firn fraction of an annual layer, ρ ice is the density of ice, and c is a dimensionless parameter depending on the mean annual firn temperature T (K) and the average mass balance b (m w.e.yr −1 ) that controls the overlying weight.It is separated by ρ f = 550 kg m −3 to respect different stages during the densification process (Herron and Langway, 1980).Here, we neglect the lower stage (ρ f < 550 kg m −3 ) because the model will be calibrated and measured autumn snow densities were not considerably lower (492 kg m −3 , see below).For ρ f ≥ 550 kg m −3 c is obtained as where R = 8.314 J K −1 mol −1 is the gas constant, E = 21 400 J mol −1 is the activation energy, and f (m −0.5 yr −0.5 ) is an empirically derived constant controlling the densification rates.This steady-state approach was used by Reeh (2008) but with variable accumulation rates, initial densities, and air temperatures.The density ρ f (t, t d ) of an annual firn layer t years after its deposition at t d is then given by solving Eq. ( 1) as an initial value problem with ρ f (0, t d ) being the layer's individual initial density.Thereby, instead of assuming steady-state conditions, an individual c is obtained for each layer from the mean annual mass balances b prev and mean annual temperatures of the t − 1 preceding years since deposition.
In the context of our study, the local annual mass balances are initially unknown and the model cannot be applied as it is.However, along the GPR profiles, the IRH travel times of firn layers and the snow cover are given.We used this information to progressively derive the unknown densities ρ f,i and water equivalents ω i of all i = 1, . .., n layers at the time of the GPR measurements, i.e. t i years after each layer was deposited at t d .We follow Reeh (2008) by substituting the long-term mean mass balance in Eq. ( 2) with the mean water equivalent of the i − 1 overlying firn layers to obtain an individual proportionality factor c i for each firn layer as This allows for estimating the density ρ f,i of a firn layer by solving Eq. (1) over t i , the time since layer deposition, if the initial density ρ f (0, t d ) and the water equivalents of overlying firn layers are given.Then, the difference in IRH twoway travel time τ i can be converted to the firn layer water equivalent ω i by with the radio wave velocity v(ρ f,i ) estimated from the density using the empirical relation by Frolov and Macheret (1999).
When the density of the uppermost layer is known, Eq. ( 5) provides its water equivalent and Eqs. ( 4) and ( 1) can be used to estimate the density of the second layer.The densities and water equivalents of all other layers can then be calculated similarly in a stepwise manner by using the previously obtained water equivalents of all overlying layers.In our case, the uppermost layer (i = 0) is the winter snow cover.We set its density to ρ f,0 = 400 kg m −3 , which is in line with observations in spring 2010-2013 (395 ± 32 kg m −3 from 14 measurements ≥ 3200 ma.s.l.).According to Reeh (2008), the deposition of a layer at time t d is attributed to the end of the respective melting season.Because the GPR measurements were taken in spring, we set the age of the first annual firn layer (i = 1) as t 1 = 0.3 yr.This also accounts for the associated compaction processes being stronger during the summer months due to higher temperatures and the refreezing of meltwater.For all subsequent layers we set t i+1 = t i + 1 yr.
To account for the temperate nature of Findelengletscher, temperature was set to T = 0 • C = 273.15K (constant).For the initial firn density we used ρ f (0, t d ) = 492 kg m −3 as the mean density measured within the conventional glaciological accumulation measurements at the end of the summer in 2010-2013.We derived a conservative uncertainty estimate of ±61 kg m −3 from the mean standard deviations of (1) multiple density measurements within the same snow pits, i.e. the measurement error (±21 kg m −3 ); (2) within single years, i.e. the spatial variability (±31 kg m −3 ); and (3) at locations with annual repeat measurements, i.e. the temporal variability (±49 kg m −3 ).The initial density is assumed to be spatially and temporarily constant.However, as the uppermost firn layer is of age t 1 = 0.3 yr, it experienced compaction.Thus, this assumption does not imply a spatially homogeneous layer density, but rather defines the initial (i.e. the preceding autumn) model condition.

Refreezing of meltwater
Due to cooling from the surface during winter, firn can periodically have temperatures below 0 • C even on temperate glaciers.The relocation and refreezing of surface meltwater alters the water equivalent of firn layers and affects their density.Because of the temperate conditions, we assume that during each summer all subfreezing temperatures are compensated by the refreezing of meltwater.This is supported by the direct observations of summer mass balance that generally indicate > 1 m w.e. of snow melt during the summer months.Thus, the amount of refreezing is given by the firn temperature profile that is generated by the surface cooling during winter.We calculate a one-dimensional end-of-winter firn temperature profile at each GPR measurement location using heat conduction from the surface in the period from 1 October to 1 May.A vertical density profile with a resolution of 0.1 m was set up, describing the winter snow cover and the underlying firn layers.The winter snow cover is built up linearly between 1 October and 1 May to reach the thickwww.the-cryosphere.net/9/1075/2015/The Cryosphere, 9, 1075-1087, 2015 ness estimated from the IRH travel time.The density was assumed to increase linearly over time from 100 to 400 kg m −3 .For the firn section the respective modelled end-of-winter density profile was used.The thermal diffusivity of snow and firn was estimated from an empirical relationship with density (Calonne et al., 2011).The heat equation could then be solved numerically on a 1 h resolution with a forward in time, central in space scheme.
We applied the model to every annual layer and all winter periods since its deposition using the daily mean air temperature at the automatic weather station at Stockhorn, corrected for elevation with a moist adiabatic lapse rate of −0.006K m −1 .In order to account for differences between the air and snow surface temperatures we applied a mean temperature offset for each month.The offset was based on the observed monthly mean difference in 2002-2010 at the Stockhorn weather station that was estimated from the outgoing long-wave radiation.The offset is the largest in December (−6.7 K) and decreases to −3.5 K in April.For each year and each layer, the amount of refrozen meltwater was calculated from the specific heat capacity of snow and the latent heat of fusion and was added to the layer water equivalent.The given layer thickness then allowed for updating of its modelled density.
However, the two models for temperature and density were not fully coupled because they run on different domains.The densification model primarily steps through the layers from top to bottom, whereas the heat conduction is solved over time.Thus, each year's density profile remains static in the temperature model, whereas the layer water equivalents stay constant within the densification model but are adjusted separately.This is in line with the model itself neglecting variations in the loading history of each firn layer.

Model calibration
Aside from the refreezing, the change in density over time is controlled by the density contrast to the density of ice and the weight of overlying layers (Eq.4).The model sensitivity to those variables is controlled by the factor f (Eq.3).Changes in f propagate through the entire model as they affect not only the densities of individual layers but also their water equivalent, their weight, and, hence, the modelled densities of the following layers.Because the initial density ρ f (0, t d ) can be constrained with conventional measurements in autumn, we used the parameter f to calibrate the model to our site characteristics.
At locations where repeat GPR measurements are available in 2012 and 2013, the density change of each layer was modelled over the time of 1 year.To reduce the effects of ice flow, a potential imprecision of the positioning and differences in the GPR footprint size, the model was applied using the mean IRH travel times within a distance of 25 m from the common location.In order to obtain the new overlying weight on each layer, the IRH travel times of the snow cover and the new annual accumulation layer were considered.The resulting layer densities for the second year were used to estimate the respective hypothetical GPR signal travel times.The model was calibrated with a scaling factor of f = 2614 m −0.5 yr −0.5 (Eq.2) that was found by minimising the root mean square deviation of the modelled and measured IRH travel times of 35 layers at 12 locations (Fig. 1).However, the temporal density change also depends on the density contrast that is linked to the initial density.Therefore, the optimal value of f is related to the choice of ρ f (0, t d ) (Fig. 3).

Local results at firn core locations
We observed a variable density profile for both firn cores that increased with depth (Fig. 4).Maximum densities corresponded to melt features of bubble-free or bubble-poor ice.They are more transparent and form when surface snow melts and the meltwater percolates into deeper layers where it fills the pores and refreezes.In many cases melt features involve enhanced concentrations of NH + 4 and δD.Due to the temperate conditions at the drilling site, refreezing of meltwater not only contributes to the densification of firn and alters the water equivalents of firn layers but also strongly influences the chemical and stable isotope composition (Ginot et al., 2010).The drilling took place in spring 2012 and thus, the snow accumulation of winter 2011/12 did not experience major melting.In contrast, the firn layers below experienced one or more cycles of melting and refreezing, leading to an enhanced smoothing of seasonal variations and potentially, a partial relocation of chemical species.Nevertheless, in agreement with studies on the nearby polythermal Grenzgletscher (Eichler et al., 2001), NH + 4 and δD seasonal cycles were still detectable, even under temperate conditions.
The measured density profile allowed for a direct estimation of the GPR signal velocity using the empirical relation by Frolov and Macheret (1999).On the travel-time domain, we observed agreement between the GPR signal and pronounced changes in the density profile (Fig. 4), suggesting that the reflectors were composed of dense layers.This underlines the high potential of GPR for the detection of summer layers.However, not all density maxima and melt layers correspond to a former glacier summer surface, but can also originate from exceptional melting and refreezing in other seasons, or from percolation of meltwater into deeper parts of the firn.In core 1 at a depth of 9.2 m (7.2 m in core 2), the accumulation layer of 2009 contained a pronounced ice layer of 0.12 m (0.15 m) thickness, corresponding to maximum measured densities of 856 kg m −3 (853 kg m −3 ) (Fig. 4).The low δD ratio implies low temperatures during the formation of the precipitation that was incorporated into the ice layer by in situ densification or refreezing.We suppose that this layer was formed by surface meltwater that was hindered from further percolation by a blocking layer, such as a boundary with a change in the grain structure or a thin ice layer (Pfeffer and Humphrey, 1996).This is supported by a high concentration of most of the measured ions right below the ice layer (e.g.SO 2− 4 ) that typically just extends over one or two samples.By preventing meltwater from percolating and refreezing underneath the ice layer, the chemical composition is preserved and does not experience enhanced smoothing.After a period of major snow falls, high wind speeds and mild temperatures were registered on 20 December 2008.It was followed by clear weather until mid-January (Etter et al., 2011).We suggest that this allowed for the formation of a dense layer at the surface of a fresh snow pack, creating a hydraulic barrier for meltwater formed later in the season.Because no defined clear ice layer was built up in the following summer of 2009 (Fig. 4), meltwater of several seasons could have penetrated down to and have been retained above the blocking layer, resulting in the thick ice layer with a considerable water equivalent within the firn layer of 2009 and a strong IRH in the GPR signal.

Spatial extraction of recent accumulation rates
The firn core profiles of density, NH + 4 concentrations, and δD ratios suggested that, aside from the ice lens inherent to the 2009 accumulation layer, the IRHs can be attributed to dense summer layers.In fact, the related reflection in the GPR signal was observed in all analysed profile sections and was more pronounced at higher altitudes.The spatial con-tinuity of the ice layer is further supported by its evidence in both firn cores.By disregarding the respective IRH, the water equivalents of annual accumulation layers were extracted.As described above, the firn densification model was applied to all GPR traces individually where multiple IRHs were found that correspond to subsequent previous summer surfaces.Thus, the analysis is limited to areas where sufficient annual accumulation allows for distinguishing the corresponding IRHs.Furthermore, a continuous tracking of IRHs is impeded due to crevasses, ice dynamics, irregular accumulation such as avalanche deposits, changes in the distance of the helicopter to the snow surface, or lateral variations in layer properties that determine the strength of IRHs.The total length of analysed IRHs was 7.3, 6.8, 5.9, 2.8, 1.3, and 0.1 km for the accumulation layers of 2011 to 2006.
Conventional glaciological measurements refer to the annual surface mass balance, as they are taken from the recent accumulation layer in the end of the hydrological year.This layer experienced melting and the meltwater, which may partially refreeze deeper in the firn column, is not covered by the measurement.In contrast, water equivalents obtained from IRH travel times contain internal accumulation originating from surface melt in multiple seasons.As our goal is to reconstruct surface mass balance, consistent with conventional measurements, the computed quantity of refrozen meltwater was subtracted from each layer and was not re-attributed to its source layer.
The water equivalents derived from GPR exceeded those obtained by spatial extrapolation of conventional mass balance measurements, albeit reproducing their temporal pat-  tern (Fig. 5a).This can be attributed to the limited number of direct observations in the mass balance monitoring programme and the extrapolation scheme applied.Layer water equivalents derived from the firn cores revealed similar discrepancies with extrapolated mass balance.These findings highlight the need for a further investigation and better representation of accumulation processes in current glacier mass balance monitoring programmes.The spatial distribution of each annual layer's water equivalent showed considerable small-and medium-scale vari-ability, following the general topography to some degree (Fig. 6).The large-scale pattern of areas with increased annual accumulation rates recurs in all analysed layers.However, with increasing depth, older IRHs were more difficult to locate in the radargrams and appear in shorter profiles sections, thus hampering a direct temporal comparison of the accumulation patterns.Modelled layer densities increased with age as expected (Fig. 5b), reaching a mean value of 728 kg m −3 4 years after deposition.The cumulative effect of spatial variations in layer water equivalent, together with the fixed initial density in the densification model, results in a higher variability of the density of deeper, i.e. older, firn layers.Refreezing is a considerable contributor to firn densification on Alpine glaciers.It accounted for approx.10 % of the density increase over time and depth (Fig. 5b).At the analysed GPR measurement locations, the annual internal accumulation in the entire firn column is about 0.028 m w.e.yr −1 , affecting several previous accumulation layers.In summer 2011 the accumulation layers of 2010, 2009, and 2008 gained 0.019 m w.e., 0.007 m w.e., and 0.001 m w.e., respectively.On average, 2 % of the firn stratigraphy is made up by refrozen meltwater.In contrast, a higher value of 10 % was measured by Miller and Pelto (1999) on a temperate glacier in Alaska.This may be due to the different climatic condition but can also be related to processes not incorporated in the model, such as the formation of blocking layers and the storage or lateral flow above such barriers.
The conservative uncertainty estimate for the initial model density that covers measurement uncertainty, spatial and temporal variability, and the related uncertainty introduced within the model calibration affects the modelled layer densities and thus, the layer water equivalents by ±8 % (Fig. 5a).For conventional glaciological accumulation measurements the density measurement error is approx.4 % (±21 kg m −3 from repeat measurements).In contrast, the combination of GPR with a firn density model provides a considerable spatial coverage for several annual accumulation layers with a small trade-off in terms of uncertainty.

Data interpretation and error analysis
Generally, the formation of ice layers is not limited to surface processes.At boundaries with changes in the grain structure or on top of a blocking layer, e.g. a thin ice layer, meltwater can be retained (e.g., Pfeffer and Humphrey, 1996).Thus, an outstanding weather event during winter is not a necessary criterion to generate an ice lens that can potentially be misinterpreted as a former summer surface in the GPR signal.In the case of Findelengletscher, such a layer could be identified within the 2009 accumulation layer with the help of reference data from firn cores.Ultimately, the presented method of deriving past accumulation rates from GPR is based on the interpretation of GPR profiles by an observer, with the IRHs relating to high-density layers (Fig. 4).Therefore, the interpretation of GPR data can be ranked similar to a visual and density-based interpretation of a firn profile.Conventional glaciological measurements, comprising snow probings and snow pits for density measurements, are likewise prone to misinterpretations of the stratigraphy.This can be either visually in a snow pit, by hitting an internal ice layer, or by penetrating the designated summer surface during probing.
In contrast to the presented approach, conventional accumulation measurements are conducted on an annual basis, cover individual years, and thus do not require interpretation of multiple layers.Similar to other approaches which rely on annual layer counting, a misinterpretation of an IRH affects not only the estimated water equivalent of the respective accumulation layer but also the chronology of all the following layers is shifted; therefore, an independent evaluation of the annual character of IRHs is essential.We assessed potential alternatives to the in situ dating through firn cores by comparing the mean GPR-derived accumulation layer water equivalents to mass balance measurements and winter precipitation sums.Despite the small number of data points when comparing annual mean values (n = 5) and the resulting low statistical significance of the results, Pearson product-moment correlation was chosen to depict a potential linear relation between the data sets.
When the correct chronology from the firn cores is used to derive accumulation from IRH travel times, correlation with extrapolated mass balance measurements is r = 0.57, but insignificant (p = 0.32).If the ice lens within the 2009 accumulation layer was misinterpreted as a former summer surface, the water equivalents for the layers of 2009 and 2008 are strongly reduced and the dating of all following layers is shifted by 1 year (Fig. 7a).This results in a reduced corre-lation with the extrapolated mass balance of r = 0.41 (p = 0.50).Because mass balance measurements are not available at locations covered by analysed GPR profile sections this comparison relies on the extrapolation scheme used.We suggest that this explains the large overall deviation between the two data sets (Fig. 5) and the insignificance of the presented correlations.
A similar comparison was performed with precipitation sums from October to April at the Zermatt weather station.To avoid the effect of variations in spatial extents of analysed layers, only sections covering all layers of 2007-2011 were used.With the correct dating from cores, correlation of winter precipitation sums with GPR-derived water equivalents was r = 0.87 (p = 0.06) (Fig. 7b).In contrast, it reduces to r = 0.10 (p = 0.87) when the chronology is disturbed by the extra layer in 2009.
Additionally, the effect of a potential further perturbation of the layer chronology was tested by simulating an absence of individual or multiple IRHs in the GPR data.If the ice lens within the 2009 accumulation layer is kept and one other IRH is disregarded, correlations with measured mass balances and winter precipitation sums are all negative.This also holds for some perturbations if two of the detected IRHs are disregarded.The only case where this yields considerable positive correlations is if the ice lens within the 2009 layer is correctly disregarded and the subsequent IRH, representing the 2008 summer surface, would be missing.This yields a larger accumulation for 2009 than with the correct chronology, but does not conflict with the precipitation sums and measured mass balances.Therefore, only the bottom of the IRH sequence is disturbed and the further decrease in the number of data points generates implausibly high correlations.
The correlation of inferred winter accumulation with precipitation sums or measured mass balances requires a sufficient number of mostly annual IRHs and preferably a wide data range.In the given case it depends on the single large accumulation layer of 2009.Although such a comparison cannot provide a validation as robust as through firn cores -particularly if multiple layers are missing or misinterpretedwe suggest that it can aid an extraction of accumulation rates from GPR data when no firn cores are available.
To convert IRH travel times to depth and water equivalent, we apply an empirical relation of density and the dielectric permittivity of firn (Frolov and Macheret, 1999).However, this relation has been determined for dry snow conditions whereas the presence of liquid water has a considerable effect on the dielectric properties (Schmid et al., 2014).In our study, the analysed IRH travel times were derived from GPR measurements in spring (May 2012 andApril 2013).Field measurements of temperature in the snow cover of the accumulation area in April 2014 reveal subfreezing conditions apart from a shallow surface layer.Thus, the presence of liquid water in upper firn can be ruled out at this time of the year.On the other hand, residual water from previous melting seasons might be present deeper in the firn.This has a potential impact on the evolution of the temperature profile over the winter season because latent heat is released when it refreezes.More importantly, it affects the GPR propagation velocity and, thus, the resulting layer water equivalents.The potential impact was assessed using a rather conservative hypothetical liquid water content of 5 % vol.and the mixing formula for dielectric permittivity by Looyenga (1965).The related change in the resulting layer water equivalent is about −16 %, comprising the direct effect on the estimated layer thickness (Eq.5), as well as the effect on the modelled densities.However, comparison with the firn cores showed good agreement in the depths of IRHs (Fig. 4) and, thus, no evidence for the presence of liquid water.On the other hand, this finding cannot be transferred to the lower parts of the accumulation area, and does not necessarily hold for all investigated years.
The firn densification model is calibrated with the observed changes in IRH travel time from spring 2012 to 2013 at 12 locations and for 35 firn layers.Therefore, a different amount of liquid water at the date of the measurements could lead to an over-or underestimation of the compaction rate.This could, to some degree, account for the compaction rate scaling f = 2614 m −0.5 yr −0.5 being higher than previously published values.When the model was originally set up, Herron and Langway (1980) obtained f = 575 m −0.5 yr −0.5 for densities above 550 kg m −3 , using ice cores in cold firn on Greenland and Antarctica.Huss (2013) calibrated the same model with a set of 19 firn cores from temperate and polythermal mountain glaciers and ice caps in the European Alps, western and Arctic Canada, Central Asia, Patagonia, and Svalbard.The larger value of f = 1380 m −0.5 yr −0.5 is in line with the different climatic conditions and different firn compaction regimes compared to the original formulation.
Aside from the potential liquid water content, the scheme used to calibrate the model is affected by horizontal ice mo-tion.At the two accumulation measurement sites (Fig. 1), flow velocities are about 35 m yr −1 (north-eastern site) and 15 m yr −1 (south-western site).For the upper accumulation area, no direct observations of flow speed are available.However, due to a smaller ice thickness, velocities are expected to be considerably lower and to remain within the 25 m radius that was used to evaluate the mean IRH travel times.
In order to obtain a reliable change rate in the model we did not run the calibration for individual locations to account for spatial variations in the densification regime.Rather, the derived parameter set was assumed to represent the mean temporal and spatial condition for the investigated site.Comparison with the densities and layer water equivalents derived from the firn cores revealed a moderate overestimation of the local densities for the 2010 and 2009 accumulation layers resulting in a difference of 9 % for the water equivalents (Fig. 8).This also holds when hypothetical GPR travel times are calculated from the measured densities to run the model.Because the density estimate for the uppermost firn layer (2011) and, thus, the initial model density is adequate, this indicates an overestimation of the compaction rate.We suggest that this discrepancy is due to spatial variations in the density regime, related to processes not incorporated in the model, such as the snow grain structure, the temperature dependence of densification processes or particular weather conditions.
The water equivalent stored as firn between IRHs also differs from the accumulation component within the mass balance term regarding the definition of a "year".Mass balance monitoring approaches typically use the hydrologic year, with 30 September marking the end of the ablation period in the Northern Hemisphere.In contrast, the analysis of the GPR data is based on the firn stratigraphy, where the buildup of a high-density layer due to particular weather conditions provides the chronology of firn layers.Thus, a comparison or combination of such data sets requires converting between the stratigraphic and the fixed-date mass balance (Mayo et al., 1972;Huss et al., 2009), but does not generally conflict with the mass balance concept.With a coupled energy-balance-snow model the temporal evolution of accumulation layers can be tracked.Van Pelt et al. (2014) calibrated the accumulation input to such a model by matching modelled and observed IRH travel times.With this inverse approach they obtain the adjusted accumulation distribution along a GPR transect on a polythermal Svalbard glacier.Their model solves the surface energy balance and computes percolation, capillary storage, refreezing, and runoff of liquid water to provide subsurface density, temperature, and water content.Storage and lateral flow of water on top of ice layers are not resolved.Post-depositional processes related to refreezing are taken into account.Thereby, uncertainties in the model initialisation and physics directly affect the obtained accumulation patterns (Van Pelt et al., 2014).In this study, we used IRH travel times as input to a simple model for firn densification and meltwater refreezing.Thus, firn layer wa-ter equivalents could be derived in a forward approach.Our model involves strong assumptions on liquid water content and percolation, temperature evolution and densification of firn.By assuming that all subfreezing temperatures are compensated by meltwater refreezing we avoid modelling melt rates.This is reasonable under temperate conditions, but precludes an application to cold firn.On the other hand, a sound representation of physical processes depends on extensive input data such as air temperature, pressure, humidity, precipitation, and radiation (Mitterer et al., 2011).Finally, both approaches share the largest source of uncertainty which is the misinterpretation of IRHs as annual layers and a careful validation with reference data remains essential.

Conclusions
We established a new approach to derive past accumulation rates on temperate Alpine valley glaciers from repeated GPR surveys.Our method is based on the interpretation of IRHs as previous summer surfaces.By linking the measured GPR travel time of each layer to a simple model for firn densification, a density estimate was derived for each layer.Thereby, we take into account the refreezing of meltwater by modelling the end-of-winter temperature profile that, for temperate firn, is entirely compensated by refreezing.The model was calibrated by comparing the modelled and measured change in IRH travel times of 35 layers at locations where GPR measurements exist in 2 subsequent years.It was applied to each GPR measurement location where multiple IRHs were found to obtain the distribution of the mass of annual firn layers along the GPR profiles.
The strongest limitation to the methodology is the nonannual character of IRHs for our study site as it was found from the analysis of two shallow firn cores.A potential misinterpretation of IRHs affects the results for all subsequent layers by changing their ages.We suggest that, if no firn cores are available, measured mass balance data or winter precipitation sums could serve as a plausibility check.However, the formation of high-density and ice layers is linked to smallscale physical characteristics of snow and firn and is driven by particular weather conditions.Thus, a general conclusion about the validity of a dating-by-counting scheme for different study sites cannot be drawn.
We showed that GPR measurements can be used to estimate multiyear accumulation rates on a temperate valley glacier.Thereby, measurements can be conducted from a helicopter, facilitating efficient data acquisition even in inaccessible areas.Because the GPR penetrates several layers of accumulation, it could be used not only as a complement to existing monitoring programmes, but also for a retrospective extension of newly initiated time series.

Figure 1 .
Figure 1.Findelengletscher with the set-up of annual mass balance measurements.The mean equilibrium line altitude (ELA), the firn core locations, and all analysed GPR profiles measured in 2012, as well as the locations used for the model calibration are shown.The overview map shows its location within the European Alps.

Figure 2 .
Figure 2. Processed GPR profile exemplarily showing IRHs within the firn in the accumulation area of Findelengletscher.Red lines highlight potential previous summer surfaces.

Figure 3 .
Figure 3. Optimal scaling f of the compaction rate (Eqs. 2 and 3) as a function of the initial model density ρ f (0, t d ) = 492 kg m −3 , found by comparing 35 modelled and measured changes in IRH travel times at locations where GPR repeat measurements are available in 2012 and 2013.The shading shows the root mean square deviation (RMSD) while the line provides the optimal relation.The horizontal error bar shows the uncertainty in the initial model density.The related uncertainty in f is indicated with black markers.

Figure 4 .
Figure 4. Ground penetrating radar profile at firn core 1 (Fig. 1), the measured vertical profiles of density, transformed to the travel-time domain, with marked melt layers (orange) and dust layers (yellow), ammonium (NH + 4 ) and sulfate (SO 2− 4 ) concentrations, and the water stable isotope ratio (δD).Grey bars represent summer surfaces.For comparison, density, NH + 4 and δD are shown for firn core 2.
Figure 5. (a) Average water equivalents of annual layers derived by GPR and by extrapolated glaciological field measurements at locations where all four summer layers could be extracted from the GPR signal.Error bars show the uncertainty from the initial density ρ f (0, t d ).(b) Mean layer densities, the related uncertainty and data range for the same locations.

Figure 6 .
Figure 6.Annual accumulation layer water equivalent (m) derived along the GPR profiles (highlighted in black) for 2008-2011.For visualisation, the water equivalent was interpolated by inverse linear distance weighting, restricted to the 400 m surrounding of all GPR data points.

Figure 7 .
Figure 7.Comparison of GPR-derived mean annual layer water equivalents with (a) extrapolated mass balance measurements and (b) winter precipitation sums (October-April) at the weather station of Zermatt (1638 m a.s.l.).With the correct layer chronology obtained from firn cores (red dots), the high accumulation and winter precipitation in 2009 are reproduced.A misinterpretation of an IRH within the 2009 accumulation layer disturbs the chronology and affects all following layers (black circles, connected to the original values by stippled lines).

Figure 8 .
Figure 8. Layer water equivalents (a) and densities (b) at firn core 1, derived from GPR, directly from the firn core, and modelled using the hypothetical travel time obtained from the firn core layer water equivalent and density.Error bars show the uncertainty derived from variations in the initial density ρ f (0, t d ).