Introduction

Temporal changes in glacier mass balance are known to be an excellent indicator for changing climatic conditions (Reference Haeberli and BenistonHaeberli and Beniston, 1998; Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser and others, 2006; Reference SolomonSolomon and others, 2007). The impacts of future changes in climate forcing on mass balance and runoff from glacierized catchments are the subject of intense research (e.g. Reference HussHuss, 2011; Reference Farinotti, Usselmann, Huss, Bauder and FunkFarinotti and others, 2012; Reference Salzmann, Machguth and LinsbauerSalzmann and others, 2012). Models are used to assess these impacts and provide general conclusions from a small subset of monitored glaciers (Reference Zemp, Hoelzle and HaeberliZemp and others, 2009).

Presently, there is a wide range of mass-balance models available that cover ablation processes to all degrees of complexity (Reference HockHock, 2005, and references therein). To compute annual glacier mass balance all models need, aside from the ablation processes, to account for the amount and spatial distribution of snow accumulation. Typically, accumulation is estimated from precipitation measurements at nearby weather stations, gridded precipitation climatologies or climate model outputs that are tuned to site characteristics, usually by applying an elevation gradient and temperature threshold for snowfall (Reference HockHock, 1999; Reference Machguth, Paul, Hoelzle and HaeberliMachguth and others, 2006a, Reference Machguth, Paul, Kotlarski and Hoelzle2009). However, this is a weak representation of the actual snow accumulation pattern. Aside from avalanches, the interplay of topography and wind leads to strong variations in the spatial distribution of snow depth (Reference Liston and SturmListon and Sturm, 1998; Reference Winstral and MarksWinstral and Marks, 2002; Reference Machguth, Eisen, Paul and HoelzleMachguth and others, 2006b; Reference Lehning, Löwe, Ryser and RadeschallLehning and others, 2008; Reference Dadic, Mott, Lehning and BurlandoDadic and others, 2010a). Due to the complexity of the physical modelling of preferential snow deposition and redistribution processes, and due to the difficulties inherent in the measurement of solid precipitation in high-alpine terrain, the accumulation part of state-of-the-art mass-balance models is much less developed than the description of melting processes (e.g. Reference Dadic, Mott, Lehning and BurlandoDadic and others, 2010a). Improving the level of understanding with regard to the accumulation of snow on alpine glaciers is a prerequisite for strengthening the physical justification of the accumulation modules of mass-balance models (Reference Machguth, Eisen, Paul and HoelzleMachguth and others, 2006b). To date, only a few studies have tackled this issue by providing the necessary extensive end-of-winter snow depth distribution measurements (Reference Machguth, Eisen, Paul and HoelzleMachguth and others, 2006b; Reference Dadic, Mott, Lehning and BurlandoDadic and others, 2010a,Reference Dadic, Mott, Lehning and Burlandob). However, evaluating the accuracy of the estimated spatial distribution requires separate reference data to be available.

The spatial distribution of snow depth is commonly measured by time-consuming in situ probings that provide point values and can have large uncertainties in accumulation zones where locating the previous summer surface can be challenging. Other possible approaches are ground-penetrating radar (GPR) and differencing of digital elevation models (DEMs). GPR provides a non-destructive method for various purposes in glaciology. This is supported by the low conductivity and, hence, the deep penetration of the signal into snow and ice at MHz frequencies. The signal is reflected from boundaries with a contrast in dielectric permittivity that can be related to a change in material properties (e.g. from ice to bedrock, variations in density, impurities such as dust or ash, or water content (Reference UlriksenUlriksen, 1982; Reference Plewes and HubbardPlewes and Hubbard, 2001)). Since the 1960s this has been used for ground-based or airborne detection of ice thickness (Reference Robin, Evans and BaileyRobin and others, 1969; Reference Bauder, Funk and GudmundssonBauder and others, 2003; Reference DammDamm, 2004; Reference Lalumiere, Prinsenberg, Chung, Prinsenberg, Hong and NagataLalumiere and Prinsenberg, 2009) or internal layers within snow and firn (Reference Kohler, Moore, Kennett, Engeset and ElvehøyKohler and others, 1997; Reference Kanagaratnam, Gogieneni, Ramasami and BraatenKanagaratnam and others, 2004; Reference Dunse, Schuler, Hagen, Eiken, Brandt and HøgdaDunse and others, 2009; Reference Heilig, Eisen and SchneebeliHeilig and others, 2010). Radar-based surveys of snow accumulation distribution have focused on ice sheets (Reference HelmHelm and others, 2007; Reference Dunse, Eisen, Helm, Rack, Steinhage and ParryDunse and others, 2008; Reference EisenEisen and others, 2008; Reference Kruetzmann, Rack, McDonald and GeorgeKruetzmann and others, 2011) and snow on the ground (Reference Marchand, Killingtveit, Wilén and WikströmMarchand and others, 2003; Reference Yankielun, Rosenthal and DavisYankielun and others, 2004; Reference Marshall and KohMarshall and Koh, 2008). For mapping snow accumulation on alpine glaciers, however, the airborne application of GPR has only been addressed by Reference Machguth, Eisen, Paul and HoelzleMachguth and others (2006b), although results show that it can be a highly effective tool, especially in remote and inaccessible areas.

