The studied glaciers include Bertilbreen (BLB), Svenbreen (SVB), Vestre Muninbreen (VMB), Austre Muninbreen (AMB), Ferdinandbreen (FDB), Elsabreen (ELB) and an unnamed glacier known as 16 819 in the Svalbard database (Hagen et al. 1993) and called No Namebreen (NNB) in this study (Fig. 1). All are situated in the Billefjorden and Austfjorden basins of Dickson Land, some kilometres north of the abandoned Russian mining town of Pyramiden. The mean annual temperature at Svalbard Airport, 50 km south of Pyramiden, is −5.1°C, and the mean summer (June–August) temperature is 4.9°C, as measured by the Norwegian Meteorological Institute in the period 1981–2010. Meteorological investigations in the study area have shown that summers are similarly warm to those at Svalbard Airport or even warmer (Rachlewicz 2003, 2009a b; Przybylak et al. 2006; Przybylak et al. 2007; Rachlewicz & Styszynska 2007; Láska et al. 2012), and according to Korjakin et al. (1985), it is the warmest region in Spitsbergen in the summer months. Previous investigations of glaciers in this region were concerned with their geomorphological impact on the landscape (e.g., Karczewski 1989; Kostrzewski et al. 1989; Gibas et al. 2005; Rachlewicz 2009a, b; Szuman & Kasprzak 2010; Ewertowski et al. 2010Ewertowski et al. 2012; Evans et al. 2012), although Soviet scientists performed detailed glaciological work on BLB in the 1970s and 1980s (e.g., Gohman et al. 1982; Žuravlev et al. 1983; Gokhman & Khodakov 1986; Gokhman 1987; Troicki 1988).

The studied glaciers represent different sizes, types, aspects and elevation range (Table 1, Fig. 1). All of them have been continuously retreating since the termination of Little Ice Age, at sometime around 1920 (Rachlewicz et al. 2007). Radar soundings on the two largest glaciers, BLB and SVB (both ca. 4 km²), show that they are relatively shallow, with maximum depth not exceeding 150 m, and that they are mostly cold, with limited zones of temperate ice in their deeper sections (Žuravlev et al. 1983; Małecki 2013). VMB and AMB may have similar thermal structure due to their size, while the other ice masses, NNB, FDB and ELB, are most likely entirely cold, as they are relatively small and probably too shallow to have any temperate ice. The thermal regime of the studied glaciers influences their dynamic activity, recognized by Hagen et al. (2005) as a crucial element in understanding elevation changes. Direct measurements performed on BLB and SVB indicate maximum flow velocities on the order of only 2–3 m a⁻¹, (Mavlyudov 2010; Małecki 2013) and none of the studied glaciers are known to surge. Therefore, overall dynamically driven elevation changes are assumed to be small compared to those driven by climate, for the study period (Nuth et al. 2007; James et al. 2012).

The 1960 glacier geometry was taken from a 1:100 000 topographic map by the Norwegian Polar Institute (NPI), constructed from 1:50 000 vertical aerial images (Norwegian Polar Institute 1988). Contour lines of 50 m vertical interval were digitized by and obtained from the NPI. In 1990In 1990a similar NPI photogrammetric survey was carried out over the study area, with images obtained at a scale of 1:15 000. From these photographs, the NPI created a DEM with 20 m resolution. The most recent survey took place in 2009. Digital images taken by the NPI with ground resolution of about 0.5 m were used to generate a new DEM by Strzelecki (2013). However, about 5% of the 1990 area of VMB is missing, in the lateral part close to the front, because of incomplete photographic coverage. Therefore, elevation change in this area is interpolated from the neighbouring surveyed zone to obtain a full-coverage map of the 1990–2009 elevation variations for VMB. DEM resolution is 20 m for most of the studied ice masses, but 8 m resolution for BLB and FDB. To check for potential misalignments of the available data, a universal co-registration correction was used, as described in detail by Nuth & Kääb (2011). The data sets were iteratively aligned to the 1990 DEM until the computed misalignments in x, y and z directions were lower than 0.2 m.

Kääb (2008) and Nuth et al. (2010) demonstrate that even though the vertical resolution of 1960s–1970s maps is relatively low, they may be successfully used for glaciological analysis of large (80–260 km²) ice caps or for comparably sized regions. Here it is shown that geometry changes of individual small glaciers (≤4 km²) located in a mountainous terrain may be assessed as well. Digitized contour lines were exported to a point format, and each vertex was interpolated to the grid of the 1990 DEM by bilinear interpolation. Hence, local elevation differences (dh) between 1960 and 1990 could be obtained only at these points on the glacier surface. To extrapolate dh from their point locations to whole glaciers, three different procedures were used, which is somewhat similar to those described in Kääb (2008).

