Spatio-temporal variability in elevation changes of two high-Arctic valley glaciers

Abstract

Uncertainties in estimates of glacier and ice-cap contribution to sea-level rise exist in part due to poor quantification of mass-balance errors, particularly those resulting from extrapolation of sparse measurements. Centre-line data are often assumed to be representative of the glacier as a whole, with little attention paid to extrapolation errors or their effect on mass-balance estimates. Here we present detailed digital elevation model (DEM) measurements of glacier-wide elevation changes over the last ~40 years at two glaciers on Svalbard, Norwegian Arctic. Austre Br0ggerbreen and Midtre Lovenbreen are shown to have lost 27.54 ± 0.98 and 9.65 ± 0.76 × 10⁷m³ of ice, respectively, between 1966 and 2005, findings that we relate to trends in average summer air temperatures and winter accumulation. These volume losses correspond to geodetic balances of -0.58 ± 0.03 and -0.41 ± 0.03 mw.e. a^-1, respectively. Our analysis revealed high spatial complexity in patterns of elevation change, varying between glaciers, between measurement intervals and within and between elevation bins. Balances from extrapolated centre-line geodetic data were the same (within errors) as those from full-coverage DEM differencing in the majority of comparisons, yet significantly underestimated balance in three instances. Additionally, field mass balance from centre-line ablation stake data underestimated balances from full-coverage geodetic measurements during three of six measurement periods. These findings may support the hypothesis that field measurements underestimate Svalbard glacier mass loss, at least partly as a result of the failure of centre-line measurements to account for glacier-wide variations in ablation. Our results demonstrate the importance of deriving accurate interpolation functions and constraining extrapolation errors from sparse measurements.

1. Introduction

Measurement of fluctuations in the mass of the Earth’s glaciers, ice caps and ice sheets is important due to their dominant role in eustatic global sea-level rise. Glacier and ice-cap (GIC;ice bodies outside the Greenland and Antarctic ice sheets) contribution to sea-level rise has accelerated over the past decade (Meier and others, 2007), currently comprising around half of the present-day sea-level budget (Cazenave and others, 2009). However, the errors associated with estimates of this contribution are large, in part due to uncertainty of existing mass change measurements. Improvement of mass-balance estimates and better understanding of their associated errors will help to reduce these uncertainties.

GIC mass balance is usually measured by summing accumulation and ablation terms from stakes distributed over the ice surface (known as the glaciological method), or by measuring changes in ice surface elevation over time, often extrapolated from repeat profiles (known as the geodetic method). Hereafter we refer to mass balance from the glaciological method as ‘net mass balance’ or ‘field mass balance’, and mass balance from the geodetic method as ‘geodetic balance’. Both these approaches suffer from errors resulting from inadequate size and spatial coverage of available datasets, and from challenges related to upscaling these sparse measurements. At the individual glacier scale, constraints imposed by working costs, logistics and access safety and the need for standardized measuring procedures has led to a glaciological sampling strategy whereby stakes are usually placed in longitudinal transects along a central ice flowline (e.g. Østrem and Brugman, 1991). Net mass balance, B _n, is then derived by assigning point (stake) measurements of specific balance, b _n, to map contours of equal balance (usually based on elevation bands), and integrating over the total area of the glacier (Østrem and Brugman, 1991;Paterson, 1994). For similar reasons, estimates derived from the geodetic method are also often extrapolated from centre-line or sparse profile measurements (e.g. Sapiano and others, 1998;Arendt and others, 2002;Abdalati and others, 2004). Typically, these measurements are averaged within an elevation bin, and that average is weighted by the area of the bin and integrated (by summation) over the total glacier area.

Both the glaciological method and the profile-based geodetic method rely on the fact that for most alpine-type valley glaciers altitude is the dominant control on accumulation and ablation (e.g. Nye, 1960;Paterson, 1994). Restricting data collection to centre-line points or profiles provides no information on variability within individual elevation bands, and by definition assumes transverse variations to be of secondary importance (Fountain and Vecchia, 1999). One result of this assumption is that nonelevation-related variations in accumulation (due to factors such as avalanches, snowfall or wind redistribution) or ablation (due to factors such as topographic shading or debris material reducing surface albedo) may not be included in mass-balance estimates. A further concern is that measurement uncertainty cannot be accurately quantified and may accumulate over time (Elsberg and others, 2001). Given that glaciers evolve towards an equilibrium condition of zero balance, error assessment may be fundamental to determining the sign, as well as the magnitude, of balance (Fountain and Vecchia, 1999).