More attention has been paid to the use of DEMdifferencing to assess changes in glacier mass balance (Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt and others, 2002; Reference Bamber and RiveraBamber and Rivera, 2007; Reference Abermann, Fischer, Lambrecht and GeistAbermann and others, 2010; Reference Kääb, Berthier, Nuth, Gardelle and ArnaudKääb and others, 2012). Light detection and ranging (lidar) scanning systems offer a high-precision tool for range measurements, that has been applied for ∼20 years. The travel time of a laser pulse from emission until reflection from a target and return to a detector, coupled with the speed of light in the air conditions, provides a measure for the distance between the sensor and the illuminated spot. In airborne applications, the exact position and orientation of the aircraft are necessary to locate the ground points, which are then interpolated to the regular grid of a DEM (Reference AckermannAckermann, 1999; Reference Wehr and LohrWehr and Lohr, 1999). In our study, repeated lidar DEMs of Findelengletscher, Switzerland (Reference Joerg, Morsdorf and ZempJoerg and others, 2012), provide spatially distributed elevation changes between autumn and spring. To obtain distributed snow depth over glaciated surfaces, this requires corrections of the vertical ice velocity and firn compaction in the surface differencing scheme. Such corrections have not been realized so far, although studies indicate their need (Reference Dadic, Mott, Lehning and BurlandoDadic and others, 2010a,Reference Dadic, Mott, Lehning and Burlandob; Reference Helfricht, Schöber, Seiser, Fischer, Stötter and KuhnHelfricht and others, 2012, Reference Helfricht, Kuhn, Keuschnig and Heilig2013).

In this study we discuss simultaneous measurements of snow accumulation distribution on an alpine valley glacier, using the traditional glaciological approach of in situ snow probings and snow density pits, helicopter-borne GPR measurements of the snow depth and surface elevation changes given by differencing of two lidar DEMs. Snow probings and DEM-differencing provide the winter mass balance of Findelengletscher. Furthermore, we show how the geodetically derived snow depth distribution can be corrected for the vertical ice velocity and firn compaction, based on the in situ probings and long-term mass-balance and geometry changes. The result is a fully distributed dataset that is validated with the GPR measurements and correlated against elevation, surface slope and curvature. Analysis of the three datasets provides new insights into the ability of the utilized methods to correctly represent snow depth, and contributes to a better understanding of processes that govern winter accumulation of alpine glaciers.

Study Site and Field Data

Findelengletscher is a large temperate valley glacier (13.0 km²) in Valais, Switzerland. It is northwesterly exposed and covers an elevation range from 2600 to 3900 m a.s.l. (Fig. 1). Situated directly below the main ridge, local precipitation largely depends on the wind direction. The accumulation area strongly benefits from south- and southeasterly conditions. A glacier mass-balance monitoring programme was started in 2004 (Reference MachguthMachguth, 2008). Annual mass balance is determined with a network of 13 ablation stakes and two snow density pits. On average, the measurements indicate a negative mass balance (Reference Joerg, Morsdorf and ZempJoerg and others, 2012), leading to a continued retreat of the glacier snout (Glaciological Reports, 2011). Glaciological measurements of seasonal mass balance from snow probings and density profiles combined with repeated terrestrial snowline photographs are available from 2009 (Reference HussHuss and others, 2013).

Fig. 1. Map of the study site: Findelengletscher, Valais, Switzerland. Dots show manual snow probings, lines represent GPR profiles, the colour code shows measured snow depth in April 2010.

The conventional method of vertically probing the underlying ice (or last summer) surface through the snowpack provides a direct measurement of snow depth as a point value (Reference Østrem and BrugmanØstrem and Brugman, 1991). The probe measurements were carried out on 10–11 April 2010. For practical reasons, they were taken along the walking paths leading from the upper accumulation area to the tongue of Findelengletscher (Fig. 1). In its accumulation area the boundary between snow and firn is typically made up of a clearly defined ice lens of variable thickness, that often can be unambiguously detected with a probe. However, careful probing is necessary to avoid misinterpretation of internal ice layers within the snow cover. The entire glacier elevation range was covered by 403 probings (Fig. 2), each of which was made up of two to five averaged probings to account for outliers (crevasses) and small-scale variability. Mean snow depth from all probings was 1.86 m, and the standard deviation within multiple measurements at a given location was 0.13 m. The measured snow depth was elevation-dependent but highly variable (Fig. 2).

Fig. 2. Measured snow depth from 403 probings in April 2010.

Additionally, the bulk density of a snow column was measured in 11 snow pits across the glacier. Values ranged from 223 to 294 kg m⁻³ on the glacier tongue to 407 kg m⁻³ at the centre of the glacier (3130 m a.s.l.). This probably represents (1) different process regimes in snowpack evolution, involving wind effects on initial deposition, erosion and redeposition of snow and/or (2) differences in the prevalent processes of snow metamorphism. While the internal vertical variability of the snow density was not covered by the sampling scheme, multiple bulk measurements were taken at ten of the snow pits. Standard deviation within measurements at these locations was 21 kg m⁻³. Density was averaged over 100 m elevation bands. For bands that were not covered by at least one measurement, the average of neighbouring bands was used. Weighting density according to the area of the respective elevation bands revealed a spatial mean snow density of 362 kg m⁻³. To follow the ongoing monitoring programme of seasonal mass balance of Findelengletscher we applied an enhanced traditional extrapolation scheme to the snow probings, to obtain snow depth on a 25 m × 25 m resolution grid. It is based on a distance-weighted average of probings within a radius of 400 m around each gridpoint. If fewer than three measurements were available, the search radius was increased. Additionally, we superimposed small-scale snow depth variability due to local effects that might not be captured by the spatially discontinuous probings. According to a procedure described by Reference Huss, Bauder, Funk and HockHuss and others (2008), snow depths were decreased over steep slopes (>40°) and convex terrain shapes, whereas they were increased for concave surfaces. These corrections are relatively small for the central parts of the glacier, but are important to yield realistic snow depths for regions prone to avalanching and wind erosion.