In Method A, a map of surface elevation changes of 20 m resolution was created using estimated dh at each point. Different interpolation schemes were tested: kriging, inverse distance weighting and the Topo to Raster tool (a discretized thin plate spline technique available in ESRI ArcGIS software). In the latter, contours of equal elevation change at the glacier fronts are drawn manually because high gradients of dh and sparse data produced unrealistic patterns, i.e., closed zones of high surface lowering along contours, with areas of modest thinning between them. The Topo to Raster procedure, with manual corrections, gave the best results, showing good agreement between gridded values of the interpolated raster and the measured dh at each point, as well as smooth transition between zones of high and low dh in contrast to the other interpolation schemes. Topo to Raster was therefore chosen as the best approach for interpolation of available data. When compared to results delivered by other interpolation schemes, it gave the most negative values due to manual improvements at the fronts, where other techniques failed to reflect high thinning. Obviously, spatial variability of glacier elevation changes described with this method is somewhat simplified as it lacks the input data from wide areas between the contour lines, but it reproduces the general pattern well.

In Method B, a DEM with 20 m resolution was generated from the 1960 contour lines using a finite difference interpolation technique, to create a hydrologically correct relief model (Hutchinson 1989). The DEM was aligned to the 1990 DEM grid, so the latter could be subtracted from 1960 DEM, resulting in a raster of elevation differences. The method's weak point is the questionable quality of the 1960 DEM, derived from sparse contour lines.

In Method C, glacier elevation changes along 1960 contours were averaged and then multiplied by the area of elevation bin centred at the contour (using the hypsometry of the larger glacier) to calculate volume change in a particular zone. Volume changes in each elevation bin were then summed to obtain dV. This method seemed to underestimate surface elevation change, as contour lines tend to be more curved (i.e., have more vertices) in the lateral zones of glaciers, where elevation change is smaller than in the glacier centre (Nuth et al. 2007).

As a measure of average surface elevation change, d h ¯ , between 1960 and 1990, Method A is used due to the superior results.

Geometry changes for the 1990–2009 epoch were estimated by subtracting the co-registered and aligned 2009 DEM from the 1990 DEM. For DEMs with 20 m grids, the subtraction is straightforward. However, for BLB and FDB a 2009 DEM with finer resolution (8 m) was available, so before subtraction the 20 m 1990 DEM is first resampled to 8 m resolution, using bilinear interpolation.

From the resulting rasters, total volume change of each glacier, dV, was first computed. Glacier-wide surface lowering, d h ¯ , was inferred by dividing dV by the average area of a glacier in the studied period. Even though the dynamics are considered as being of low importance in elevation changes, interpretation of the results needs to be undertaken with caution. However, when dh information is integrated over the whole glacier into d h ¯ , the dynamic component of surface motion in x and y directions cancels, if steady state is assumed (Paterson 1994).

To evaluate errors for Method A, it was used to calculate d h ¯ in the second studied period (1990–2009), for which DEMs are available. The results of the tests were in good agreement with those obtained by DEM differencing. Residuals ranged from −0.62 to 0.55 m. If weighted by area of the studied glaciers, the mean difference in d h ¯ between Method A and DEM differencing was −0.24 m, with a standard deviation of 0.42 m. In terms of percentage, Method A gave on average more negative values by 2%, with differences ranging from −5% to 4% for individual glaciers, in relation to DEM-obtained results. Therefore, it is assumed that overall performance of Method A is similar to DEM differencing.

To estimate overall errors, elevation differences from non-glacier covered terrain were investigated using 18 000 (1960–1990) and 35 000 (1990–2009) test-points. Because this work focuses on very small glaciers, where relatively steep sections (potentially causing large elevation errors) comprise significant portions of their area, it was important to pay particular attention to slope distribution. In this error assessment, it is assumed that elevation errors may be described by slope and surface type. Hence, conventional elevation error histograms (Fig. 2a) only partially represent uncertainties of the available data sets because they are inferred from non-glacier covered mountainous terrain. Ice and snow surfaces are characterized by greater elevation errors than rock terrain due to low radiometric contrast on aerial images used for map/DEM construction. In this study, glacier ice with good contrast has error characteristics similar to that of non-glacier covered terrain, while for snow-covered areas the error is assumed to be doubled (Nuth et al. 2010).