The error budget of geodetic balance estimates from centre-line data commonly focuses on profile and map errors, particularly the quality of digital elevation model (DEM) or contour plots used to extrapolate elevation changes, rather than any uncertainty associated with the use of profile data to represent all changes within an elevation band (e.g. Rignot and others, 2003;Abdalati and others, 2004). Usually this is simply due to a lack of data with which to quantify extrapolation uncertainty. Exceptions include studies with access to validation data, such as full- coverage DEMs created independently from the primary geodetic datasets (e.g. Sapiano and others, 1998;Arendt and others, 2006). These studies calculate extrapolation or ‘profile-to-glacier’ errors as the standard deviation of elevation differences within elevation bins, averaged across all bins across the glacier. In the absence of full-coverage independent validation data, it is impossible to accurately quantify profile-to-glacier errors or their impact on geodetic balance estimates. With the exception of Berthier and others (2010), few studies to date have provided both a detailed picture of glacier surface elevation-change variability and a quantification of the difference between geodetic balance estimates derived from both centre-line extrapolation and full-coverage DEM differencing.

In this paper, we combine airborne remote-sensing datasets following the methods of James and others (2006) and Barrand and others (2009) to construct high-quality, full- coverage surface elevation models of two Svalbard valley glaciers. These models are then differenced to examine spatio-temporal variability of glacier elevation changes over the past four decades. We calculate 39years of geodetic balance at the two sites, deriving estimates using both glacierwide (full-coverage) and down-sampled (sparse centre-line) profile datasets. A comparison of these two approaches allows us to quantify the difference between centre-line extrapolated and full-coverage geodetic balances.

2. Study Area

The sites chosen for this analysis were Austre Brøggerbreen (AB) and Midtre Lovenbreen (ML), two northward-flowing valley glaciers located on the Brogger peninsula, northwestern Spitsbergen, Svalbard (Fig. 1). For the most recent year of mapping (2005), AB and ML had surface areas of 10.2 and 5.1 km² and central ice flowline lengths of 5.8 and 4.4 km, respectively. Although the glaciers neighbour one another and share both a similar altitudinal range (~60- 650 m a.s.l.) and climatological setting, they differ in thermal regime, with AB being predominantly cold-based and ML being polythermal (Bjornsson and others, 1996; Rippin and others, 2005). The two glaciers also have dissimilar surface geometries, with AB characterized by greater width and perimeter resulting from three substantial tributary basins, and ML more closely approximating a typical valley-glacier configuration of a single central ice-flow unit draining several smaller tributary basins (Fig. 1). Proximity to the settlement of Ny-Ålesund has resulted in long-term (>40year) field mass- balance records for both glaciers (from accumulation maps and centre-line ablation stakes) and meteorological data from a nearby weather station.

Fig. 1. Location map of Austre Brøggerbreen (AB, left) and Midtre Love´nbreen (ML, right), Brøggerhalvøya, northwest Svalbard (inset). Ice surface-elevation contours (2005) and centre-line mass-balance stake locations are shown.

3. Data and Methods

3.1. Topographic datasets