The GPR system was a Noggin 500 (Sensors & Software Inc., Canada) that operated a shielded antenna with a centre frequency of 500 MHz. It was directly mounted on the bottom of a helicopter. Traces were recorded with a constant time increment of 0.1 s at a waveform sampling interval of 0.2 ns. Their position was provided by a synchronized GPS with a horizontal accuracy of 1–5 m. On 10–11 April, ∼12 000 traces were recorded with the helicopter flying at ∼10 m s⁻¹ and 5–10 m above ground. The angle of beam spread is 45° and 60° perpendicular and along the flight direction, respectively. Hence, the footprint size at the snow surface is less than 8 m × 12 m. The survey covered 12.7 km of linear tracks in ten profiles between 2900 and 3550 m a.s.l. across the glacier (Fig. 1). In an airborne application of the GPR system the risk of interference with other essential and safety-of-life radio services is high. A special licence for the campaign was granted by the Swiss telecommunication regulator (OFCOM).

Lidar surveys were carried out on 4 October 2009, while the ablation area of Findelengletscher (<3200 m a.s.l.) was snow-free, and on 10 April 2010 at the end of the accumulation season. The Optech ALTM 3100 scanning system was mounted to a Pilatus Porter fixed-wing aircraft equipped with differential GPS and an inertial measurement unit to derive its position and orientation. Average point density was 7.6 and 8.1 m⁻², respectively. Reference Joerg, Morsdorf and ZempJoerg and others (2012) describe the geo-referencing, point cloud interpolation and generation of the grid (1 m × 1 m cell size) in further detail and provide an extensive uncertainty assessment. Standard deviation on reference surfaces was 0.09 m (October 2009) and 0.19 m (April 2010, reference surfaces snow-covered). Additionally, longer-term geometry change is given by two lidar DEMs from similar surveys on 28–29 October 2005 and 29 September 2010.

Methods

Processing of GPR data

Many processing techniques for GPR data have been adopted from seismic reflection surveys that have a similar principle (Reference Davis and AnnanDavis and Annan, 1989; Reference Fisher, Stewart and JolFisher and others, 1996). Since the selection of processing steps and their parameter settings depend strongly on field conditions, survey design and intention, no universal procedure is available (Reference UlriksenUlriksen, 1982; Reference Fisher, Stewart and JolFisher and others, 1996; Reference AnnanAnnan, 1999). In the given order, the following processing steps were found to sufficiently improve the helicopter-borne GPR data quality:

1. 1. Spatial interpolation; correction for variations in the helicopter’s velocity. Traces were interpolated linearly to an equidistant spacing of 1 m.

2. 2. Frequency bandpass 300–1000 MHz; the high-pass filter removed bending and constant offset of the signal while the low-pass filter removed high-frequency speckle and noise (Reference JolJol, 2009).

3. 3. Background removal; subtracting the average over the distance domain from all traces in the profile reduced the effects of signal ringing and the airwave.

4. 4. Gain function, linear and exponential; to account for signal attenuation with depth.

Migration did not improve the data quality and was not applied. For all measured radar profiles these processing steps allowed the the air/snow interface and the snow/firn or snow/ice interface to be distinguished without the help of separate snow depth data (Fig. 3). The reflectors were marked in the processing software (Reflexw, Sandmeier Scientific Software) either manually or semi-automatically with a built-in phase-following tool. To convert the difference in two-way travel time between the two marked reflectors from the time to the depth domain, an assumption for the radar wave velocity is required. Density and liquid water content strongly affect the dielectric properties of snow and several empirical formulas exist to compute the permittivity (Reference Kovacs, Gow and MoreyKovacs and others, 1993; Reference Frolov and MacheretFrolov and Macheret, 1999). According to Reference Frolov and MacheretFrolov and Macheret (1999), we obtained the relative permittivity, , of dry snow from the mean measured density, ρ = 362 kg m⁻³, as

(1)

Fig. 3. First 1500 m of GPR profile 1 which lie in the ablation area (Fig. 1) after processing. (a) The reflection of the interface between air and snow. (b) The reflection of the interface between snow and glacier ice. The green markers indicate a crevassed area in the centre of the glacier.

The mean radar wave velocity, u, was then computed as . This constant value was used to derive snow depth from the reflectors in the GPR profiles. However, uncertainties and variations in snow density directly affect u and, therefore, the conversion from travel time to snow depth. Further, we assume that the snowpack was dry. This is in line with observations during fieldwork, except for the glacier tongue that was not covered by GPR profiles. An analysis of the related errors is given in the Discussion.

The 500MHz antenna also resolved some internal layering of the snowpack. In the accumulation area of Findelengletscher, several firn layers from previous years could be distinguished (Fig. 4). While the 2009 surface could still be determined consistently, the successive dating of the earlier summer surfaces remained speculative. The penetration depth of the GPR signal depends on the height of the helicopter above ground and the physical properties of the subsurface. In the accumulation area we found reflections from down to 20 m below the snow surface.