Test-points were grouped in 1° slope bins, for which root mean square error was calculated, representing mean elevation error of an individual point at a given slope (E _p ). E _p increases with slope and is well described by a third order polynomial function (Fig. 2b). Slope histograms for the glaciers were extracted, and the polynomial functions were used to convert these to E _p histograms. Then, for each glacier and epoch, a weighted mean of E _p was calculated to obtain error values representative for the whole population of points lying within each glacier's boundaries (E _g ). The total error of the elevation change measurement of each glacier ( E d h ¯ ) was calculated by dividing E _g by the square root of the sample size (N): 1 E d h ¯ = E g N Elevation errors are spatially autocorrelated, so it was not appropriate to use the number of sample points; rather, 1000 m as a correlation scale was used (e.g., Nuth et al. 2007; James et al. 2012), and N becomes the glacier area in km². Total elevation change errors computed with this method ranged from ±3 to ±15 m for individual glaciers.

Elevation changes in the first analysed period, 1960–1990, are calculated using different procedures and compared in Table 2. The final maps of elevation changes from both periods are shown in Fig. 3, and the results are quantified in Table 3, together with volume changes and error estimates. Mean annual change in 50 m elevation bins for both study intervals is compared in Fig. 4.

The results show a typical pattern for glaciers with a negative mass balance: high thinning at fronts with limited thickening (1960–1990) or no thickening at all, even in the highest parts (1990–2009). The greatest thinning was on VMB, with −0.66 and −0.96 m a⁻¹, respectively, for the first and second periods, for a total of −39.3 m elevation change over the entire glacier for 1960–2009. NNB thinned in total by 35.3 m with high rates of −0.68 and −0.81 m a⁻¹. Other glaciers experienced less thinning, ranging from −0.32 to −0.53 m a⁻¹ for 1960–1990 and from −0.61 to −0.82 m a⁻¹ for 1990–2009. For all glaciers, the latter period was characterized by significantly higher thinning rates compared to the first period. The area-weighted mean of d h / d t ¯ in 1960–1990 was −0.49 m a⁻¹, with some thickening above 575 m, while after 1990 it increased by 54% to −0.78 m a⁻¹ and no elevation increase occurred even at higher elevations.

In terms of volume, all glaciers have experienced significant losses (Table 3). VMB volume change (dV) since 1960 was the most prominent, with −134 · 10⁶ m³, and a ca. 20% acceleration in volume loss rate (dV/dt) between the two epochs. Total volume loss for BLB was 133 · 10⁶ m³, with a similar acceleration in dV/dt (26%). SVB lost 84 · 10⁶ m³ of ice in total, experiencing the greatest acceleration in volume loss rate (80%) in the second period. Only two of the smaller glaciers, NNB and ELB, show a deceleration in the volume loss rate, decreasing by +5% and +2%, respectively, most likely due to decrease of glacier area exposed for melt. The total dV for all the study glaciers is −517 · 10⁶ m³ since 1960, with a 52% increase in dV/dt after 1990. This represents a roughly 39% volume loss, based on 1960 volume estimates by Hagen et al. (1993).

For the largest study glaciers, the increase in thinning rates after 1990 shows a clear trend with altitude (Fig. 5). In the lower and middle elevation bands, dh/dt was relatively stable, with thinning increased by 0.2–0.4 m a⁻¹. At higher elevations, the thinning rates increased the most, by up to 1.1 m a⁻¹. For the three smallest glaciers, no such trend is visible.

Earlier studies of other small Svalbard glaciers concluded that climate shift is the main driver of elevation changes (e.g., Kohler et al. 2007; Barrand et al. 2010; James et al. 2012). From the beginning of the study interval, in 1960, summer temperature at Svalbard airport (50 km from the study area) has been rising by 0.02°C a⁻¹; since 1990, the trend is even more pronounced, at 0.07°C a⁻¹ (James et al. 2012). There is no significant trend in snow precipitation and accumulation, although there may be a slight recent decrease (Pinglot et al. 1999; Pohjola et al. 2002; Førland & Hanssen-Bauer 2003), so glaciers receive an increasing energy available for melt which is not balanced by higher accumulation input. Therefore, a strong link between climate and negative volume changes of the glaciers studied in this paper is anticipated.