DEMs of AB and ML were generated for the years 1966, 1977, 1990 and 2005 from a combination of airborne laser scanning (lidar) and digital photogrammetry (Table 1). Lidar datasets were collected over both glaciers in late summer 2005 with an Optech ALTM3033 scanning system using a scan rate of 28 Hz, a laser pulse rate of 33 kHz and a scan angle of ±18°, resulting in mean along- and across-track point spacing of 1.38 and 1.33 m, respectively. Lidar GPS baselines were differentially post-processed from base- station data collected via an antenna splitter at the International GPS Service (http://igscb.jpl.nasa.gov/) known point ‘NYA1’ located at the geodetic observatory at Ny-Ålesund airport. Point cloud lidar data were processed into 20 m DEMs using a linear triangulation interpolation routine. Historical DEMs at 20 m resolution were produced by digital photogrammetry using SOCET SET software from archival stereo imagery collected by the Norwegian Polar Institute (NPI) near the end of the balance year in 1966, 1977 and 1990. Ground control information was extracted from 2005 lidar datasets to perform three-dimensional triangulation of historical stereo imagery, thus circumventing the need for field-measured control. James and others (2006) and Barrand and others (2009) provide complete descriptions of this methodology.

Table 1. AB and ML digital elevation data sources and their vertical errors. Photogrammetric DEMs of AB and neighbouring ML were generated from the same photo block in each year (1966, 1977, 1990), hence their identical error statistics. Aerial photographs were scanned from panchromatic (P) or infrared false-colour (IR) negatives. The quality of each bundle adjustment is expressed in terms of the root mean square (RMS) of adjusted ground control point (GCP) positions, and the total image unit weight (a good indicator of the quality of the overall solution). Elevation accuracy is expressed as the standard deviation, σ, of elevation residuals between topographically stable test sites and 2005 lidar models (photogrammetric DEMs), and between differential GPS check data and 2005 elevations (lidar DEMs)

3.2. Surface elevation change

Glacier surface elevation changes were derived by differencing DEMs for the periods 1966-77, 1977-90 and 19902005. Following the naming convention of symbols and related measurements described by Arendt and others (2006), individual surface elevations, h were differenced to obtain elevation change, Δh. Annual elevation-change rates, , were calculated by dividing Δh by the time between DEMs. All elevation- and volume-change rates are therefore time-averaged over the measurement period, rather than representing instantaneous rates of change.

3.3. Volume change and geodetic balance

Total volume change, B (m³), was obtained by pixel summation of difference DEMs (e.g. Etzelmuller and others, 1993), whereby the i pixels of each difference DEM, Δh_i , containned within the larger glacier surface, A, were summed and multiplied by the area of each pixel, l _p ², where l _p is the grid spacing, expressed as

(1)

The geodetic balance rate, Ḃ (m³w.e. a^-1), was determined by multiplying B by 0.918 (the ratio of the density of ice to water, p _i/p _w; Paterson, 1994) and dividing by the time between epochs. This approach assumes Sorge’s law (Bader, 1954), that the density-depth profile at the initial (t ₁) and final (t ₂) time was the same over the entire glacier. Given that the datasets were collected over predominantly snow- free terrain, we consider this a reasonable assumption. To compare glaciers of different sizes, Ḃ was divided by the average of the old and new glacier areas to give the area-averaged net geodetic balance, (mw.e. a^-1) (e.g. Echel-meyer and others, 1996;Arendt and others, 2002).

To examine the ability of sparse centre-line geodetic measurements to estimate glacier-wide geodetic balance, we down-sampled full-coverage DEM difference models to centre-line-point and centre-line-profile only datasets. Three approaches were therefore used to estimate geodetic balance: (1) full-coverage difference models, (2) extrapolated centre-line points and (3) extrapolated centre-line profiles;hereafter referred to as ‘measured’ (1) and ‘extrapolated’ (2 and 3) approaches. Point and profile elevation changes were extracted from AB and ML difference model centre-line transects following the location of mass-balance stake lines (Fig. 1). Geodetic balances derived from extrapolated centre-line points (2) were calculated by extracting a single Δh measurement from the intersection of the centre line and each 50m bin median elevation (e.g. 225 m for the 200-250 bin). These values were extrapolated across the glacier by area-averaging (i.e. applied to all pixels of each respective bin). Geodetic balances were then calculated as above. Averaged-profile geodetic balances (3) were determined by calculating the mean Δh of all difference DEM pixels along the centre line of each 50 m elevation bin and applying this value to each pixel of each respective bin.

Full-coverage DEM differencing was considered fully representative of all elevation changes, given that every Δh value from the area of each bin was included in the calculation. Experiments using point and profile extrapolations were designed to simulate the extrapolation of sparse geodetic data, similar to those commonly available from airborne or satellite altimetric profiling instruments (e.g. Abdalati and others, 2004;Bamber and others, 2005; Sauber and others, 2005).

3.4. Error analysis

Quantification of photogrammetric DEM error can be difficult when modelling a dynamic surface. DEM error depends on results of both the photogrammetric bundle adjustment and the DEM generation procedure (image correlation). The former is affected by factors including the quality of ground control points (GCPs), their measurement and the camera calibration. The latter is affected by factors including image scale, image quality and surface texture. Thus, overall model solution and DEM quality are typically expressed in terms of both bundle adjustment and stereomatching results. The quality of the bundle adjustment is reported in the form of the root-mean-square error (RMSE) of GCP residuals (in x,y,z components) and image coordinate residuals (Table 1). While these statistics confirm the quality of the adjustment, they may not be indicative of DEM quality which also depends on the quality of photographs and the results of the stereo-matching algorithm and manual editing. DEM quality is therefore usually assessed by comparison to some independent check data that are regarded as ‘truth’.

We quantified lidar elevation accuracy as the standard deviation, σ, of elevation residuals between DEM surfaces and ground-based differential GPS (DGPS) check data collected on the surface of AB (n = 1121) and ML (n = 529), on days either side of airborne surveys (Barrand and others, 2009). This approach, however, was not feasible for historical photogrammetric DEMs given the absence of time-coincident ice surface check data. Instead we estimated errors by selecting a 1.64km² topographically stable test site of comparable surface texture and relief within the forefield between both glaciers (and within the stereo coverage of each historical image pair). This test site contained stable terrain (bedrock, vegetated surfaces) of similar surface texture and relief to the glacier surface and was therefore a more appropriate DEM comparison site than the active forefield (subject to thermo-erosion, fluvial reworking and mass movements) and the very steep surrounding mountainsides (where small positional errors translate to large vertical errors). We describe the error of historical photogrammetric DEMs by deriving elevation- change histograms and descriptive statistics (mean and standard deviation) between each historical surface and elevations from the 2005 lidar DEM (which we assume to be the ‘truth’ model) (Fig. 2).

Fig. 2. Frequency distribution and statistics of DEM difference models calculated between AB and ML historical photogrammetry (1966, 1977, 1990) and contemporary lidar (2005), over a topographically stable forefield test site.

We calculated the errors of measured and extrapolated volume changes and geodetic balances as follows. The error of an individual point elevation change, E _PT, was estimated as the root sum of squares (RSS) of DEM error

(2)

where σ_DEM are uncertainties calculated for individual DEMs from analysis of non-glacier terrain (Table 1). The error of the mean elevation change of the entire glacier for measured (non-extrapolated) DEMs, E, assuming errors to be random and normally distributed over the glacier surface, was then calculated by

(3)

where n is the number of independent measurements (or pixels).

Given that continuous surface DEMs have some degree of spatial autocorrelation (Etzelmuller and others, 1993; Rolstad and others, 2009), we calculated n assuming a correlation scale of 1000 m (e.g. Nuth and others, 2007). Rolstad and others (2009) used geostatistical methods to estimate three scales of spatial correlation in photogram- metric difference DEMs of the western Svartisen ice cap, Norway: hundreds of metres (260 and 430m);kilometres (3100m);and tens of kilometres (17 000 m). The authors considered the smallest scale of correlation to result from matching errors, and the larger scales of correlation to result from errors in photogrammetric block instability and absolute orientation. Given the quality of our photogrammetric adjustments (Table 1), we consider 1000 m to be a conservative estimate of correlation scale in our difference DEMs. We calculated volume-change errors, E _VOL, by multiplying E by the area, A, of each glacier (assuming any errors in A to be insignificant at the spatial scale of our volume changes).

In the case of extrapolated DEM measurements, we calculated volume-change and mass-balance errors following the approach of Nuth and others (2010). We calculated profile-to-glacier (extrapolation) errors, E _EXT, as the standard deviation of elevation change, Δh, within 50 m bins (e.g. Berthier and others, 2004; Arendt and others, 2006; Nuth and others, 2010). Mean elevation-change errors, E, were calculated by combining point elevation errors, E _PT, and extrapolation errors, E _EXT. As elevation changes are averaged by elevation, errors are reduced by the square root of the number of independent measurements within each bin, assuming errors to be spatially correlated at a scale of 1000 m (as above). As both these error sources were assumed to be random and normally distributed, they are combined to produce total elevation-change errors for each bin, Z, by

(4)

Volume-change errors, E _VOL, were then estimated by the RSS of all elevation bin errors, E _Z, multiplied by bin area, A _Z, assuming that errors are independent between bins (e.g. Nuth and others, 2010),

(5)

4. Results

4.1. Surface elevation change, volume change and geodetic balance

Difference DEM maps of AB and ML during the periods 1966-77, 1977-90 and 1990-2005 show widespread surface lowering indicative of glacier mass wastage, with only small areas of positive elevation change at the highest elevations (Figs 3 and 4). At AB, positive Δh of 0-10 m was evident at all the highest elevations between 1966 and 1977 (Fig. 3a), and at some lower elevations between 1977 and 1990 (Fig. 3b), and yet was almost completely absent by 1990-2005 (Fig. 3c). At ML, the largest areas of positive Δh occurred between 1966 and 1977 (Fig. 4a). These areas were slightly reduced in size by 1977-90 (Fig. 4b), and by 1990-2005 were restricted to just one location on the southwestern cirque (Fig. 4c). Positive Ah in this cirque resulted from a small patch of seasonal snow cover likely causing an elevation blunder in the 1990 DEM. This positive elevation change is masked (Δh = 0) for the purposes of further volume-change and geodetic-balance calculations. As expected, all Δh maps showed a broad trend to greater thinning at lower elevations, with less thinning towards progressively higher elevations. This trend, however, was noticeably non-uniform, varying in both space and time. At AB, the region of greatest thinning between 1966 and 1977 was the northeast part of the glacier snout within the lowest elevation band (50-100 m a.s.l.). By 1990-2005, a large section on the west side of the main glacier trunk between 100 and 200 m had thinned the most. This example illustrates the difficulty of assuming that the centre line is representative of elevation-band-wide Δh. In this case, a single point or averaged profile from the centre of the band would likely underestimate band-wide mean Δh.

Fig. 3. Surface elevation change, Δh (m), at AB, 1966–77 (a), 1977–90 (b) and 1990–2005 (c). Plots are projected in Universal Transverse Mercator (UTM) World Geodetic System 1984 (WGS84) (m), zone 33 north. White gaps within the glacier margin represent nunatak locations.

Fig. 4. Surface elevation change, Δh (m), at ML, 1966–77 (a), 1977–90 (b) and 1990–2005 (c). Plots are projected in UTM WGS84 (m), zone 33 north. The semi-transparent box in the southwestern cirque (b, c) delineates the region of blunder masking in the 1990 DEM.

Comparison of Figures 3 and 4 shows variations in spatial patterns of thinning between the two glaciers during all three measurement periods. The regions of greatest thinning at ML were confined predominantly to the lower glacier tongue below 200 m a.s.l. In contrast, AB thinned by similar amounts over a larger area, and up to the 250 m contour (particularly in the most recent period; Fig. 3c). We examined surface elevation-change variability in space and time at AB and ML by plotting elevation-change rates at each point in the DEM, Δh, as a function of elevation (Fig. 5). Between 1966 and 1977, AB thinned ~0.5 m a^-1 more than ML within the elevation range 150-400 m. Above 400 m, the distribution of Δh is similar between glaciers. Comparison between periods reveals that during the most recent measurement interval (1990-2005), Δh at both glaciers was predominantly negative, even at the highest elevations.

Fig. 5. Elevation-change rates, Δh, as a function of surface elevation at AB (top, blue) and ML (bottom, red).

The total volume loss of AB between 1966 and 2005 was approximately three times that of ML within errors (Table 2). Respective rates of geodetic balance were -0.58 ± 0.03 and -0.41 ± 0.03 mw.e. a^-1 for AB and ML over the entire measurement period. Within individual measurement periods, AB and ML show a similar pattern of geodetic balance, with negative rates between 1966 and 1977, slightly less negative balance between 1977 and 1990, followed by the most negative balances during the most recent period, 1990-2005 (Table 2). Rates of geodetic balance were greater (more negative) at AB than at ML during all measurement periods.

Table 2. Total volume change, B, and (area-averaged) geodetic mass balance, , of AB and ML at measurement periods between 1966 and 2005, from lidar and lidar-controlled photogrammetric DEM differencing

4.2. Comparison of extrapolated and measured geodetic balances

Comparison of point- and profile-extrapolated geodetic balances with measured rates (from full-coverage DEM differencing) shows that for 9 of 12 comparisons extrapolated balances were the same as measured balances within stated error bounds (Table 3). At both AB and ML during the periods 1966-77 and 1977-90, point- and profile- extrapolated balances were not significantly different from measured balances (i.e. outside the range of stated errors) (Table 3). Likewise, between 1990 and 2005 at AB, profile- extrapolated geodetic balance was the same (within errors) as that from full-coverage geodetic measurements. We may therefore state that, within the range of uncertainty of our measurements, in these nine comparisons geodetic balance from extrapolation of centre-line point and averaged-profile measurements was no different than that from full-coverage (no extrapolation) geodetic measurements. During the period 1990-2005, however, the point extrapolation approach at AB and the point and profile extrapolation approaches at ML significantly underestimated measured geodetic balance. In these three instances, centre-line extrapolated balance underestimated measured balance by 12% (AB 1990-2005; point), 22% (ML 1990-2005;point) and 22% (ML 1990-2005;profile) (Table 3).

Table 3. Profile-to-glacier errors, E _EXT, and geodetic balance, of AB and MB at measurement periods between 1966 and 2005, derived from point extrapolation, profile extrapolation and fullcoverage DEM differencing. Extrapolated geodetic balances that differ significantly from full-coverage balances are marked in bold. Errors of point- and profile-extrapolated approaches differ from fullcoverage DEM differencing as they include a parameterization of extrapolation errors, yet are identical as they each have the same number of independent measurements per bin (n = 1) at our chosen scale of spatial correlation

Profile-to-glacier (extrapolation) errors were similar between AB and ML, despite differences in the size and geometry of the two glaciers (Table 3). At both glaciers, extrapolation errors increased in subsequent measurement periods, a finding that may be broadly related to observed differences between extrapolated and measured geodetic balances. The three instances where extrapolated balances significantly underestimated measured balances (point extrapolation at AB between 1990 and 2005, and point and profile extrapolation at ML between 1990 and 2005) were coincident with the largest profile-to-glacier errors (Table 3). Both measurement intervals with the largest variance in Δh therefore coincide with significant differences between extrapolated and measured geodetic balances.

5. Discussion

We have shown that despite the similarities between AB and ML, elevation changes at the two glaciers were variable in both space and time. Geodetic balances were negative at both glaciers between 1966 and 2005, with a period of slightly less negative balance between 1977 and 1990, and the most negative balances at ML occurred during the most recent period, 1990-2005. Experiments designed to simulate the calculation of geodetic balances from extrapolation of sparse centre-line data showed that balances from extrapolated data were the same (within errors) as balance rates derived from full-coverage measurements in the majority of instances (9 of 12). However, balance estimates from extrapolated geodetic datasets underestimated measured rates by up to 22% in three cases. We now compare these results with previously published geodetic balances and with field-measured mass balances covering the same time period. Finally, we place our results in the context of local meteorological conditions.

5.1. Comparison with previous geodetic change estimates

An annual balance of-1.06 ± 0.30 m w.e. a^-1 was calculated between 1970 and 1990 from geodetic data covering the lower ~3 km of AB (Pope and others, 2007). The equivalent geodetic balance rate for the same time period from our data was -0.55 ± 0.05 m w.e. a^-1. The discrepancy between these two values is likely to be a result of patchy data and overrepresentation of lower-elevation thinning in the estimate of Pope and others (2007). Their estimate was based on DEMs with no manual editing where poor image or scanning quality resulted in correlation failures in areas of shadow and at high elevations.

Previous studies at ML have estimated a geodetic balance of -0.61 m w.e. a^-1 between 1977 and 1995 (Rippin and others, 2003), and an acceleration of glacier thinning (equivalent to geodetic balances of -0.14, -0.18, -0.28, -0.30, -0.46 and -0.62 m w.e. a^-1 between the years 1936, 1962, 1969, 1977, 1995, 2003 and 2005) from comparison of a range of geodetic data (Kohler and others, 2007). Our geodetic balance rate between 1977 and 1990 was -0.34 ± 0.04 m w.e. a^-1, around half the 1977-95 estimate of Rippin and others (2003). The large measurement uncertainty of the Rippin and others (2003) estimate (-0.61 ± 0.7 m w.e. a^-1) implies that their input DEM datasets were of insufficient quality to measure geodetic balance at this glacier. Our longer measurement period and high- accuracy, full-coverage datasets provide a strong alternative to Rippin and others (2003) for DEM-based (geodetic) mass balance of ML. While direct comparisons with the results of Kohler and others (2007) are difficult to make due to their different period lengths and lack of error bounds, our findings are similar (within errors) between 1977 and 2005. Although we do not identify constant acceleration in thinning/mass loss (geodetic balances were slightly less between 1977 and 1990 than between 1966 and 1977, although within errors), our results confirm the finding of accelerating mass loss, with the most recent period (19902005) experiencing the most negative rate of balance.

5.2. Comparison with field-measured mass balance

We plotted complete time series of field-measured summer, winter and net area-averaged surface mass balances at AB and ML between 1967 and 2006 from data collected by the NPI. Polynomial-fit trend lines indicate a consistently negative net balance of ~-0.4ma^-1 between 1967 and the mid-1990s at both glaciers, with a trend toward more negative balance since then (Fig. 6). Although net balances were generally more negative at AB than at ML, this trend was only very pronounced in a few years, including 198990, the record low year of 1998 and 2004-06. In contrast, our geodetic results suggest that balances at AB were consistently more negative than at ML (Table 2). This discrepancy may be a result of the failure of centre-line stake measurements to account for ice losses within the larger area and more complicated geometry of AB.

Fig. 6. Field-measured mass balance of AB (blue) and ML (red), 1968–2006. Separate plots are shown for summer (square symbols), winter (triangles) and net (dots) balances.

To directly compare field and geodetic measurements, we calculated cumulative mass balances within each measurement period (1966-77, 1977-90 and 1990-2005) (e.g. Elsberg and others, 2001). Field measurements at AB and ML began in 1967 and 1968 respectively, a fact we account for in our geodetic comparisons by cumulating 1966-77 balances from 1967 and 1968 onwards, respectively. We accept that this approach may introduce some bias in the early-period comparisons, the level of which is dependent on the state of balance at AB in 1966 and at ML in 1966 and 1967. While we were able to compare full-coverage geodetic measurements at ML with field measurements (as data from ablation stakes are applied to a hypsometry of the entire glacier surface), at AB we limited our geodetic measurements to the eastern wing only, for consistency with the 6.12 km² hypsometry applied by the NPI to their centre-line stake measurements (personal communication from J. Kohler, 2010). For comparison purposes therefore, the AB geodetic balances in this subsection refer only to the hypsometry of the eastern wing of the glacier. AB geodetic balances presented elsewhere in this paper refer to the entire geometry of the glacier.

At AB, cumulative mass balances from field measurements during the periods 1966-77, 1977-90 and 1990-2005 were -4.98, -5.77 and -9.09 m w.e., respectively. The pattern is similar to our results from geodetic methods, yet these were more negative for all three periods (significant during 1966-77 and 1990-2005; respectively -6.45 ± 0.33, -5.95 ± 0.42 and -9.50 ± 0.22 m w.e.). At ML, cumulative mass balances from field measurements during the periods 1986-77, 1977-90 and 1990-2005 were-3.59, -4.59 and -7.45 m w.e., respectively. Cumulative balances from geodetic methods during the first two periods were the same (within errors) as field measurements (-3.70 ± 0.40 m w.e. for 1967-77, and -4.76 ± 0.56 m w.e. for 1977-90), yet were significantly more negative than field measurements during the most recent period (1990-2005;-8.22 ± 0.28 m w.e.). While it is possible that the discrepancy between these values may be explained in part by the uncertainty of field measurements, no such formal error analysis of these data has yet been undertaken (personal communication from J. Kohler, 2010).

A number of previous studies on Svalbard have suggested that field measurements may underestimate mass loss from glaciers (Fleming and others, 1997; Hagen and others, 2000;Rippin and others, 2003), a hypothesis that is partially supported by our results. Possible explanations for this underestimation may include stakes melting into the firn resulting in overestimation of accumulation (Braithwaite, 2002) or higher net ablation towards the margins (away from the centre line) due to advection of turbulent energy from valley sides (Rippin and others, 2003). Given that centre-line extrapolated geodetic data were either the same within errors or underestimated full-coverage geodetic balances (Table 3), we contend that underestimation of field-measured mass balance is at least partly a result of glacier-wide spatial complexity in patterns of ablation not accounted for by measurements at the centre line (e.g. Arnold and others, 2006).

5.3. Drivers of mass-balance trends

Previously published work has identified strong statistical relationships between the net (field-measured) mass balance of AB (R ² = 0.90) and ML (R ² = 0.83), and combined time series of climate-station variables (mean summer air temperatures and summed winter precipitation) (Lefauconnier and Hagen, 1990; Lefauconnier and others, 1999). We therefore view our results in the context of mean summer temperature and summed winter precipitation, from climate-station data at nearby Ny-Ålesund.

Our results at AB and ML coincide with a period of recent summer (June-August) warming yet no clear increase in winter (October-May) precipitation (Fig. 7). Between 1966 and 1977, average summer air temperatures were 3.8°C, during 1977-90 they fell slightly to 3.6°C, and between 1990 and 2005 they rose to 4.1 °C. This trend is consistent with our geodetic balance results of slightly less negative balance during 1977-90 (compared to 1966-77), and more negative balance during 1990-2005 (Table 3). Less negative geodetic balances between 1977 and 1990 may be at least partly explained by reduced summer temperatures in 1977, 1982 and 1987, accompanied by two years of near-zero balance in 1977 and 1982, and positive net balance in 1987 (Fig. 6). Increasing rates of mass loss at AB and ML may be a result of not just warmer summer (June-August) temperatures causing additional melt, but also decreases in winter accumulation. Although winter precipitation remained relatively stable between 1966 and 2005 (Fig. 7), winter balances at both glaciers have been decreasing since the early 1990s (Fig. 6).

Fig. 7. Mean summer (June–August) temperature (a) and summed winter (October–May) precipitation (b) data from the climate station at Ny-Ålesund. Thick black curves represent 5 year moving averages. Bars in (a) show geodetic data measurement periods for AB and ML. Data courtesy of the Norwegian Meteorological Institute.

6. Summary and Concluding Remarks

We have presented an analysis of spatio-temporal variability in ice surface elevation changes over the last ~40 years at two small valley glaciers in northwest Svalbard. Our data show widespread surface lowering with complex spatial patterns of elevation change, varying within and between elevation bins, between measurement periods and between neighbouring glaciers. Our long measurement period and accurate, full- coverage datasets provide a benchmark for geodetic mass balance of Austre Broggerbreen and Midtre Lovenbreen, two high-Arctic valley glaciers which were shown to have lost 27.54 ± 0.98 and 9.65 ± 0.76 × 10⁷ m³ of ice between 1966 and 2005. Experiments designed to simulate extrapolation of sparse geodetic data showed that balance estimates from centre-line measurements were typically the same (within errors) as those from full-coverage measurements, yet in several instances underestimated full-coverage measurements by as much as 20%. Likewise, during three of six measurement periods, mass balance from field measurements underestimated values from full-coverage geodetic methods. These results may support the hypothesis that field measurements underestimate Svalbard glacier mass loss, most likely due to complex variations in ablation not accounted for by measurements at the centre line.