Fig. 4. GPR profile 2 which lies in the accumulation area (Fig. 1), after processing with reflectors of the snow surface (marked) and (a–d) annual summer surfaces.

Geodetic elevation change from lidar DEMs

In a first-order approach, the elevation change of each gridcell was interpreted as snow accumulation (Reference Grünewald, Schirmer, Mott and LehningGrünewald and others, 2010; Reference Egli, Jonas, Grünewald, Schirmer and BurlandoEgli and others, 2012; Reference Helfricht, Schöber, Seiser, Fischer, Stötter and KuhnHelfricht and others, 2012). The resulting pattern of maximal snow depths on the glacier tongue and small snow depths in the accumulation area (Fig. 5) was, however, inconsistent with in situ snow probings that revealed increasing snow depth with elevation (Fig. 2). The principle of comparing two DEMs of different dates makes this approach sensitive to secondary processes that alter the glacier surface elevation. This involves (1) the vertical displacement of the ice surface as part of the glacier flow, hereafter referred to as the emergence velocity, which is positive (emergence) in the ablation area and negative (submergence) in the accumulation area, (2) surface lowering due to firn compaction in the accumulation area and (3) surface lowering from melt after the acquisition of the October 2009 DEM. The formation of superimposed ice on the summer surface is assumed to be negligible in our case, due to the temperate nature of Findelengletscher.

Fig. 5. Raw lidar-derived elevation change from October 2009 to April 2010. Grey shading indicates negative values.

The ice emergence velocity affects the measured surface elevation change as the glacier moves during the time-span between the two generated DEMs (Reference Helfricht, Schöber, Seiser, Fischer, Stötter and KuhnHelfricht and others, 2012, Reference Helfricht, Kuhn, Keuschnig and Heilig2013). For a glacier that is in balance with climatological and all other forcings, this vertical component of glacier flow is equal to the reversed local annual mass-balance rate, −b _a, divided by the density of ice, ρ _ice. In the accumulation area, the compaction of firn, v _firn, adds to the glacier surface elevation change. A potential imbalance of mass balance and ice flow can be accounted for with the observed annual geometry changes, Δz. The annual vertical displacement, v _z, of the last summer surface is derived for fixed points, (x, y), on the glacier surface. Thus, it is not affected by horizontal flow and slope (Reference Cuffey and PatersonCuffey and Paterson, 2010) and can be written as

(2)

where positive values point upwards (emergence), and negative values indicate a lowering of the surface (submergence). Over the cycle of one hydrological year and integrated over the whole firn column, one annual accumulation layer is transformed from snow to ice. At the top of the firn column this results in a surface lowering of

(3)

with ρ ₀ as the initial firn density. Equations (2) and (3) can be combined to give

(4)

with ρ = ρ _ice in the ablation area and ρ = ρ ₀ in the accumulation area. For the initial firn density we set ρ ₀ = 550 kg m⁻³ as an average from repeated measurements in previous years. Density of glacier ice was set to ρ _ice = 900 kg m⁻³.

The annual mass balance of every cell of a 25 m × 25 m resolution grid for the hydrological years 2005/06 to 2009/10 was evaluated using a distributed accumulation and melt model (Reference Huss, Bauder and FunkHuss and others, 2009), that was tuned to the in situ mass-balance measurements of each year individually. In this model, the spatial variations in annual mass balance for unmeasured regions of the glacier are given by a description of physical processes (Fig. 6a). The geodetic geometry change was available for the same period from differencing of two lidar DEMs. By annually averaging the mass balance and geometry changes for all years of the geodetic survey, and were obtained. The spatially integrated estimation of the emergence velocity is equal to zero when the effect from firn compaction is omitted (Fig. 6b). This is a result of the cumulative annual mass balance, b _a, corresponding to the geodetic mass change, Δz, in Eqn (2). This estimation refers to the displacement of the summer surface, so a spatial average different from zero would imply a glacier volume change due to ice flow.

Fig. 6. (a) Mean annual mass balance 2005/06 to 2009/10. (b) Vertical component of ice flow (emergence velocity), derived from the mean annual mass balance and observed geometry changes. A positive velocity is directed upwards (emergence). (c) Modelled melt that occurred after the lidar DEM was generated in October 2009.

Because the time-step from October 2009 to April 2010 corresponds to a fraction rather than a full year, the annual average emergence velocity must be scaled to fit this time-span. However, the choice of this scaling must integrate temporal variations of the respective processes. Changes in the glacial drainage system affect ice-flow characteristics and are known to show seasonality (Reference Iken and TrufferIken and Truffer, 1997). Furthermore, the densification of firn depends on its temperature and the accumulation in subsequent years. Thus, we introduced a tuning factor, f = 0.795, to account for seasonal and annual variations in ice-flow velocity and firn compaction rates. This factor f was derived by minimizing the square error of the resulting snow depth to the in situ snow probings (see Discussion).