The observed trends of general glacier wastage are similar to those found in previous work in other Svalbard regions (e.g., Dowdeswell et al. 1997; Hagen, Kohler et al. 2003; Hagen, Melvold et al. 2003; Nuth et al. 2007; Rachlewicz et al. 2007; Sobota 2007Sobota 2011; Zagórski et al. 2008; Moholdt et al. 2010; Nuth et al. 2010). Unfortunately, the number of papers which could be used for a direct interregional comparison is limited, mainly due to differences in temporal coverage of earlier works and methods involved. The published long-term (1960s − 2000s) elevation changes of individual Svalbard ice masses from different regions were −0.55 m a⁻¹ and −0.61 m a⁻¹ (Kääb 2008), −0.41 and −0.58 m a⁻¹ (Barrand et al. 2010), −0.25, −0.32, −0.44, −0.64 and −0.87 m a⁻¹ (James et al. 2012, excluding glaciers with partial data coverage). In this paper, the average surface lowering in central Dickson Land for the period 1960–2009 was −0.61 m a⁻¹ (Table 3). While it is true that the mass loss of glaciers in this quasi-continental study area is faster than on the larger, more coastally influenced Svalbard ice masses (−0.36 m a⁻¹; Nuth et al. 2010), there is no evidence of more negative than on small valley glaciers elsewhere in the archipelago. It may be explained by the higher elevation of Dickson Land ice masses, which partly offsets high summer temperature and limited precipitation.

An important finding of this study is the apparent acceleration of thinning after 1990. In the first studied period, the mean zero-change line was situated at 575 m. In the more recent epoch, the average dh/dt did not approach zero at all, implying that the geodetic equilibrium has shifted by about 150 m up-glacier. In effect, overall glacier d h / d t ¯ in the period 1990–2009 was on average 58% more negative than between 1960 and 1990. The thinning increase is particularly marked in the upper zones of the study glaciers, in agreement with Kohler et al. (2007) and James et al. (2012). Fig. 5 shows that the most negative changes occur above 500–600 m a.s.l., corresponding to the approximate position of zero-change line between 1960 and 1990. While the traditionally defined equilibrium line altitude may be situated below this elevation (if ice submergence is not balanced by mass gain), it is assumed here that above the geodetic equilibrium elevation the only source of surface elevation increase would be net accumulation. Therefore, the observed post-1990 acceleration of surface lowering in areas that previously have been experiencing net mass surplus may indirectly support the hypothesis that this process is caused by albedo change due to faster snow removal under the warming climate (Kohler et al. 2007; James et al. 2012).

The highest zones of glaciers used to be covered with snow or firn throughout the summer season, but warming after 1990 has most likely been progressively removing snow to uncover ice. When ice is exposed, its lower albedo increases energy absorption and therefore melting. Albedo in the lower ablation zone remains at the same levels, so thinning rates there are more stable. This mechanism is in agreement with direct observations by the author from the last decade when little or no snow remains during summers, so that firn and ice melt, even in the highest areas of glaciers. The role of strong dh/dt acceleration cannot be neglected in future projections of Svalbard contribution to SLR, but further detailed observations are needed, as there are other possible mechanisms. As outlined in James et al. (2012), those may include a shift in atmospheric stratification, albedo change due to debris incorporation from the rock walls, increased ablation due to lack of internal accumulation, or change in energy balance conditions of glacial valleys due to ice area decay.

This study revealed that neighbouring glaciers of similar size and type may show very different dh/dt curves, illustrating the problem of sample size and quality when assessing regional mass balance trends. Important factors modulating general mass loss are aspect and area-altitude distribution. Between 1960 and 2009, the greatest changes occurred either on glaciers receiving the highest amounts of solar radiation or on glaciers with more area at lower elevations. Therefore, VMB, NNB, BLB and AMB all had very high overall mass loss, while SVB, FDB and ELB have been thinning at slower rates. This underscores the importance of careful site selection for region-wide mass balance assessments. Ideally one should measure elevation change on all glaciers in a particular area to obtain a convincing regional mass balance.

This study focuses on the little investigated central Svalbard glaciers of Dickson Land. Between 1960 and 2009, the overall d h ¯ of seven glaciers analysed here ranged from −20.9 to −39.3 m, equivalent to a roughly 39% volume loss. In the first period, 1960–1990, the study glaciers have thinned at fronts with some limited thickening in their upper parts, with d h / d t ¯ =−0.49 m a⁻¹. In response to rising summer temperatures, the zero-change elevation increased by about 150 m, and the thinning rates in all elevation bands increased by 58% for d h / d t ¯ = − 0.78 m a⁻¹. Between neighbouring and somewhat similar glaciers, significant differences in mass change are observed, underscoring the importance of choosing a representative sample when investigating mass balance of whole regions, which should particularly consider aspect and hypsometry. Glaciers in the study area have been losing mass at a similar rate as other small, land-terminating glaciers in different Svalbard regions, even though the local climate is more continental than in other areas. This can be explained by the relatively high elevation of local ice-masses, which partly offsets the unfavourable climatic conditions. The results also show a strong increase in thinning rates in the upper elevations of glaciers, as has been found for other areas of Svalbard, and which has been attributed to albedo decreasing in the accumulation zones.