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. Reference James, Murray, Barrand and BarrJames and others (2006) and Reference Barrand, Murray, James, Barr and MillsBarrand and others (2009) provide complete descriptions of this methodology.

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 Reference ArendtArendt 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. Reference Etzelmuller, Vatne, Odegard and SollidEtzelmuller 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

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; Reference PatersonPaterson, 1994) and dividing by the time between epochs. This approach assumes Sorge’s law (Reference BaderBader, 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. Reference EchelmeyerEchel-meyer and others, 1996;Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt 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. Reference AbdalatiAbdalati and others, 2004;Reference Bamber, Krabill, Raper, Dowdeswell and OerlemansBamber and others, 2005; Reference Sauber, Molnia, Carabajal, Luthcke and MuskettSauber 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 (Reference Barrand, Murray, James, Barr and MillsBarrand 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

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

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

Given that continuous surface DEMs have some degree of spatial autocorrelation (Reference Etzelmuller, Vatne, Odegard and SollidEtzelmuller and others, 1993; Reference Rolstad, Haug and DenbyRolstad and others, 2009), we calculated n assuming a correlation scale of 1000 m (e.g. Reference Nuth, Kohler, Aas, Brandt and HagenNuth and others, 2007). Reference Rolstad, Haug and DenbyRolstad 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 Reference Nuth, Moholdt, Kohler, Hagen and KaabNuth 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. Reference Berthier, Arnaud, Baratoux, Vincent and RemyBerthier and others, 2004; Reference ArendtArendt and others, 2006; Reference Nuth, Moholdt, Kohler, Hagen and KaabNuth 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

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. Reference Nuth, Moholdt, Kohler, Hagen and KaabNuth and others, 2010),

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.

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).

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.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 (Reference Pope, Murray and LuckmanPope 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 Reference Pope, Murray and LuckmanPope 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 (Reference RippinRippin 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 (Reference KohlerKohler 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 Reference RippinRippin and others (2003). The large measurement uncertainty of the Reference RippinRippin 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 Reference RippinRippin and others (2003) for DEM-based (geodetic) mass balance of ML. While direct comparisons with the results of Reference KohlerKohler 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. Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg 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 (Reference Fleming, Dowdeswell and OerlemansFleming and others, 1997; Reference Hagen, Etzelmuller and NuttallHagen and others, 2000;Reference RippinRippin 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 (Reference BraithwaiteBraithwaite, 2002) or higher net ablation towards the margins (away from the centre line) due to advection of turbulent energy from valley sides (Reference RippinRippin 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. Reference Arnold, Rees, Hodson and KohlerArnold 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) (Reference Lefauconnier and HagenLefauconnier and Hagen, 1990; Reference Lefauconnier, Hagen, Örbæk, Melvold and IsakssonLefauconnier 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.