The melt occurring after the lidar survey in October 2009 was estimated using a distributed temperature-index model (Reference HockHock, 1999; Reference Huss, Bauder and FunkHuss and others, 2009). The model was driven by daily temperature and precipitation obtained from a weather station at Zermatt (6 km distance from the glacier tongue, 1638 m a.s.l.). The resulting cumulative melt, m (m ice eq.), was added to the difference, Δh, of the two lidar DEMs in October and April. On the glacier tongue it had the highest impact of up to 0.3 m (Fig. 6c). Together with the scaled correction for emergence velocity and firn compaction, v _z, this provides a corrected snow depth, d, (Fig. 7) as

(5)

Fig. 7. Snow depth derived from lidar DEMs corrected for emergence velocity, firn compaction and autumn melt (Eqn (5)). Grey shading indicates negative values.

The elevation-dependent distribution of snow depth derived from the extrapolated snow probings and lidardifferencing underlines the need for this correction (Fig. 8). The uncorrected lidar-derived snow depth is strongly affected by emergence velocity and firn compaction. After application of the correction the two methods show good agreement in most of the elevation bands.

Fig. 8. Vertical distribution of mean snow depth in 100 m elevation bands derived from extrapolated snow probings, raw lidardifferencing and lidar-differencing corrected for emergence velocity, firn compaction and autumn melt (Eqn (5)). The distribution of the glaciated area is shown in grey.

The presented method to derive the snow accumulation distribution from lidar-differencing requires, in addition to the two DEMs, long-term annual mass-balance and geometry changes and in situ measurements of snow depth. This expense can be reduced if the geometry changes can be omitted, implying that the glacier is in equilibrium. Furthermore, information about the vertical displacement of the glacier surface could be derived, for example, from GPS measurements or triangulation of the positions of the commonly used mass-balance stakes. If the vertical displacement of the last summer surface during the winter is known, at least for a few points, the scaling of the correction and its calibration with in situ snow probings are redundant. Otherwise, an estimation of snow depth is necessary to scale the correction to fit the time-span between the generation of the two DEMs. To provide a stable fit, these estimations should preferably exist for areas where the annual emergence velocity is high, i.e. on the glacier tongue and in the upper accumulation area, but are not required to cover the glacier representatively.

Results

The three datasets differ in the type of spatial information they present. GPR profiles provide snow depth in linear transects and do not cover the entire glacier representatively. To compare GPR with in situ probings and lidar-derived snow depth, the probings must be extrapolated to achieve a least common spatial coverage of all three datasets. None of them is considered to be the ground reference, to allow full advantage to be taken of their simultaneous availability.

All three methods reveal the increase in snow depth with elevation. They disagree in crevassed areas, where the determination of a meaningful snow depth value is hampered by its large variability within the GPR footprint and movement of the crevasses with time that affects the local lidar-derived snow depth (Fig. 9). Both lidar and GPR are able to resolve the small-scale variability of snow depth, to a certain degree, due to their spatial resolution. Obviously, this is not the case for the spatially extrapolated snow probings (Fig. 9). Along the GPR profiles, with a constant spacing of 1 m between measurements, the rootmean-square error (rmse) and bias in the datasets is calculated using the closest gridpoints of the snow depth maps from lidar and extrapolated probings. In this comparison, the extrapolated snow probings reveal a higher, and lidar-differencing a lower, snow depth than GPR (Table 1). The generally high rmse of 0.39–0.62 m between the methods can partly be explained by the extrapolation scheme for the snow probings and the mass balance underlying the correction of lidar snow depth. Further, due to the high resolution of the lidar DEMs, the horizontal ice flow adds a local stochastic uncertainty by shifting small-scale surface features, including crevasses, to neighbouring gridcells. The smallest rmse was found to be 0.39 m between the extrapolated probings and lidar-derived snow depth. This is probably because of the fitting of the correction to the snow probings. Differences of the in situ snow probings (not extrapolated) to GPR and lidar will be highlighted in the Discussion (Figs 12 and 13).

Fig. 9. Snow depth obtained from GPR and at the nearest gridcells of the extrapolated snow probings and corrected lidar DEM-differencing along the first 1500 m of GPR profile 1, ranging from 2880 to 3080 m a.s.l. (Fig. 1).

Table 1. Root-mean-square error (rmse) and bias of extrapolated snow probings (ESP), lidar and GPR-derived snow depth along the GPR profile lines (crevassed areas masked, profile lines shown in Fig. 1)

The traditional glaciological approach, using extrapolated snow probings and snow density in the elevation bands, provides a winter mass balance of Findelengletscher of 0.78 m w.e. The lidar-differencing yields a similar result of 0.72 m w.e. Although the point snow probings were used for the extrapolation and the scaling of the correction for ice emergence velocity, the good agreement in winter mass balance underlines the high potential of deriving snow depth by lidar DEM-differencing.

The extrapolated snow probings allow a spatial comparison with the lidar-derived snow depth distribution to be made. Because the in situ snow probings were used to find the optimal scaling factor, f, of the lidar correction, the two datasets are not fully independent. However, the linear scaling of the correction (Eqn (5)), which has a strong altitudinal gradient (Fig. 6b), affects the distribution of snow between the ablation and accumulation areas and, only to a lesser extent, its variability at a smaller scale. Furthermore, a similar value of f would have been found using the GPR snow depth instead of the in situ probings (see Discussion). Thus, the dependence of the two datasets plays a minor role when comparing the spatial snow accumulation distribution. The two methods show good agreement for large parts of the glacier (Fig. 10). They disagree on snow depth for marginal regions of the glacier, especially in areas that were not covered by snow probings or mass-balance measurements and, thus, were not adequately represented in the respective spatial extrapolation. Particularly, this involves parts of the accumulation area and the uppermost part of the glacier, where wind erosion of snow is most pronounced.

Fig. 10. Difference of corrected lidar-derived to extrapolated snow depth from probings. Blue indicates an underestimation, and red an overestimation, of lidar-derived snow depth compared with the extrapolated probings.

The availability of distributed snow depth data and high-resolution DEMs allowed us to investigate snow depth distribution with surface topography, such as elevation, slope and curvature. While elevation affects the amount of precipitation, slope and curvature were taken as proxies for local variations in snow depth due to preferential deposition and redistribution (Reference Huss, Bauder, Funk and HockHuss and others, 2008; Reference Lehning, Löwe, Ryser and RadeschallLehning and others, 2008). We followed previous studies that used multiple linear regression (Reference Schmidt, Weber and WinigerSchmidt and others, 2009; Reference Farinotti, Magnusson, Huss and BauderFarinotti and others, 2010) in the form

(6)

to estimate snow depth, d. All grids were upscaled to a 10 m × 10 m cell size to reduce computation time. The October 2009 lidar DEM provided elevation, x ₁, and allowed computation of slope, x ₂, and curvature, x ₃, (Reference Zevenbergen and ThorneZevenbergen and Thorne, 1987). The coefficients β_i were found by least-squares estimation from the lidar-derived distributed snow depth (Fig. 7). Curvature was derived within a radius of 70 m around each gridcell. Although not supported by any physical process, this radius provided the highest correlation with snow depth. A positive sign indicates a convex surface. Crevassed areas were masked because the lidar-derived snow depth gives unreliable results in these regions.

Over the entire glacier, correlation of snow depth with elevation is r = 0:40, while slope (95% central range: 3–31°) has little influence (r = 0:20). Higher snow depths are found in areas with a concave curvature (r = −0:27). This multiple linear regression model can only partially explain snow depth variability (r = 0:51, p < 0:001 for all given correlations). This is in line with the findings of Reference Farinotti, Magnusson, Huss and BauderFarinotti and others (2010). They point out that this cannot be interpreted as absence of an influence of these factors on snow depth, as it can be nonlinear or masked by other factors, such as variations in wind speed and its prevailing direction (Reference Purves, Barton, Mackaness and SugdenPurves and others, 1998; Reference Plattner, Braun and BrenningPlattner and others, 2006).

From previous studies it is known that deposition and redistribution of snow is largely determined by the interplay of wind and surrounding topography (Reference Liston and SturmListon and Sturm, 1998; Reference Winstral and MarksWinstral and Marks, 2002; Reference Lehning, Löwe, Ryser and RadeschallLehning and others, 2008; Reference Dadic, Mott, Lehning and BurlandoDadic and others, 2010a). In the regression this was incorporated as curvature on the small scale of a 70 m radius. Obviously, this approach does not respect actual local wind fields that are also formed by surface features on larger scales. To account for spatial variations in the influence of the selected variables, we introduced a zonal regression using the same variables but restricted to a 2 km × 2 km moving box around each gridcell of the lidarderived snow depth map. By running a separate regression for the surroundings of each gridcell, we derived a smoothed map of the strength of the linear dependence of snow depth on the selected variables (Fig. 11).

Fig. 11. Correlation of lidar-derived snow depth with (a) elevation, (b) curvature (70 m scale) and (c) regression-derived snow depth for the surrounding 2 km × 2 km box of each gridcell.

The zonal regression shows a better performance than the regression for the entire glacier and allows a more detailed interpretation. The multiple linear correlation coefficient ranges from r = 0:22 on the upper glacier tongue (2900 m a.s.l.) and in the lower accumulation area (3250–3350 m a.s.l.) to r = 0:83 in the upper ablation area (3100–3200 m a.s.l.) and upper accumulation area (>3450 m a.s.l.) (Fig. 11c). If instead the same number of gridcells is randomly chosen from across the glacier, the multiple linear correlation coefficient is r = 0:45 ± 0:01. The size of the moving box determines the smoothness and the data range of the correlation coefficient maps in the expected way (e.g. 0:18 ≤ r ≤ 0:94 for 1 km × 1 km, 0:30 ≤ r ≤ 0:83 for 3 km × 3 km). Thus, the variable correlation coefficients for the zonal regression can be attributed to spatial variations in the relation of snow depth to the examined topographic variables. While the upper glacier tongue might show overestimation of lidar-derived snow depth, as indicated by GPR and snow probings (Figs 9 and 10), the weak performance in the lower accumulation area is related to the changing sign in the correlation with elevation. Zonal correlation with elevation (Fig. 11a) highlights that a linear elevation gradient is not representative for the entire glacier. Instead, the correlation of elevation and snow depth changes sign at high altitudes. For slope, the correlation is less than lrl = 0:38 for 75% of the glaciated area and also changes to a negative sign at high altitudes (p < 0:001 for all correlations). At the same time, a regression over the whole of this area of negative zonal correlation with slope gives r = −0.03. Curvature yields strictly negative correlations that are generally stronger at higher altitudes (Fig. 11b). However, the effect of curvature on snow depth is not restricted to the 70 m scale that is used here and can thus be underestimated by the regression.

For elevation, slope and curvature, the effect on snow depth is not constant over the entire glacier. This is likely to be linked to variations in wind speed and exposure that are controlled by the surrounding topography. Therefore, elevation, slope and fixed-scale curvature alone cannot sufficiently explain snow depth distribution in a linear regression model for the entire glacier, but can within subregions of similar influence of wind fields.

Discussion

Systematic bias in the snow probings dataset

In contrast to the snow probings, GPR provides snow depth as continuous profiles, instead of discrete point values. Thus, we evaluated a potential difference of the snow probings from the GPR data as a function of distance between the measurements. Therefore, all possible pairs of one measurement from each method were identified. Mean and standard deviation of the difference between every two measurements were obtained for all pairs with less than a critical distance, increased in 10 m steps. This allows a comparison of the two datasets to be made, despite the variable distance between the snow probings and the GPR profiles.

A good agreement between the two datasets is found where the distance between snow probings and GPR profiles is small (Fig. 12), and variability is of the same order within multiple probings as at individual locations (0.13 m). However, variability not only increases with distance, but the snow probings exceed the GPR-derived snow depth. This becomes important at distances >30 m and indicates that the analysed subsets of the two methods cover different regions of snow accumulation. We interpret this distance-dependent deviation between the datasets as a systematic bias in the in situ measurements that might be explained by either one or both of these processes: (1) the GPR profiles cover parts of the glacier surface that are inaccessible and are therefore not represented in the snow probings and/or (2) the walking path corresponding to the snow probings inadvertently avoided ridges and steeper parts of the glacier that are prone to wind erosion and thus shallower snow depth, and rather was attracted by depressions, shallow troughs or flat areas in which snow tends to be deposited.

Fig. 12. Distance-dependent deviation of GPR-derived snow depth for probings, shown as mean and standard deviation of all measured deviations within each 10 m class. Each class contains all possible measurement pairs of the two datasets with a distance shorter than the upper class limit.

Uncertainties in GPR-derived snow depth

GPR provides estimates of snow depth that are not affected by glacier dynamics, such as firn compaction or ice flow. Nevertheless, uncertainties that arise from the choice of radio-wave velocity and the practical vertical resolution of the radargrams need to be taken into account. Also, it should be noted that the accuracy of the observer when digitizing the reflectors could play a role that cannot be exactly quantified.

The resolution of the GPR data limits the accuracy of recognized reflector travel times. It depends on the signal bandwidth, that can be approximated by the centre frequency for most GPR systems. It is typically estimated to be about half the wavelength, λ, but depends on the signal-to-noise ratio (Reference JolJol, 2009). Thus, in the given case with a frequency of 500 MHz and a velocity of 0.229 m ns⁻¹, the vertical resolution is ∼ λ/2 = 0.23 m.

As shown above, the velocity of the radar pulse depends on the characteristics of the snowpack. For snow density, the standard deviation of 21 kg m⁻³ from snow pits where more than one density profile is available is used as an estimate of the measurement accuracy. The total standard deviation of measured density from the area-weighted mean (68 kg m⁻³) accounts for uncertainties in the density extrapolation method and the use of a constant density. Using Eqn (1), the related total standard error in radio-wave velocity is ±0.011 m ns⁻¹. By linearity, the effect on the resulting snow depth then is ±4.7%.

For the computation of radio-wave velocity we assumed the snowpack to be dry along the GPR profiles. Although this was in line with observations during fieldwork, from guideline values (e.g. Reference FierzFierz and others, 2009) a water content of <3%vol. cannot be ruled out. Due to the high permittivity of water, this could have a considerable influence on the radio-wave velocity in the snowpack. Based on numerous sets of measurements, Reference Frolov and MacheretFrolov and Macheret (1999) compute the increase in dielectric permittivity, Δε′, due to the relative water content, θ, as

(7)

Thus, liquid water content of up to 3% could reduce the radio-wave velocity from 0.229 m ns⁻¹ to 0.199 m ns⁻¹ by causing a higher dielectric permittivity of the snowpack. For the resulting snow depth this relative change of up to −13.1% could make a significant difference but cannot be determined exactly. Furthermore, the water content varies not only with depth in the snowpack but also spatially (Reference Techel and PielmeierTechel and Pielmeier, 2011). This makes it difficult to incorporate this factor in the analysis. During fieldwork, the snowpack in the lower parts of the glacier (i.e. the tongue) showed evidence of some liquid water content, whereas the snow pits in the elevation range covered by GPR profiles were dry, due to lower temperatures. Furthermore, the velocity of 0.229 m ns⁻¹ which is used in this study indicates good agreement with the in situ snow probings (Fig. 12). A reduced velocity would even increase the already negative deviation from the extrapolated snow probings (Table 1). Thus, the effect is likely to be small.

Ice-flow dynamics affect lidar-derived snow depth

From the vertical distribution of mean snow depth (Fig. 8) it is obvious that the correction for emergence velocity and firn compaction is necessary when evaluating snow depth distribution from repeated lidar DEMs. As described, the resulting vertical velocity of the autumn glacier surface was estimated as an annual average from long-term mass-balance and geometry changes. Since the temporal difference of the DEMs did not correspond to a full year and to account for interannual variations, the correction was scaled by a factor f (Eqn (5)).

Earlier studies interpreted seasonal variations in the ice-flow dynamics of Findelengletscher to be a result of changes in the subglacial drainage system (Reference Iken and BindschadlerIken and Bindschadler, 1986; Reference Iken and TrufferIken and Truffer, 1997). The horizontal ice-flow velocity was measured by triangulation at several points on the glacier tongue and revealed a counter-intuitive behaviour. Peak velocity in May and June was followed by a rapid decrease. Starting in October and over the whole winter, horizontal velocity increased (Reference Iken and TrufferIken and Truffer, 1997). Although these results cannot be directly transferred to the current situation or individual years, using the half-year fraction f = 0.5 of v _z to cover the time-span from October 2009 to April 2010 is inappropriate. Instead, the optimal fraction of f = 0.795 is found iteratively by minimizing the rmse of snow depth, d, to the available in situ data from snow probings. The GPR data were not chosen as a reference here because they do not cover the glacier representatively. Nonetheless, the optimal f would correspond well if it were instead obtained with the GPRderived snow depth (f = 0.676; Fig. 13a).

Fig. 13. (a) Determination of the optimal scaling factor, f, of the annual vertical surface velocity by minimizing the rmse of lidar-derived snow depth to GPR and in situ probings. (b) Elevation-dependent bias in lidar- to GPR-derived snow depth (red) and probings (blue) for different choices of the scaling factor, f, shown as mean deviation from n measurements within 100 m elevation bands.

Comparison of lidar-derived snow depth with GPR and probings for different values of f shows the strong elevation-dependent bias in the uncorrected data (f = 0; Fig. 13b). With the correction scaled by f = 0.795 to best fit the snow probings, the imbalance of snow depth overestimation in the ablation area and underestimation in the accumulation area is compensated. Underestimation of snow depth on the lower glacier tongue and in the upper accumulation area might be explained by a local overestimation in the correction for the vertical surface velocity.

The high value of f compared with the expected half-year fraction cannot be explained by seasonal variations in ice-flow velocity alone. Since the correction is based on average mass-balance and geometry changes from a 5 year period, short-term fluctuations are not adequately represented. These can be due to changes in ice-flow velocity but also to variations in the mass balance in recent years that affect firn compaction rates in the accumulation area. The glaciological time series shows a mass balance of +0.83 m w.e. for the accumulation area of Findelengletscher in 2009. For the other years that were used for the calculation of the correction (2005–10), values range from +0.58 to +0.65 m w.e. Thus, firn compaction played a larger role in the period October 2009 to April 2010 than is estimated in the correction, leading to a higher value of f when the correction is fitted to the snow probings.

By calibrating the scaling factor, f, with the snow depth from probings, the correction of the lidar-derived snow depth takes the interannual and seasonal variations in emergence velocity and firn compaction into account. However, the two processes are independent of each other, do not affect the same area and probably do not vary synchronously. Therefore, the scaling factor, f, represents a tuning parameter that integrates all influential processes and allows reconciling the correction of raw lidar-derived snow depth with the individual conditions during the time-span from October 2009 to April 2010.

Conclusions

We have discussed measurements of snow accumulation distribution in April 2010 on Findelengletscher from 403 snow probings, 12.7km of GPR profiles and differencing of two lidar DEMs. The lidar-derived elevation change must be corrected for the glacier emergence velocity and firn compaction. We successfully obtained this correction from long-term mass-balance and geometry changes. All three datasets reveal the general increase in snow depth with elevation. In contrast to the snow probings, GPR and lidar-differencing are able to capture its small-scale variability. Over the entire glacier, snow depth can only partially be explained by elevation, slope and curvature, since deposition and redistribution of snow is strongly determined by local wind fields. We have shown that the explanatory power of these topographic variables is significant within subregions on the glacier. Comparison between snow probings and GPR-derived snow depth indicates a bias in the dataset of snow probings that probably originates from a non-representative distribution of sampling sites.

Single-offset GPR is an effective tool for snow depth measurements that is not affected by glacier dynamics and firn compaction in the accumulation area. It is prone to few uncertainties if estimates for snow density and liquid water content exist. However, the GPR profiles should cover the glacier representatively (e.g. on a grid). Airborne data acquisition is efficient and flexible, but the Swiss telecommunication regulator (OFCOM) allows airborne use of the GPR system only under strict conditions and in uninhabited regions, such as on alpine glaciers. Thus, further technical and regulatory development is needed to fulfil the requirements of national and international authorities for future studies. Lidar provides a fully distributed set of snow depth data with very high spatial resolution, which is not achievable by the application of GPR or conventional glaciological measurements. When corrected for glacier dynamics, lidar-differencing is a powerful tool to assess snow depth distribution on alpine glaciers, but it is expensive in terms of external data requirements. In the presented comparison, the traditional approach of using snow probings remains the only method that is independent of separate input data. However, snow probings typically do not reproduce the small-scale variability of snow depth. We conclude that the choice between the presented methods should depend on the availability of in situ measurements and the desired balance between accuracy and spatially distributed data. GPR provides a compromise in both respects. It requires snow density and wetness estimations that need not necessarily be gained from in situ measurements but that determine the main uncertainties. At the same time, because it provides snow depth along linear profiles, its ability to reproduce the spatial variability lies between the probings and the fully distributed lidar-differencing scheme.

The time-matched availability of the three datasets accounts for the need to have separate data sources as a basis to evaluate the performance of individual methods. Their comparison and combination illustrates a detailed image of winter snow accumulation that benefits from each method’s advantages and will provide a valuable dataset in the future to validate accumulation models of varying complexities.