Geodetic mass balance record with rigorous uncertainty estimates deduced from aerial photographs and lidar data – Case study from Drangajökull ice cap, NW Iceland

Abstract

In this paper we describe how recent high-resolution digital elevation models (DEMs) can be used to extract glacier surface DEMs from old aerial photographs and to evaluate the uncertainty of the mass balance record derived from the DEMs. We present a case study for Drangajökull ice cap, NW Iceland. This ice cap covered an area of 144 km${}^{\mathrm{2}}$ when it was surveyed with airborne lidar in 2011. Aerial photographs spanning all or most of the ice cap are available from survey flights in 1946, 1960, 1975, 1985, 1994 and 2005. All ground control points used to constrain the orientation of the aerial photographs were obtained from the high-resolution lidar DEM. The lidar DEM was also used to estimate errors of the extracted photogrammetric DEMs in ice- and snow-free areas, at nunataks and outside the glacier margin. The derived errors of each DEM were used to constrain a spherical semivariogram model, which along with the derived errors in ice- and snow-free areas were used as inputs into 1000 sequential Gaussian simulations (SGSims). The simulations were used to estimate the possible bias in the entire glaciated part of the DEM and the 95 % confidence level of this bias. This results in bias correction varying in magnitude between 0.03 m (in 1975) and 1.66 m (in 1946) and uncertainty values between $\pm$0.21 m (in 2005) and $\pm$1.58 m (in 1946). Error estimation methods based on more simple proxies would typically yield 2–4 times larger error estimates. The aerial photographs used were acquired between late June and early October. An additional seasonal bias correction was therefore estimated using a degree-day model to obtain the volume change between the start of 2 glaciological years (1 October). This correction was largest for the 1960 DEM, corresponding to an average elevation change of $-$3.5 m or approx. three-quarters of the volume change between the 1960 and the 1975 DEMs. The total uncertainty of the derived mass balance record is dominated by uncertainty in the volume changes caused by uncertainties of the SGSim bias correction, the seasonal bias correction and the interpolation of glacier surface where data are lacking. The record shows a glacier-wide mass balance rate of$\dot{B}=-0.26\pm 0.04$ m w.e. a${}^{-\mathrm{1}}$ for the entire study period (1946–2011). We observe significant decadal variability including periods of mass gain, peaking in 1985–1994 with $\dot{B}=0.27\pm 0.11$ m w.e. a${}^{-\mathrm{1}}$. There is a striking difference when $\dot{B}$ is calculated separately for the western and eastern halves of Drangajökull, with a reduction of eastern part on average $\sim$ 3 times faster than the western part. Our study emphasizes the need for applying rigorous geostatistical methods for obtaining uncertainty estimates of geodetic mass balance, the importance of seasonal corrections of DEMs from glaciers with high mass turnover and the risk of extrapolating mass balance record from one glacier to another even over short distances.

Introduction

Mountain glaciers and ice caps accounted for more than half of the land ice runoff contribution to global mean sea-level rise during the 20th century (Vaughan et al., 2013). Understanding how these glaciers respond to a changing climate is essential to close the budget of the sea-level rise over the last decades and project the sea-level rise in the near future. In recent years an increased part of our knowledge about how these glaciers are changing has been based on remote sensing. The majority of these studies describe current or recent glacier changes in different parts of the globe by applying geodetic methods (e.g., Gardelle et al., 2012; Berthier et al., 2010). Others have presented results on the geodetic mass balance extending further back (e.g., Fischer et al., 2015; Nuth et al., 2007; Soruco et al., 2009); these studies are particularly important since they indicate how the glaciers responded to 20th century climate variability. Such observations can be used to constrain or correct glacier mass balance models that are used to estimate how the glaciers will respond to future climate changes (e.g., Clarke et al., 2015).

Studies on long-term geodetic mass balance are generally based on digitized contour maps, with some exceptions where mass balance records have been derived from digital elevation models (DEMs) extracted from old archives of aerial photographs by applying digital photogrammetry (e.g., James et al., 2006, 2012). The applicability of geodetic mass balance records as a key to predicting future glacier changes depends on the accuracy of such records and their resolution. To maximize both the accuracy and the resolution we should rather focus, if possible, on archives of aerial photographs, because

• i.

  these archives often span more epochs than the published topographic maps;

• ii.

  with new and rapidly improving tools in digital photogrammetry the potential to produce more accurate and detailed DEMs than those deduced by interpolating elevation contours from old maps has increased significantly;

• iii.

  the availability of high-resolution DEMs has opened a new source of ground control points (GCPs) for constraining the orientation of photogrammetric DEMs (James et al., 2006; Barrand et al., 2009). Like (ii), this will lead to more accurate DEMs from aerial photograph archives in future studies. New spaceborne sensors such as Worldview and Pléiades may allow such studies in remote areas without conducting expensive field campaigns to survey GCPs (Papasodoro et al., 2015).

In order to maximize the value of geodetic mass balance records, realistic uncertainty assessments are required. If the uncertainty is overestimated, the value of the information that we can extract from the geodetic data will be diminished, and the results will be neglected by the scientific community or not even be published. If, however, the uncertainty is underestimated, geodetic mass balance records with significant errors will be interpreted as reliable observations. When extracting volume change from two different DEMs, a common approach is to use the standard deviation of the DEM difference in the unglaciated part of the DEMs as a proxy for the uncertainty of the average elevation change (e.g., Cox and March, 2004). This method corresponds to an extreme case, assuming that the errors of the surface elevation change are totally correlated between all grid cells within the glacier. The opposed extreme case assuming that the errors of surface elevation change are totally uncorrelated between all grid cells has also been applied in the literature (e.g., Thibert et al., 2008). This approach results in an estimated uncertainty reduced by a factor $\sqrt{n}$ compared to the totally correlated uncertainty where n is the number cells for which the difference is calculated. The third alternative, where the spatial dependence of the DEM errors is estimated and inherent in the uncertainty estimate, was described by Rolstad et al. (2009). This method results in uncertainty somewhere between the two extremes and has been adopted in several studies (e.g., Trüssel et al., 2013; Zemp et al., 2013; Fischer et al., 2015). This method includes some simplifications, which so far have not been validated with other geostatistical methods.

Location of study area. Blue lines in (a) are the outline of the larger glaciated areas in Iceland. The triangles in (b) indicates the locations of the meteorological stations at Æðey and Gjögur. Image (c) shows a lidar DEM of Drangajökull (glacier margin shown with blue line) and vicinity obtained in 2011 (Jóhannesson et al., 2013) represented as shaded relief image and contour map (100 m contour interval). The names and locations of the three main outlet glaciers are shown.

[Figure omitted. See PDF]

Here, we present a case study of Drangajökull ice cap in NW Iceland (Fig. 1) based on seven sets of aerial photographs in 1946–2005 and a lidar DEM obtained from an airplane in 2011 (Jóhannesson et al., 2013). The glacier covered an area of 144 km${}^{\mathrm{2}}$ in 2011 and is the fifth largest glacier in Iceland. This study describes an alternative method to estimate uncertainties of the average elevation change derived by differencing DEMs, applying geostatistical methods. The approach, which uses the DEM difference from ice- and snow-free areas as input, allows for a simultaneous estimate of a bias correction for the glaciated part of the DEMs. Both the estimated uncertainty and the bias correction are compared with results from conventional methods. We also interpolate volume changes in areas where data are lacking and inspect how much of the derived volume change may be caused by seasonal variation. The study results in a seasonally corrected mass balance record of Drangajökull ice cap with estimates of possible errors contributing to the record as well as the derived net uncertainty.

Dates, main parameters and notes describing the data sets used in the study.

Date	No. images	Average GSD${}^{\ast }$ (m)	Notes
12 Oct 1946	15	0.94	Not covering the southernmost part of Drangajökull. Over-exposed areas.
Summer 1960	40	0.42	Divided in three flights: 14 Jun 1960, 8 Jul 1960 and 12 Jul 1960
5 Sep 1975	18	0.77	Not covering the lower part of Leirufjarðarjökull
27 Jul 1985	32	0.70	Lowermost part of Reykjarfjarðarjökull obscured by clouds
4 Aug 1986	5	0.70	Used for filling the gaps of 1985 on Reykjarfjarðarjökull
29 Aug 1994	21	0.53	The southern part of Drangajökull is not covered
27 Jul 2005	57	0.53	Complete coverage
20 Jul 2011	–	–	Complete coverage (lidar)

${}^{\ast }$ GSD: ground sampling distance.

Data and methods

In this study, seven sets of aerial photographs covering Drangajökull ice cap in 1946, 1960, 1975, 1985, 1986 and 1994 from the archives of the National Land Survey of Iceland and in 2005 from Loftmyndir ehf were used. Negative films were scanned with a photogrammetric scanner in a resolution of 15 and 20 $\mathrm{\mu }$m. The aerial photographs have an average scale between $\sim$ $\mathrm{1}:\mathrm{30\hspace{0.17em}000}$ and $\sim$ $\mathrm{1}:\mathrm{40\hspace{0.17em}000}$, which result in a ground sampling distance (GSD) of $\sim$ 0.4 to $\sim$ 1 m. Complete camera calibration information is available for the surveys of 1975, 1985, 1986, 1994 and 2005, but calibration information is lacking for the oldest flights (1946 and 1960). Only the focal length is available for the photographs of 1946, and focal length and radial distortion are available for the photographs of 1960. Table 1 summarizes the main characteristics of each series.

During the International Polar Year 2007–09, a major effort was initiated to produce accurate DEMs of all the major Icelandic glaciers and ice caps (Jóhannesson et al., 2013). In July 2011 Drangajökull ice cap was surveyed with airborne lidar model Optech ALTM 3100. The lidar DEM covers the entire ice cap as well as the close vicinity of the glacier, which provides a useful reference to constrain and validate the other DEMs produced in this study. Specifications of the survey are described in Jóhannesson et al., 2013. The average density of the point cloud measured with the lidar corresponded to 0.33 hits m${}^{-\mathrm{2}}$. The high density facilitates a well-constrained bilinear interpolation of the point cloud into a grid with 2 m $\times$ 2 m cell size. Cells where the distance to nearest lidar hit exceeds 4 m were masked out. A comparison of differential GPS profiles and a 5 m $\times$ 5 m grid derived from identical lidar survey in the Snæfelljökull ice cap in western Iceland indicated vertical accuracy well within 0.5 m (Jóhannesson et al., 2011).

Creation of DEMs and orthorectified photographs

The DEMs were created from the aerial photographs using the software bundle IMAGINE Photogrammetry (^©Intergraph). The photogrammetric processing is carried out in four steps: orientation of the images, automatic stereo matching, manual editing of the DEMs and orthorectification of aerial photographs.

The coverage of aerial photographs (black line boxes) at different epochs with the lidar DEM as background. The GCPs used for orientation of each series of aerial photographs are marked with triangles.

[Figure omitted. See PDF]

Each series of aerial photographs was oriented individually by means of a rigorous bundle adjustment (Wolf and Dewitt, 2010). The glacier is covered by a single series of images for all years except in 1960 when the glacier was covered by three tiles, one per date (Table 1). Tie points were automatically measured in the images and semiautomatically revised, ensuring a good connection between all the adjacent photographs and between strips. The exterior orientation was constrained by using series of GCPs extracted from the lidar DEM (2 m $\times$ 2 m cell size), applying a similar approach to the one carried out by James et al. (2006). The lidar DEM was viewed as a hillshade with approximately the same sun position as during the acquisition of the photographs. This allowed recognition and extraction of GCPs from stable features such as boulders and sharp edges in the ice-free areas. To ensure stability in the orientation, a fairly regular distribution of GCPs over the photographed area as well as over the elevation span of the terrain is required (Kraus, 2007; Nuth and Kääb, 2011). Artificial dip or rise in the DEM due to insufficient coverage of GCPs would skew the geodetic mass balance record and make its uncertainty estimate as explained below less valid (further explained in Sect. 2.2). In our case the nunataks of Drangajökull ice cap allow fairly even spatial and vertical distribution of GCPs for all epochs (GCP locations shown for each DEM in Fig. 2). The photogrammetric orientations performed in this study never span more than two photographs without having constraints from a GCP. This is considered to be sufficient coverage of GCPs for a reliable orientation (Kraus, 2007). The assigned uncertainty of the GCPs used in the DEM processing was 2 m standard deviation in $XY$ and 0.5 m for $Z$ corresponding, respectively, to the lidar DEM resolution and expected vertical accuracy. A significant part the large-scale errors in the derived photogrammetric DEM may be related to errors in the GCPs 3-D locations, particularly for DEMs from 1975 and later when other data constraining the geometric model are relatively accurate. For the 1946 and 1960 DEMs a lack of camera calibration information is likely to be a more important source of errors.

The orientation of the 1960 images was carried out using the focal length and lens distortion information obtained from the calibration report of the DMA Cameras (Spriggs, 1966). The 1946 images included information of the focal length written at the margin of the first image of each strip. Both cases needed auxiliary pre-calibration, and therefore pseudo-fiducial marks were created to allow the location of a pseudo-principal point (see, e.g., Kunz et al., 2012, for details). The orientation of both sets of images included additional parameters in the bundle adjustment for refinement of the camera geometry. Bauer's model (Bauer and Müller, 1972) was used for the images of 1946 and Jacobsen's model (Jacobsen, 1980) was used for the images of 1960.

Once oriented, we produced the elevation point clouds from stereo-matching of the images. The routine eATE (enhanced Automatic Terrain Extraction) of the software allows for a pixel-wise evaluation in the matching process, thus obtaining a high density of points. The low contrast in firn- and snow-covered areas caused failures in the matching process. The point clouds for low-contrast areas were therefore created from reduced resolution of the stereo images and a larger window size and lower correlation coefficient of the stereo matching. This resulted in an improved coverage of points automatically measured in the snow-covered areas. A first edition of the point clouds was carried out with the software CloudCompare (GPL Software); automatic outlier removal was performed using the routine “Statistical Outlier Removal” (Rusu and Cousins, 2011). The dense point clouds were then subsampled in regular density of points corresponding to $\sim$ 10 m $\times$ 10 m spacing for all epochs except 1946 and 1985, to which density equivalent to $\sim$ 20 m $\times$ 20 m spacing was applied. This was done to reduce the size of the point clouds and remove double points that could introduce noise when interpolating the point clouds as a grid with fixed cell size (Sect. 2.2 and 2.3). The lower subsampled point density was due to large GSD in the case of 1946 and a high level of noise in the images in 1985, resulting in large amount of outliers. Finally, a thorough manual revision of the results in stereoscopic vision was carried out, editing the DEMs in the glacier areas where the automatic matching failed and surface details were still perceptible.

To delineate the glacier margin and mask out snow-covered areas (Sect. 2.2 and 2.4) orthorectified photographs were required. The orthorectification was carried out using preliminary DEMs linearly interpolated from the point clouds as grids with 10 m $\times$ 10 m (DEMs of 1960, 1975, 1994 and 2005) and 20 m $\times$ 20 m cell size (DEMs of 1946 and 1985). The series of 1975 included two strips spanning the glacier (Fig. 2) without covering it completely in stereo. These images were orthorectified using the lidar DEM, revealing the location of the glacier margin at its intersection with the bare ground, free of ice and snow in both 1975 and 2011 (resulting in insignificant elevation change at the 1975 margin location in this area). The orthorectification of all the series of photographs was performed in resolution corresponding to a 2 m $\times$ 2 m pixel size.

DEM error assessment and bias correction

We use the high-resolution lidar DEM obtained in 2011 to assess the quality of the photogrammetric DEMs. The photogrammetric DEMs are expected to be of significantly worse quality in terms of accuracy than the lidar data and we therefore assume for simplicity that statistical parameters derived from the difference between the photogrammetric DEM and the lidar DEM (in areas assumed stable) describe errors in the photographic DEM. This is likely to produce a minor underestimate of the actual quality of the photographic DEMs. As described below, all photogrammetric DEMs were bias corrected relative to the lidar DEM. A possible bias in the absolute location of the lidar DEM does not affect our result since the bias is canceled out when calculating the difference between the DEMs.

The first step in estimating the quality of a DEM derived from the aerial photographs was calculating the difference between the photogrammetrically derived point clouds (Fig. S1 in Supplement) and the lidar DEM with 2 m $\times$ 2 m cell size. This was calculated using the residual operation in Surfer 12 (^©Golden Software, Inc.). From this a digital model of the difference between the DEMs was linearly interpolated for a grid with 20 m $\times$ 20 m cell size. All cells with snow or glacier cover at either or both dates (photograph and lidar acquisitions) were masked out as well as cells where distance to the next element of the point cloud exceeds 40 m. The glacier outlines were delineated manually (see Sect. 2.4) and the snow-covered areas were derived with semiautomatic classification of the orthorectified aerial photographs and the intensity images derived from the lidar scanning. The mean and the standard deviation ($\mathit{\sigma })$ of the derived difference (photogrammetric DEM$-$lidar DEM) for the remaining data after snow and glacier masking is shown in Table 2.

Extraction of geodetic mass balance requires co-registered DEMs prior to calculation of glacier volume changes. This usually includes estimates of relative vertical and horizontal shift between the DEMs, using areas where the elevation change is expected to be insignificant (Kääb, 2005; Nuth and Kääb, 2011; Guðmundsson et al., 2011). In this study the GCPs used during the orientation of the photographs were extracted from the lidar DEM in maximum resolution (2 m $\times$ 2 m cell size). We were able to extract several GCPs at nunataks near the glacier center. The distribution of GCPs is therefore fairly regular over the survey area in all cases both spatially (Fig. 2) and with elevation. The orientation of aerial photographs resulted in horizontal RMSE of the GCPs < 3 m in all cases and typically 1–2 m (Table 2). These values are obtained from least square adjustment resulting in residual mean equal to 0. The horizontal shift relative to the lidar DEM is likely to exceed the derived horizontal RMSE locally for a given photogrammetric DEM. It is, however, unlikely that the average horizontal shift relative to the lidar DEM exceeds the derived horizontal RMSE of the GCPs. We therefore concluded that horizontal shift corrections are not required for the photographic DEMs.

The series of DEMs of Drangajökull ice cap created from the aerial photographs. The shaded relief images and contour maps indicate the glaciated part of each DEM. The elevation difference off ice (after masking out outliers and areas with slope > 20${}^{\circ }$) are shown as color images. The color scale is extended for the DEM in 1946 and reduced for the 1994 and 2005 DEMs. A vertical histogram next to the scale bar shows the error distribution.

[Figure omitted. See PDF]

The elevation difference between DEMs covering stable areas is commonly used to estimate zero order (bias correction, see, e.g., Nuth and Kääb, 2011; Guðmundsson et al., 2011) or higher order correction (e.g., Rolstad, 2009; Nuth and Kääb, 2011) to compensate for slowly varying errors in DEM difference over glaciated areas. The result from such approach is, however, sensitive to the area chosen as the reference area. One can choose to use the entire area covered by both DEMs outside the glacier or an area limited by a certain distance from the glacier. In this study we apply geostatistical methods for deriving bias correction of each photogrammetric DEM within the glacier and an estimate of the uncertainty in the derived bias correction. These calculations consisted of the following five main steps.

• Preparation of DEM error input data (derived from the comparison with the lidar), explained below. Resulting error data from ice- and snow-free areas are shown in Fig. 3.

• Transformation of the derived DEM errors into a new variable with the nscore function (Deutsch and Journel, 1998) in WinGSlib V.1.5.8 (^©Statios LLC). The histogram of the new variable fits a normal distribution, with zero mean and $\mathit{\sigma }=\mathrm{1}$. This step is a recommended preparation of a data set for valid sequential Gaussian simulation (SGSim) carried out in step 5 particularly if the histogram of the DEM error does not closely resemble a normal distribution.

• Calculation of semivariogram for the nscored input data, in which the semivariogram describes the variance, $\mathit{\gamma }$, of a given coordinate-based variable as a function of distance, $d$, between sampled locations.

• Calculation of a spherical semivariogram model, fitting the derived semivariogram.

• Use of the derived spherical model and the nscored data that constrain the semivariogram to run 1000 SGSims of the nscored errors in the glaciated areas using the sgsim function (Deutsch and Journel, 1998) in WinGSlib. The sgsim function includes reversed transformation from the nscored variable to the derived DEM error. SGSims are commonly applied in errors assessments of geostatistical studies (e.g., Lee et al., 2007; Cardellini et al., 2003). The results from the sgsim runs were used to estimate both the most likely bias of each photogrammetric DEM within the glacier and 95 % confidence level of this bias, as explained further below.

The approach adopted here requires that the statistics of the DEM errors outside the glacier be descriptive for the errors in the photogrammetric DEM within the glacier margin. This should be kept in mind both during the photogrammetric processing and in the preparation of input data (step 1) used in the geostatistical calculation. The photogrammetric processing requires fairly even spatial distribution of GCPs; otherwise artificial dip or rise in the photogrammetric DEM are likely to be produced in areas far from a GCP (Kraus, 2007). Such errors would not be represented in a semivariogram based on DEM error in areas where distribution of GCPs is adequate.

The horizontal RMSE of the GCPs (no. of GCPs within brackets), glacier coverage and error assessment of the photogrammetric DEMs, using four different approaches: (i) direct comparisons of ice-free areas (mean and standard deviation ($\mathit{\sigma }$)); (ii) comparisons in ice-free areas after masking out outliers and areas with slope > 20${}^{\circ }$ (see Sect. 2.2); (iii) SGSim, in which $z$_bias corresponds to the mean elevation bias from 1000 simulation, $z$_bias${}_{\mathrm{u}}$ and $z$_bias${}_{\mathrm{l}}$ to the upper and lower 95 % confidence level and $\mathrm{\Delta }$$z$_bias $=$ ($z$_bias${}_{\mathrm{u}}-z$_bias${}_{\mathrm{l}})/\mathrm{2}$; (iv) method described by Rolstad et al. (2009). To derive uncertainties with 95 % confidence level we assume normal probability function and therefore $\mathrm{\Delta }z$_bias${}_{\text{Rols}}=1.96\times {\mathit{\sigma }}_{z\text{\_bias\_Rols}}$.

Year	RMSE XY	Glacier cover-	Error mean	$\mathit{\sigma }$ ice-	Error mean ice-	$\mathit{\sigma }$ ice-free	$z$_bias	$z$_bias${}_{\mathrm{l}}$	$z$_bias${}_{\mathrm{u}}$	$\mathrm{\Delta }z$_bias	$\mathrm{\Delta }z$_bias_Rols.
	GCPs (m)	age (%)	ice-free (m)	free (m)	free masked (m)	masked (m)	(m)	(m)	(m)	(m)	(m)
1946	2.99 [43]	75.3	$-$0.95	5.09	$-$0.86	4.80	1.66	0.12	3.27	1.58	3.41
1960 W	2.87 [25]	31.0	0.37	2.23	0.49	1.84	0.48	$-$0.34	1.34	0.84	1.05
1960 C	2.54 [31]	30.5	$-$0.31	2.08	$-$0.26	1.52	0.34	$-$0.29	1.02	0.66	1.04
1960 E	2.21 [47]	35.6	0.03	2.26	0.09	1.51	0.20	$-$0.45	0.93	0.69	0.96
1975	1.22 [44]	96.5	0.48	2.05	0.39	1.52	0.03	$-$0.47	0.48	0.48	0.62
1985	1.37 [33]	87.2	$-$0.67	1.97	$-$0.60	1.15	$-$0.48	$-$0.80	$-$0.17	0.32	0.47
1994	0.84 [40]	66.3	$-$0.09	1.04	$-$0.09	0.80	0.22	$-$0.03	0.47	0.25	0.72
2005	1.14 [55]	100.0	$-$0.24	1.30	$-$0.26	0.87	0.22	0.01	0.42	0.21	0.78

The low contrast of snow-covered glacier surface may also result in a difference in error statistics between the glacier and the ice- and snow-free areas (Rolstad et al., 2009). The low contrast should mostly produce high-frequency errors, whereas low-frequency errors are mostly caused by an inaccurate orientation. The eATE configuration used resulted in fewer but better matching points in the low-contrast areas (Sect. 2.1) and the thorough manual 3-D revision likely removes most of the high-frequency noise in the resulting DEM. A semivariogram of the difference between the point cloud in 1946 at low-contrast glacier areas and the lidar DEM (blue crosses in Fig. 4c) reveals the variance with distance for the elevation error plus the elevation changes in 1946 to 2011. The variance of elevation changes over a short distance should be small for smooth glacier surface. At short distances the semivariogram should therefore mainly represent the DEM errors. For $d$ < 200 m the low-contrast areas show variance at similar level as for the DEM error data outside the glacier (Fig. 4c), indicating similar level of high-frequency error for the two area types. This supports that the errors in low-contrast areas are unlikely to significantly skew our geostatistical analyses.

A difference in terrain slope between areas can produce a significant difference in the calculated semivariogram (Rolstad et al., 2009). Local horizontal shift between DEMs can produce significant artificial elevation difference in steep areas. The average slope on the glacier in 2011 was 6.2${}^{\circ }$ whereas the unglaciated area in the 2011 lidar DEM had an average slope of 9.8${}^{\circ }$. The preparation of our data (step 1) therefore includes exclusion of all data where slope exceeds 20${}^{\circ }$; unglaciated area in the 2011 lidar DEM below this slope limit has an average slope of 7.2${}^{\circ }$.

The glaciated parts of the photogrammetric DEMs were all manually revised using 3-D vision, securing removal of significant outliers within the glacier. A thorough revision was not carried out for the unglaciated areas. Instead we apply automatic removal of outliers. This was carried out by calculating the standard deviation of the DEM error (photogrammetric DEM$-$lidar DEM), ${\mathit{\sigma }}_{\mathit{\epsilon }h}$ (after masking out snow-covered, glacier-covered and steep areas) and filtering the DEM error with a 500 m $\times$ 500 m median filter. Values in the photogrammetric DEM where the difference between the unfiltered and the median filtered error exceeded ${\mathit{\sigma }}_{\mathit{\epsilon }h}$ were then masked out. The mean DEM error and ${\mathit{\sigma }}_{\mathit{\epsilon }h}$ after the slope and outlier masking is shown in Table 2.

The semivariograms obtained with (step 3) and without the nscore transformation of the 1946 DEM error in ice- and snow-free areas are shown in Fig. 4a–b. The spherical semivariogram model calculated in step 4 is given as a function of $d$ (distance between sampled locations): $$\begin{array}{cccc}\mathit{\gamma }\left(d\right)& =\mathrm{0}& & d=\mathrm{0}\\ & =c_{\mathrm{0}}+c_{\mathrm{1}}\left[{\displaystyle \frac{\mathrm{3}d}{\mathrm{2}r}},-,{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},{\left({\displaystyle \frac{d}{r}}\right)}^{\mathrm{3}}\right]& & \mathrm{0}<d\le r\\ & & =c& & d>r.\end{array}$$ $c=c_{\mathrm{0}}+c_{\mathrm{1}}$ and $\mathit{\gamma }\left(\mathrm{0}\right)$ describe the correlation of a point with itself. The main parameters in the model, nugget ($c_{\mathrm{0}})$, range ($r)$ and sill ($c)$ are shown in the Fig. 4b. We expect $c$ to equal approximately the global variance of the data set; hence $c=\sim \mathrm{1}$ for the nscored data. The shape of the semivariograms that we obtain (Fig. 4a–b and Fig. S2) indicates a reasonable fitting with a single spherical model, unlike in the study by Rolstad et al. (2009) where two spherical models describing the variance at different ranges of distances were required.

The size of the DEM error grid in full resolution (20 m $\times$ 20 m cell size) was too large for the sgsim function to operate (step 5). The data size was reduced by picking out every fifth column and line in the DEM error grid. In areas where data were sparse, at nunataks and where few data points remained due to the snow mask near the glacier margin, the 20 m $\times$ 20 m data were used. Tests with smaller study areas indicated that this reduction of the input data only has minor effects on the results derived from the simulation.

The semivariograms of the 1946 DEM error before (a) and after (b) nscoring the data. The DEM error data are derived from the elevation difference compared to the lidar DEM in ice- and snow-free areas. Outliers in the elevation difference and areas with slope > 20${}^{\circ }$ were also masked out. The spherical semivariogram model (red line) used in the SGSim and the parameters defining it ($c$, $c_{\mathrm{0}}$ and $r)$ are shown in (b). Graph (c) shows comparison between semivariograms for the deduced error (same as in (a)) and the difference compared to the lidar DEM in low-contrast areas within the glacier. (d)–(f) show the results of the SGSim for the 1946 glacier DEM. Images (d) and (e), respectively, show the mean and standard deviation of 1000 simulations at each 100 m $\times$ 100 m pixel. Graph (f) shows histogram (0.2 m bins) of the mean vertical bias values of the glacier DEM deduced from each simulation.

[Figure omitted. See PDF]

Each SGSsim, constrained by the input data and the spherical semivariogram model and calculated in resolution corresponding to 100 m $\times$ 100 m cell size, reveals possible errors in the measured glaciated area of the examined photogrammetric DEM. From each simulation the mean error of the glaciated area was calculated. From the 1000 simulations a histogram was derived and used to approximate a probability function of the likely bias in glaciated part of the DEM. Figure 4f shows the derived histogram for the 1946 DEM. Figure 4d and 4e, respectively, show the mean and the standard deviation of the derived error from 1000 simulations at each cell of the simulated area within the glacier. The latter reveals how the uncertainty in the derived error increases with distance from the input data. This should reach a maximum at a distance corresponding approximately to the range ($r)$ in the spherical semivariogram model, but all points on glacier in the 1946 DEM are at distance < $r$ from input data. The spatially varying mean error (Fig. 4d) could be used directly for correction of the photographic DEM, but instead we subtract the mean of the derived probability function to bias correct the area of interest in the photogrammetric DEM. Both approaches would lead to same result when deducing volume changes from the DEM differencing. The derived bias, $z$_bias, used to correct each DEM, and the corresponding 95 % upper ($z$_bias${}_{\mathrm{u}})$ and lower confidence limits ($z$_bias${}_{\mathrm{l}})$, is tabulated in Table 2. For comparison purposes the table also shows error bars derived by calculating analytically the expected variance (${\mathit{\sigma }}_{z_{\text{bias}}}^{\mathrm{2}}$) in the DEM error averaged over circular region corresponding to the size of Drangajökull, using a spherical semivariogram model (Rolstad et al., 2009), which fits the semivariogram without nscoring the error input data (Fig. 4a).

Finalizing the glacier DEMs

The photogrammetrically derived point clouds are typically much less dense for the snow-covered glacier surface than for bare ice or ground (see Supplement). The typical distance between points on the snow-covered glacier surface in the 1946 point cloud (the worst data set in terms of noise and point density) is $\sim$ 100 m, corresponding approximately to the resolution of the SGSim carried out. The point density is poorer for limited areas and in some regions gaps in the point clouds are caused by lack of contrast. Interpolating the elevation point clouds directly over long distances can be risky due to the spatial variability of the elevation. The spatial variability of the elevation changes derived from the difference between the point cloud and the lidar DEM is expected to be much lower (Cox and March, 2004). Therefore the bias-corrected difference was interpolated (Sect. 2.2) and added to the lidar DEM. The kriging function in Surfer 12 (^©Golden Software, Inc.) was used to interpolate the data by applying default linear variogram model and data search radius of 500 m. Even though the elevation changes compared to lidar are expected to be spatially smooth, interpolation over longer distances would reduce the reliability of the uncertainty assessment carried out for the photogrammetric DEMs. The different interpolation methods used within (kriging) and outside (linear) the glacier produces minor difference in the error statistics. For the 1946 bedrock data (after slope and outlier masking) $\mathit{\sigma }$ is 4.80 and 4.79 m for the linear and kriging methods, respectively, but 4.77 m derived directly from the point cloud difference compared to the full-resolution lidar DEM.

The resulting grids of elevation changes relative to lidar contained some larger gaps due to lack of contrast, cloud cover or incomplete coverage of aerial photographs for all data sets except the one of 2005 (Table 2). To complete the difference maps, two main interpolation methods were used. For relatively small gaps, spanning a short elevation range, kriging interpolation with data search radius > 500 m was applied using the derived elevation difference at the boundary of the data gap as input. For larger areas spanning significant elevation range we estimated a piecewise linear function for the elevation change as a function of the 2011 elevation (at 100 m elevation intervals) using the elevation difference between the point cloud and the lidar DEM as input (see Supplement). For data gaps covering an area at both the east and west side of the glacier the two different interpolations were carried out, one for the area west of the ice divides and another for the area east of it. In four cases neither of the above interpolation methods were considered applicable. The approaches adopted for each of these cases is described in the Supplement. The location of data gaps are shown in Fig. S1 and the interpolation method applied in each case is shown in the Supplement.

The uncertainties associated with interpolation of data gaps in the DEMs were approximated independently of the uncertainties of measured photogrammetric DEMs (Sect. 2.2). It is difficult to quantify these errors, but since these areas are generally small relative to the measured areas we adopted a generous approximation of the uncertainty roughly based on the scatter of the elevation change with altitude (point clouds compered to lidar DEMs). We assign three values of elevation uncertainty (95 % confidence level) to the interpolated areas, $\pm$7.5, $\pm$10 and $\pm$15 m, depending on the quality of the input data used for the interpolation and the applicability of the interpolation method (for further details see Supplement). The interpolated areas with the highest uncertainties were adopted for the lowermost part of Leirufjarðarjökull that was not covered in the 1975 survey flight (see Supplement). Also a relatively large area in the southernmost part of Drangajökull in 1946 where the interpolated area is poorly constrained by data. Clusters of data gaps are considered as a single area with assigned elevation uncertainty. We, however, assume that the error in one area is independent of the elevation error in other areas due to the distance between them.

Delineating glacier margins and nunataks

The glacier margin and nunataks at each time were delineated manually using the orthorectified aerial photographs at given time as well as the derived elevation difference compared to the lidar DEM. For 2011 the glacier outlines were drawn based on a shaded relief image of the 2011 DEM in maximum resolution and the intensity image of the lidar measurements. All glacier margins were delineated by the same person. The glacier margin was therefore interpreted in similar manner for all years, in areas where the outlines are uncertain. This working procedure minimizes variations in relative area changes of the ice cap. Due to numerous firn patches in the vicinity of Drangajökull, some of which are connected to the ice cap, it is actually a matter of definition whether these connected patches should be included as part of Drangajökull or not. We follow the approach of Jóhannesson et al. (2013) and exclude these patches. In a few areas the aerial photographs do not always reveal the glacier margin. This includes the southernmost part of Drangajökull in 1946. In this area the location of the glacier margin has been very stable since 1960. We therefore adopted the outermost glacier margin in the 1960–2011 data sets at each location as the 1946 margin in this area. Data used to approximate the location of the glacier margin in other areas where data are absent are described in the Supplement. The evolution of the glacier area is shown in Fig. 5. Also shown in Fig. 5 is the area of the eastern and western sections of the glacier, when Drangajökull is divided in two along the ice divides from north to south (see Fig. 6).

The relative area change of the entire Drangajökull ice cap, the western and eastern sections of the ice cap. The purple lines in Fig. 6 show the ice divides; they are used to define the east and west sections of the glacier. Labels give the glacier area in square kilometers at each epoch.

[Figure omitted. See PDF]

The average annual elevation change of Drangajökull during six intervals since 1946. Red indicates thinning and blue thickening.

[Figure omitted. See PDF]

Calculating volume changes

To derive the volume change, $\mathit{\delta }V(t_{\mathrm{s}},t_{\mathrm{f}})$ of the ice cap during a period $t_{\mathrm{s}}-t_{\mathrm{f}}$, the elevation difference DEM${}_{\mathrm{f}}$$-$DEM${}_{\mathrm{s}}$ (Fig. 6), was integrated over the area covered by glacier at either or both DEM dates. Continuous DEMs and glacier outlines had been completed for all years except for the year 1994, but this data set covered only $\sim$ $\mathrm{2}/\mathrm{3}$ of Drangajökull with the southernmost third of the ice cap missing. In order to estimate volume changes for this part of the glacier in the periods 1985–1994 and 1994–2005 the volume changes for the southernmost third of the glacier were plotted as a function of deduced volume changes in the other about two-thirds of the glacier for the periods 1960–1975, 1975–1985, 1985–2005 and 2005–2011 (Fig. 7). Linear fit describing the relation between the volume changes in the two areas estimated with least squares was used to estimate volume changes for the southern part of the glacier in the period 1985–1994 and 1994–2005. Errors in these volume change estimates were approximated using the 95 % confidence level of the linear fit (estimated in Grapher 10 ^©Golden Software, Inc.). Instead of approximating the position of the 1994 glacier margin, we only approximated the area covered by this part of the glacier. The volume change for the southernmost part of Drangajökull in the periods 1975–1985 was approximately the same as the estimated volume change in 1985–1994. We therefore extrapolate the 0.7 km${}^{\mathrm{2}}$ area increase of this glacier part of 1975–1985 to the period 1985–1994 to estimate the area of this glacier part in 1994.

The volume change of the southernmost of Drangajökull, which is missing in the 1994 DEM (Fig. 3), plotted as a function of the volume change in the northern part covered by the 1994 DEM, for the periods available (shown with black labels). The thick dashed line shows linear fit for the data points with the 95 % confidence area shown as light red. The red dots are the corresponding volume change estimates for the southern part in 1985–1994 and 1994–2005.

[Figure omitted. See PDF]

Seasonal correction of volume change between DEMs

The DEMs of Drangajökull were extracted from data acquired at different dates during the summer or the autumn (Table 1). Deriving mass balance records from DEM difference without taking this into account will skew the results, particularly if the acquisition time of the DEMs differs much from one DEM to another. Seasonal correction (sometimes referred to as date correction) has been applied and discussed in numerous studies (e.g., Krimmel, 1999; Cox and March, 2004; Cogley, 2009). In this study the derived volume change between DEMs ($\mathit{\delta }V(t_{\mathrm{s}},t_{\mathrm{f}})$, in Sect. 2.5) was seasonally corrected by compensating for the expected volume change of the ice cap from the acquisition date of each DEM until the end of the glaciological year (1 October). The end of the glaciological year was chosen because it makes comparison with both mass balance records and meteorological data easier and more eligible. This choice results in a larger magnitude of seasonal correction (and consequently larger uncertainty estimates) when compared to the average acquisition date of the DEMs. The seasonally corrected volume changes is given by $$\mathit{\delta }{V}^{\ast }\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)=\mathit{\delta }V\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)+\mathit{\delta }{V}_{S\text{\_cor}}\left(t_{\mathrm{s}}\right)-\mathit{\delta }{V}_{S\text{\_cor}}\left(t_{\mathrm{f}}\right).$$ The expected volume changes, $\mathit{\delta }{V}_{S\text{\_cor}}$ from the time of data acquisitions $t_{\mathrm{a}}$ until the end of the glaciological year $t_{\mathrm{end}}$ was estimated using a positive degree-day ($T_{+})$ model (e. g. Jóhannesson et al., 1995) with a constant degree-day factor (ddf) for the whole ice cap: $$\mathit{\delta }{V}_{S\text{\_cor}}={\displaystyle \frac{\mathrm{1}}{c_{\mathit{\delta }{V}_{S\text{\_cor}}}}}\cdot \text{ddf}\sum _{t_{\mathrm{a}}}^{t_{\text{end}}}\stackrel{\text{glacier}}{\int }{T}_{+}\left(t,,,x,,,y\right)\mathrm{d}A,$$ where $c_{\mathit{\delta }{V}_{S\text{\_cor}}}$ is the conversion factor from the glacier volume change during the period $t_{\mathrm{a}}-t_{\mathrm{end}}$ to the melt water draining from the ice cap in the same period. For seasonal volume correction of the DEMs in 1960 to 2011 we use daily grids of temperature at 2 m height above ground available back to 1949 (Crochet and Jóhannesson, 2011; unpublished data of the Icelandic Meteorological Office, for the year 2011). The grids were derived in two steps: (i) applying tension spline interpolation of measured temperature at meteorological stations corrected with fixed lapse rate to represent temperature at sea level; (ii) lapse rate adjustment (6.5 ${}^{\circ }$C km${}^{-\mathrm{1}})$ of interpolated temperature to compensate for the effects of topography. The temperature grids were in 1 km $\times$ 1 km cell size, but we linearly interpolated the grid in same resolution as the DEMs we are working with (20 m $\times$ 20 m cell size). Published values of ddf for Langjökull, Hofsjökull and Vatnajökull ice caps in Iceland (Jóhannesson et al., 2007; Guðmundsson et al., 2009) using comparable temperature data span the range from 4.45 mm w.e. ${}^{\circ }$C${}^{-\mathrm{1}}$ (minimum value for snow using lapse rate of 0.56 ${}^{\circ }$C km${}^{-\mathrm{1}}$; Jóhannesson et al., 2007) to 7.5 mm w.e. ${}^{\circ }$C${}^{-\mathrm{1}}$ (maximum value of firn/ice using lapse rate of 0.6 ${}^{\circ }$C km${}^{-\mathrm{1}}$; Guðmundsson et al., 2009). If these values had been obtained with lapse rate equal to the value, applied when creating the temperature grids used here (6.5 ${}^{\circ }$C km${}^{-\mathrm{1}})$, the resulting ddf would have been slightly higher. We therefore use ddf $=$ 6.5 $\pm$ 1.5 mm w.e. ${}^{\circ }$C${}^{-\mathrm{1}}$ instead of ddf $=$ 6.0 $\pm$ 1.5 mm w.e. ${}^{\circ }$C${}^{-\mathrm{1}}$ (the span of published values). Assuming that our conversion factor $c_{\mathit{\delta }{V}_{S\text{\_cor}}}=0.75\pm 0.1$ (where $c_{\mathit{\delta }{V}_{S\text{\_cor}}}=0.65$ corresponds to volume change mostly due to melting of snow and $c_{\mathit{\delta }{V}_{S\text{\_cor}}}=0.85$ corresponds to volume change mostly due melting of ice) is independent of ddf will result in seasonal corrections from Eq. (3) with 28 % uncertainty (95 % confidence level). The value of ddf is actually lower for snow than firn/ice; hence this assumption should lead to overestimation of the ratio ddf $/c_{\mathit{\delta }{V}_{S\text{\_cor}}}$ and consequently the uncertainty of $\mathit{\delta }{V}_{S\text{\_cor}}$ derived from Eq. (3). The seasonal volume correction of each DEM is given in Table 3. The aerial photographs used to produce the 1946 DEM were taken at the beginning of October before the start of winter snow fall. No seasonal correction was therefore required.

The average elevation change during periods defined by the DEMs before ($\mathit{\delta }h)$ and after $\left(\mathit{\delta }{h}^{\ast }\right)$ the seasonal correction, the seasonal correction $\left(\mathit{\delta }{h}_{S\mathit{\_}\text{cor}}\right)$ corresponding to DEM at time $t_{\mathrm{s}}$ and $t_{\mathrm{f}}$ (the correction at $t_{\mathrm{f}}$ is shown with minus sign since this correction term has minus in front of it in Eq. 2), the uncertainties (95 % confidence level) of seasonally corrected elevation change $\left(\mathrm{\Delta }\mathit{\delta }{h}^{\ast }\right)$ and the uncertainty contribution from the seasonal corrections $\left(\mathrm{\Delta }\mathit{\delta }{h}_{S\mathit{\_}\text{cor}}\right)$, DEM errors $\left(\mathrm{\Delta }\mathit{\delta }{h}_{\mathrm{m}}\right)$ and interpolation of data gaps $\left(\mathrm{\Delta }\mathit{\delta }{h}_i\right)$. All values were originally calculated in terms of volumes but are here averaged over the area $\overline{A}=\left(A\right(t_{\mathrm{f}})+A(t_{\mathrm{s}}\left)\right)/\mathrm{2}$.

$t_{\mathrm{s}}$	$t_{\mathrm{f}}$	Average	Average	Average	Average	Average	Average	Average	Average
		$\mathit{\delta }h$ (m)	$\mathit{\delta }{h}^{\ast }$ (m)	$\mathit{\delta }{h}_{S\text{\_cor}}\left(t_{\mathrm{s}}\right)$ (m)	$-\mathit{\delta }{h}_{S\text{\_cor}}\left(t_{\mathrm{f}}\right)$ (m)	$\mathrm{\Delta }\mathit{\delta }{h}^{\ast }$ (m)	$\mathrm{\Delta }\mathit{\delta }{h}_{S\text{\_cor}}$ (m)	$\mathrm{\Delta }\mathit{\delta }{h}_{\mathrm{m}}$ (m)	$\mathrm{\Delta }\mathit{\delta }{h}_i$ (m)
1946	1960	$-$7.36	$-$10.89	0	$-$3.53	2.73	0.99	1.28	2.19
1960	1975	$-$4.73	$-$1.27	3.69	$-$0.22	1.29	1.03	0.62	0.39
1975	1985	2.06	0.86	0.22	$-$1.42	0.95	0.40	0.54	0.62
1985	1994	2.15	2.84	1.42	$-$0.74	1.14	0.44	0.32	1.08
1994	2005	$-$7.11	$-$8.29	0.74	$-$1.92	1.17	0.58	0.26	0.96
2005	2011	$-$2.38	$-$3.23	1.94	$-$2.79	1.00	0.95	0.21	0

Deriving the geodetic mass balance and its uncertainty

The glacier-wide mass balance rate, $\dot{B}$ (the UNESCO IACS mass balance terminology (Cogley et al., 2011) is adopted) is estimated during the period $t_{\mathrm{s}}-t_{\mathrm{f}}$, using the equation: $$\dot{B}\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)={\displaystyle \frac{\mathit{\delta }{V}^{\ast }\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)}{\overline{A}\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)\cdot \mathit{\delta }t}}{c}_{\mathit{\delta }V},$$ where $\mathit{\delta }t=t_{\mathrm{f}}-t_{\mathrm{s}}$ and $\overline{A}(t_{\mathrm{s}},t_{\mathrm{f}})=(A\left(t_{\mathrm{s}}\right)+A\left(t_{\mathrm{f}}\right))/\mathrm{2}$ approximate the average area of the glacier during the period. It is reasonable to assume that the variables in Eq. (4) are independent of one another; hence the uncertainty in $\dot{B}$ can be approximated as $$\begin{array}{cc}& \mathrm{\Delta }\dot{B}\approx \sqrt{{\left(\mathrm{\Delta },\mathit{\delta },V^{\ast },{\displaystyle \frac{\partial \dot{B}}{\partial \mathit{\delta }{V}^{\ast }}}\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta },\overline{A},{\displaystyle \frac{\partial \dot{B}}{\partial \overline{A}}}\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta },c_{\mathit{\delta }V},{\displaystyle \frac{\partial \dot{B}}{\partial c_{\mathit{\delta }V}}}\right)}^{\mathrm{2}}}\\ & & ={\displaystyle \frac{\mathrm{1}}{\mathit{\delta }t}}\sqrt{{\left(\mathrm{\Delta },\mathit{\delta },V^{\ast },{\displaystyle \frac{c_{\mathit{\delta }V}}{\overline{A}}}\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta },\overline{A},{\displaystyle \frac{\mathit{\delta }{V}^{\ast }{c}_{\mathit{\delta }V}}{{\overline{A}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta },c_{\mathit{\delta }V},{\displaystyle \frac{\mathit{\delta }{V}^{\ast }}{\overline{A}}}\right)}^{\mathrm{2}}}.\end{array}$$ $\mathrm{\Delta }\overline{A}=\mathrm{4}$ km${}^{\mathrm{2}}$ is applied in all cases corresponding to $\sim$ 2.5 % of the glacier area, which is considered a generous estimate of the uncertainty in the glacier area for the given definition (Sect. 2.4). We used $c_{\mathit{\delta }V}=0.85\pm 0.06$ (Huss, 2013).

When estimating $\mathrm{\Delta }\mathit{\delta }{V}^{\ast }$ the error budget of $\mathit{\delta }{V}^{\ast }$ was examined. The error, $\mathit{\epsilon }$, of the seasonally corrected volume change, $\mathit{\delta }{V}^{\ast }(t_{\mathrm{s}},t_{\mathrm{f}})$, is the sum: $$\begin{array}{cc}& \mathit{\epsilon }\left\{\mathit{\delta },V^{\ast },\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)\right\}=\mathit{\epsilon }\left\{V_{\mathrm{m}},\left(t_{\mathrm{s}}\right)\right\}+\mathit{\epsilon }\left\{V_i,\left(t_{\mathrm{s}}\right)\right\}\\ & +\mathit{\epsilon }\left\{\mathit{\delta },V_{S\mathit{\_}\text{cor}},\left(t_{\mathrm{s}}\right)\right\}+\mathit{\epsilon }\left\{V_{\mathrm{m}},\left(,t_{\mathrm{f}},\right)\right\}+\mathit{\epsilon }\left\{V_i,\left(t_{\mathrm{f}}\right)\right\}\\ & & +\mathit{\epsilon }\left\{\mathit{\delta },V_{S\mathit{\_}\text{cor}},\left(t_{\mathrm{f}}\right)\right\},\end{array}$$ where the error in the measured volume at time $t$ is $$\mathit{\epsilon }\left\{V_{\mathrm{m}},\left(,t,\right)\right\}=A_{\mathrm{m}}\left(t\right)\cdot \overline{\mathit{\epsilon }}\left\{h,\left(t\right)\right\},$$ where $A_{\mathrm{m}}$ is the area of measured DEM within the glacier and $\overline{\mathit{\epsilon }}\mathit{\left\{}h\mathit{\right\}}$ the mean error of the glaciated area. The error in volume for the interpolated glacier sections lacking measurement (Sect. 2.3) is $$\mathit{\epsilon }\left\{V_i,\left(,t,\right)\right\}=\sum _{j=\mathrm{1}}^N{A}_j\left(t\right)\cdot \overline{\mathit{\epsilon }}\left\{h_j,\left(t\right)\right\},$$ where $A_j$ is the area of the interpolated section,$j$, and $\overline{\mathit{\epsilon }}\mathit{\left\{}{h}_j\mathit{\right\}}$ is the corresponding mean elevation error. Assuming that the individual errors contributing to Eqs. (6) and (8) are independent of one another the probability function of the error in $\mathit{\delta }{V}^{\ast }\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)$ is given by the multiple convolutions: $$\begin{array}{cc}& f_{\mathit{\epsilon }\left\{\mathit{\delta },V^{\ast },\left(t_{\mathrm{s}},,,t_{\mathrm{f}}\right)\right\}}=f_{\mathit{\epsilon }\left\{V_{\mathrm{m}},\left(t_{\mathrm{s}}\right)\right\}}\times f_{\mathit{\epsilon }\left\{V_i,\left(t_{\mathrm{s}}\right)\right\}}\times f_{\mathit{\epsilon }\left\{\partial ,V_{S\text{\_cor}},\left(t_{\mathrm{s}}\right)\right\}}\\ & & \times f_{\mathit{\epsilon }\left\{V_{\mathrm{m}},\left(t_{\mathrm{f}}\right)\right\}}\times f_{\mathit{\epsilon }\left\{V_i,\left(t_{\mathrm{f}}\right)\right\}}\times f_{\mathit{\epsilon }\left\{\mathit{\delta },V_{S\text{\_cor}},\left(t_{\mathrm{f}}\right)\right\}}.\end{array}$$ The probability function $f_{\mathit{\epsilon }\mathit{\left\{}{V}_{\mathrm{m}}\right(t\left)\mathit{\right\}}}$ was derived directly from Eq. (10) and by approximating $f_{\overline{\mathit{\epsilon }}\mathit{\left\{}h\right(t\left)\mathit{\right\}}}$ using the histogram of the corresponding elevation bias correction (Sect. 2.2) minus its mean. All other errors are assumed to be normally distributed with zero mean, hence $$f\left(\mathit{\epsilon }\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{\epsilon }}\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{{\mathit{\epsilon }}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}}}.$$ The probability distribution $f_{\mathit{\epsilon }\left\{V_i\right\}}$ is thus also a normal distribution with $${\mathit{\sigma }}_{\mathit{\epsilon }\left\{V_i\right\}}=\sqrt{\sum _{j=\mathrm{1}}^N{\left(A_j,\cdot ,{\mathit{\sigma }}_{\overline{\mathit{\epsilon }}\left\{h_j\right\}}\right)}^{\mathrm{2}}}.$$ The uncertainty in the volume change $\mathrm{\Delta }\mathit{\delta }V$ (95 % confidence level) was now derived from the probability distribution given by Eq. (9). Table 3 shows the 95 % confidence level of $f_{\mathit{\epsilon }\mathit{\left\{}{V}_{\mathrm{m}}\mathit{\right\}}},f_{\mathit{\epsilon }\mathit{\left\{}{V}_i\mathit{\right\}}}$ and $f_{\mathit{\epsilon }\mathit{\left\{}\mathit{\delta }{V}_{S\text{\_cor}}\mathit{\right\}}}$ for each year of acquisition, revealing the main source of error in the derived volume changes.

Results

Bias corrections and uncertainty estimates deduced from the DEM errors

Table 2 gives values of several error estimation parameters for the photogrammetric DEMs deduced by comparison with the 2011 lidar DEM in ice- and snow-free areas. Some of these parameters can be used both to correct the DEMs and to estimate the uncertainty of geodetic mass balance results. In some cases significant difference is observed between the mean DEM error, commonly used to correct for bias (zero order correction) of the DEM (e.g., Guðmundsson et al. 2011), and the bias derived from the SGSim. The greatest difference is for the 1946 DEM, which after removal of outliers and steep slopes the ice- and snow-free part of it has a mean error of $-$0.86 m whereas the SGSim results in bias of 1.66 m. The difference would presumably be lower if we would only calculate the mean error using areas within certain distance from the glacier margin, but it is not straightforward to select this distance without using geostatistical approaches.

The parameters in Table 2 that can be used to estimate the uncertainty of geodetic mass balance show even more diversity. The crudest parameter would be the standard deviation of the DEM error derived from ice- and snow-free areas. Standard deviation is commonly interpreted as 68 % confidence level assuming normal error distribution and should therefore be multiplied by 1.96 to obtain 95 % confidence level as derived for the other two approaches shown in Table 2. This interpretation of the standard deviation as uncertainty proxy of the volume change implies the assumption that the DEM errors at different locations within the glacier are totally correlated (Rolstad et al., 2009). Since the confidence level of geodetic mass balance results is typically not mentioned in studies using the standard deviation as their uncertainty proxy, the conversion of the standard deviation to 95 % confidence level is omitted in Table 2. The values of standard deviation for the ice-free DEMs are 5–45 % lower after removal of outlier and steep slopes. The lower standard deviation values are, however, still by far higher than the uncertainty (95 % confidence level) of the bias correction derived with SGSim. The SGSim results in uncertainty between 0.21 m (in 2005) and 1.58 m (in 1946). The SGSim uncertainties correspond to 24–46 % of the standard deviation (after slope and outlier removal). If we exclude the three DEMs from 1960, covering only about one-third of Drangajökull each, the range is 24–33 %. The SGSim uncertainties correspond to 27–80 % of the uncertainties derived with method described by Rolstad et al. (2009) and the percentage seems to depend strongly on the range of the spherical semivariogram model used in both calculations (Fig. 8).

The ratio between uncertainties (95 % confidence level) from the methods demonstrated in this work and the method demonstrated by Rolstad et al. (2009) as a function of the range, $r$, in the deduced spherical semivariogram model. The DEM epoch corresponding to each point is shown with a label.

[Figure omitted. See PDF]

The glacier-wide mass balance rate ($\dot{B})$ during six different periods since 1946, calculated for the entire ice cap (a), and split into the western and the eastern sections of Drangajökull (b) using the ice divides shown as purple lines in Fig. 6. The dotted line in (a)–(b) represents 95 % confidence level. Graph (c) shows the average summer temperature at the meteorological stations Gjögur, since 1949, and Æðey since 1954 (see Fig. 1b for locations). Close circles indicate data from manned station, open circles from automatic station. The dotted lines show the average summer temperature at each location filtered with 11-year triangular weighted running average.

[Figure omitted. See PDF]

DEM seasonal corrections and contribution of different error sources to the geodetic mass balance

The effects of seasonal correction and the estimated contribution of each type of error to the total volume change are summarized in Table 3. The importance of seasonal correction for Drangajökull is clearly revealed, particularly for the first two periods, 1946–1960 and 1960–1975, due to the early acquisition of the 1960 aerial photographs. The sum of the two seasonal corrections for these periods corresponds to a larger value than the derived total uncertainty of $\mathit{\delta }{V}^{\ast }$. The correction effectively increases the difference in $\dot{B}$ between the periods by 0.42 m w.e. a${}^{-\mathrm{1}}$ ($\sim$ 0.21 m w.e. a${}^{-\mathrm{1}}$ absolute change for each period). With the inferred correction the period 1946–1960 is the period of highest mass loss rate, along with 1994–2005, whereas the period 1960–1975 differs insignificantly from equilibrium (Fig. 9). For other periods the net seasonal correction changed the derived $\dot{B}$ by 0.06–0.12 m w.e. a${}^{-\mathrm{1}}$.

The main source of uncertainties is different from one period to another, but in no case is the highest contribution from the estimated uncertainty of the DEM elevation ($\mathrm{\Delta }\mathit{\delta }{h}_{\mathrm{m}}$ in Table 3). For periods when volume change is based on the 1946 or 1994 DEM, we have relative high uncertainties due to interpolations of large gaps in the derived DEMs (Fig. 3 and Table 2). The derived value of $\mathit{\delta }{V}^{\ast }$ for the period 2005–2011, obtained from the best two DEMs in terms of accuracy and coverage, has significant uncertainty due to large seasonal correction for both DEMs. The 2005 and 2011 data were acquired in late July, and the summer remainder for both years was relatively warm. The sum of seasonal corrections (which have opposite signs) is actually smaller than the uncertainty related to the seasonal corrections for the period 2005–2011.

The uncertainty percentage of $\mathit{\delta }{V}^{\ast }$ is typically significantly higher than the uncertainty percentage of A (2.5 %) and $c_{\mathit{\delta }V}$ ($\sim$ 7 %). Uncertainty of the derived $\dot{B}$ (Fig. 9) produced by the uncertainty of the latter two variables is therefore generally minor compared to the uncertainty contribution of $\mathit{\delta }{V}^{\ast }$.

The geodetic mass balance of Drangajökull

Figure 9 shows the derived $\dot{B}$ for Drangajökull during six intervals since 1946. During the period 1946–1960 relatively high mass loss rates of $\dot{B}=-0.66\pm 0.17$ m w.e. a${}^{-\mathrm{1}}$ are estimated. The glacier was near equilibrium in 1960–1985 with $\dot{B}=-0.07\pm 0.07$ m w.e. a${}^{-\mathrm{1}}$ and $\dot{B}=0.07\pm 0.08$ m w.e. a${}^{-\mathrm{1}}$ in 1960–1975 and 1975–1985, respectively. The mass balance was significantly positive in 1985–1994 with $\dot{B}=0.26\pm 0.11$ m w.e. a${}^{-\mathrm{1}}$. In the period 1994–2005 again, as in the mid-century, there is high rate of mass loss with $\dot{B}=-0.64\pm 0.10$ m w.e. a${}^{-\mathrm{1}}$ and then slightly less negative mass balance rate in 2005–2011, with $\dot{B}=-0.46\pm 0.15$ m w.e. a${}^{-\mathrm{1}}$. The glacier-wide mass balance rate for the entire period 1946–2011 is $\dot{B}=-0.26\pm 0.04$ m w.e. a${}^{-\mathrm{1}}$. In the same period Drangajökull was reduced in area by $\sim$ 11 % from 161 to 144 km${}^{\mathrm{2}}$ (Fig. 5).

The middle panel of Fig. 9 shows $\dot{B}$ for the western and eastern half of Drangajökull ice cap, as defined by the ice divides from north to south shown in Fig. 6. The results are derived in the same manner as the result for the entire glacier, where the steps taken to correct for bias of the DEM, derive seasonal correction and derive uncertainties were carried out focusing specifically on either the western or the eastern part. The bias correction of each half may vary up to a few decimeters from the correction of the entire ice cap and the uncertainty limits of the bias correction is generally slightly higher.

Figure 9 shows different evolution of the west and east glaciers. Both parts suffered significantly negative mass balance rate in 1946–1960 and 1994–2011. The period in between was significantly negative on the east side, apart from the period 1985–1994, when the upper 95 % confidence level is slightly above 0, whereas the western part had $\dot{B}$ near 0 in 1960–1975 and significantly positive mass balance rate with $\dot{B}=0.23\pm 0.10$ m w.e. a${}^{-\mathrm{1}}$ and $\dot{B}=0.52\pm 0.15$ m w.e. a${}^{-\mathrm{1}}$ in 1975–1985 and 1985–1994, respectively. A mean mass balance rate of $\dot{B}=-0.16\pm 0.05$ m w.e. a${}^{-\mathrm{1}}$ is derived for the period 1946–2011 on the western part. The mass loss rate is on average $\sim$ 3-fold higher for the eastern part with $\dot{B}=-0.41\pm 0.04$ m w.e. a${}^{-\mathrm{1}}$. This is also reflected in the area change but in 1946–2011 the eastern part decreased in area 21 %, while the western part shrank only by 3 % (Fig. 5).

Discussion

The high precision of the geodetic mass balance results presented can be primarily explained by (i) the use of the high-resolution and high-accuracy lidar DEM to extract evenly distributed GCPs for constraining the orientation of photogrammetric DEMs (obtaining equivalent distribution of GCPs in the field was not possible within the financial frame of this study) and (ii) the thorough uncertainty assessment of the results where the lidar data from ice- and snow-free areas are also key data since they enable assessment of geostatistical parameters of the photogrammetric DEMs. Both (i) and (ii) highlight the need for high-resolution and high-accuracy DEMs from the present in areas of interest to conduct studies of geodetic mass balance using aerial photographs from the past. The third important use of the lidar data in this study is the creation of DEMs from the photogrammetric point clouds within the glacier. Rather than interpolating the elevation point clouds directly we interpolate the difference between the point cloud and lidar DEM (much less high-frequency variability, and the difference is a smoother surface; Cox and March, 2004) and add the interpolated product to the lidar DEM. This results in more accurate DEMs in areas where the density of the photogrammetric point clouds is low.

Other state-of-the-art high-resolution elevation data sets obtained with airborne or spaceborne sensors are also suitable to replace the lidar data in the work procedure described here. This probably includes Worldview and Pléiades high-resolution stereo images, allowing extraction of DEM with < 5 m cell dimensions and orthorectified photographs with < 1 m $\times$ 1 m cell size (e.g., Berthier et al., 2014; Howat et al., 2015). Part of the work procedure described has already been carried out using such satellite data as replacement for the airborne lidar. In a recent study by Papasodoro et al. (2015), Pléiades data were used to collect GCPs for constraining DEMs from aerial photographs. Even though the absolute accuracy of data from spaceborne sensors does not match data from airborne lidar, it does not make the satellite data inadequate. Each photogrammetric DEM from the past can be fixed into the reference frame of the high-resolution DEM through the extraction and usage of GCPs and implementation of proposed bias correction. The relative elevation change between DEMs should therefore be fairly accurate despite lower absolute accuracy of the DEMs, and shifts and tilts of the reference frame cancel out in DEM differencing.

In this study, the derived bias correction of the glaciated DEM section and the uncertainty of volume changes related to DEM errors are obtained from the probability distribution calculated by using SGSim. The bias correction corresponds to the probabilistic mean of the average error within the glacier. As shown in Table 2 the difference between the mean error in snow- and ice-free areas and the bias derived from the SGSim (the estimated probabilistic mean of the glacier DEM error) was up to 2.5 m (in 1946). This difference would presumably be lower if we would only calculate the mean error using areas within certain distance from the glacier margin, but it is not straightforward to select this distance without using some geostatistical approaches. The relation is also not obvious between the probabilistic mean of an average DEM error within the glacier and higher-order corrections of a glacier DEM obtained with least square fit (or similar) using deduced DEM errors in ice- and snow-free areas. If the average correction does not correspond to the probabilistic mean, the results of geodetic mass balance will be incorrectly centered even if the width of the error bars is realistic.

When comparing different proxies used for estimating the uncertainty of DEM difference derived volume change, it is no surprise that using the standard deviation of the DEM error in snow- and ice-free areas leads to great overestimate of the uncertainty (Table 2). This has been shown before by Rolstad et al. (2009). Other estimators that ignore information of the spatial dependency of the DEM errors, such as the NMAD value (Höhle and Höhle, 2009), should also be considered as incomplete for this purpose.

The difference in uncertainty estimates between the method described here and the method of Rolstad et al. (2009) is especially noteworthy (Table 2 and Fig. 8). Rolstad et al. (2009) provided a simple and logical method to estimate the uncertainty of derived volume change. The DEM errors (or difference) in ice- and snow-free areas are used to calculate a semivariogram that constrains a spherical semivariogram model. From the spherical semivariogram model alone the expected variance of the DEM error (${\mathit{\sigma }}_{z\text{\_bias}}^{\mathrm{2}})$ averaged over circular region corresponding to the size of the glacier is calculated analytically. The method compensates for the spatial dependency of the DEM error at different locations within the glacier. The method does not, however, take into the account how the DEM error within the glacier depends on the DEM errors outside the glacier, unlike the method proposed here utilizing SGSim. This is most likely the explanation why the ratio between the two uncertainty estimates ($\mathrm{\Delta }z\mathit{\_}{\text{bias}}_{\text{SGSim}}/\mathrm{\Delta }z\mathit{\_}{\text{bias}}_{\text{Rols}}$) appears to be strongly dependent on the range, $r$, in the spherical semivariogram model, which is common for both approaches (Fig. 8). If $r$ is small compared the size of the glacier, meaning that large proportion of the glacier has DEM error independent of DEM error outside the glacier, the uncertainty derived by SGSim is only slightly smaller than the uncertainty derived analytically from the spherical semivariogram model alone. If $r$ is large, however, meaning that a large proportion or even the entire glacier has DEM error dependent on the DEM errors outside the glacier, the SGSim results in much lower uncertainty. This interpretation implies that the method of Rolstad et al. (2009) gives a good approximation of the uncertainty when most of the glaciated area is at distance > $r$ from ice- and snow-free areas providing measurements of the DEM errors, but it can otherwise result in great overestimation of the uncertainty in the derived volume change. The main disadvantage of the SGSim approach compared to the approach of Rolstad et al. (2009) is that it is more time consuming. The tool applied here (WinGSlib) also has problems with data sets larger than those worked with in this study. New tools enabling the SGSim approach for large data sets should be developed in order to facilitate the usage of this methodology.

Our study emphasizes the importance of including seasonal correction of DEMs for glacier with high mass turnover to avoid wrong interpretation of derived volume change. The most extreme case is the negative volume change derived from the difference between the 1960 and 1975 DEMs. The seasonal correction results in about three-quarters of this negative volume change being effectively transferred in to the period 1946–1960 due to large seasonal correction of the 1960 DEM resulting from relatively early acquisition of the aerial photographs (Table 1). The seasonally corrected volume change revealing the volume change between the start of different glaciological year obviously has higher uncertainty than the uncorrected volume change. We, however, consider this trade-off important for easy comparison with other data records, including meteorological data and in situ mass balance measurements. The uncertainty due to the seasonal correction as well as the uncertainty related to the interpolation of the data gaps should be considered as cautious estimates of the 95 % confidence level of the error associated with these two error sources. Effort should be made to constrain these uncertainties further, which could narrow the uncertainty estimates of this study and other similar even further, but it is beyond the scope of this paper.

The presented geodetic mass balance records indicate slower volume decrease for Drangajökull ice cap since the 1940s than for most other glaciers in Iceland with geodetic mass balance record extending back to that period. While we observe $\dot{B}=-0.26\pm 0.04$ m w.e. a${}^{-\mathrm{1}}$ for Drangajökull in the period 1946–2011 the corresponding values for Langjökull ice cap in 1945–2011 is $\dot{B}\approx -0.5$ m w.e. a${}^{-\mathrm{1}}$ (Pálsson et al., 2012, with extension from traditional mass measurements in 2004–2011 from Björnsson et al. (2013) until 2010 and unpublished Institute of Earth Sciences (IES) data for 2011). Two outlets of S-Vatnajökull, Kvíárjökull and Skaftárjökull have similar rate of mass decrease in 1945–2010 or $\dot{B}\approx -0.25$ m w.e. a${}^{-\mathrm{1}}$ (Hannesdóttir et al., 2015). Other outlets of S-Vatnajökull ice cap show ${\overline{B}}_n$ between $-$0.3 and $-$0.8 m w.e. a${}^{-\mathrm{1}}$ in 1945–2010 (Hannesdóttir et al., 2015; Aðalgeirsdóttir et al., 2011). For the relatively warm period in 1994–2011 we obtain $\dot{B}=-0.58\pm 0.08$ m w.e. a${}^{-\mathrm{1}}$, which is in good agreement with the study of Jóhannesson et al. (2013), which indicated $\dot{B}\approx -0.5$ m w.e. a${}^{-\mathrm{1}}$ for Drangajökull ice cap in the period 1996–2011. Comparisons of Drangajökull mass balance in 1994–2011, with results from traditional in situ mass balance measurements from Langjökull (in 1996–2011) and Vatnajökull ice caps, show that the reduction rate has been $\sim$ 140 % faster on Langjökull ($\dot{B}\approx -1.4$ m w.e. a${}^{-\mathrm{1}}$, from Björnsson et al. (2013) until 2010 and unpublished IES data for 2011) and $\sim$ 20 % faster on Vatnajökull ($\dot{B}\approx -0.7$ m w.e. a${}^{-\mathrm{1}}$, from Björnsson et al. (2013) until 2010 and unpublished IES data for 2011).

The difference in the geodetic mass balance results between the east and west part of Drangajökull highlights how difficult it is to extrapolate mass balance records from one glacier to another, even over short distances. The results, showing $\sim$ 3 times more negative mass balance rate for the eastern part of Drangajökull than the western part for the entire period 1946–2011, are not reflected in changing spatial trends of summer temperature during the period. The summer temperature measured east of Drangajökull is typically $\sim$ 1 ${}^{\circ }$C lower than revealed by measurements west of Drangajökull (Fig. 9c) and this is rather consistent throughout the survey period. Daily precipitation maps (1 km $\times$ 1 km cell size) in 1958–2006 deduced from ERA-40 (Uppala et al., 2005) by dynamic downscaling with linear model of orographic precipitation (an update of Crochet et al. (2007) described in Jóhannesson et al., 2007) do not indicate a strong trend in winter precipitation from east to west. The modeled winter precipitation may, however, not be representative for winter accumulation due to excess lee drying in the modeled precipitation or transport of snow from east to west by snow drift; the most common wind direction on Drangajökull is from the northeast. Most of the precipitation also falls on the glacier when the wind blows from the northeast. Ongoing geodetic mass balance studies of Drangajökull on a seasonal timescale may reveal further answers.

The geodetic mass balance record on Drangajökull ice cap is the first record revealing glacier volume change in Iceland on a decadal timescale of the past $\sim$ 70 years. Other records spanning this period have a coarser resolution particularly over the period 1945–1985, which is typically assigned a single mass balance value (e.g., Pálsson et al., 2012; Hannesdóttir et al., 2015). However, accurate and detailed studies pertaining to this period are of particular interest as they may reveal how the Icelandic glaciers responded to the change from a relatively warm climate in 1925–1965 to a significantly colder climate in 1965–1990 and subsequently to a warming with a short setback around 1995 (cf. Figs. 2.6 and 3.1 in Björnsson et al., 2008). We consider this study the first step in filling this gap in our knowledge. The key data set to continue this work is the archive of aerial photographs at the National Land Survey of Iceland, covering the Icelandic glaciers in the 1940s–1990s. Similar archives covering other glaciated parts of the world should be fully utilized using the new processing techniques and recent and future availability of high-resolution DEMs of the present state of the glaciated areas and its vicinity.

Conclusions

This paper highlights the opportunities that new high-resolution DEMs are opening to improve the procedure carried out to obtain geodetic mass balance records. Such DEMs are key data in three aspects of this study: (a) extracting GCPs from recent airborne lidar DEM to constrain photogrammetric DEMs at six different epochs; (b) interpolating over glacier surface the elevation difference of derived photogrammetric point cloud relative to the lidar DEM; (c) applying new geostatistical approaches based on comparison with the lidar data to estimate simultaneously a bias correction for the glacier DEMs along with its 95 % confidence level. The latter reveals the uncertainty associated with DEM errors in geodetic mass balance records.

The new geostatistical method applies SGSim using the DEM errors in ice- and snow-free areas and a spherical semivariogram model constrained by the DEM errors as input data. The resulting bias correction may differ considerably (in our case up to 2.5 m in 1946) from the simple approach of applying bias correction using the mean DEM error outside the glacier. The resulting uncertainty of the DEM (95 % confidence level) was typically estimated 20–35 % of the standard deviation derived from the DEM errors in ice- and snow-free areas after outliers and high slopes were masked out. The uncertainty contribution from DEM errors obtained with SGSim was 25–80 % of the uncertainty estimate obtained with the geostatistical method of Rolstad et al. (2009). We argue that methods typically carried out in uncertainty assessments of geodetic mass balance generally overestimate the uncertainty related to DEM errors, while the geostatistical approach described here results in more realistic uncertainty estimates.

This study also reveals the importance of seasonal corrections of geodetic mass balance for glaciers with high annual turnover; Drangajökull is a good example. The highest correction in our study was $-$3.5 m (in 1960), which corresponds to approx. three-quarters of the average elevation change between the 1960 and the 1975 DEMs.

During the whole period 1946–2011 we obtain $\dot{B}=-0.26\pm 0.04$ m w.e. a${}^{-\mathrm{1}}$ for the entirety of Drangajökull. This is among the lowest retreat rate reported for glaciers in Iceland spanning approximately this period. Only two outlet glaciers in S-Vatnajökull have similar reported retreat rates. When calculating this for the western and eastern half of Drangajökull specifically we obtain $\dot{B}=-0.16\pm 0.05$ m w.e. a${}^{-\mathrm{1}}$ and $\dot{B}=-0.41\pm 0.04$ m w.e. a${}^{-\mathrm{1}}$, respectively. This difference between the eastern and western parts of the glacier varies significantly during the survey period and does not seem to be related to relative changes in summer temperature. This great difference between east and west shows how difficult it is to extrapolate mass balance records from one glacier to another even over short distances. No glacier unit in Iceland has been reported as close to equilibrium on average since the 1940s as the western part of Drangajökull ice cap.

The Supplement related to this article is available online at doi:10.5194/tc-10-159-2016-supplement.

The writing of this paper and the research it describes was mostly carried out by the first two authors of this paper, with input from the other three co-authors. All photogrammetric processing and revision of the resulting point clouds was carried out by J. M. C. Belart. The interpolation of the point cloud differences compared to the lidar DEM and construction of glacier DEM based on that, interpolation of data gaps, delineation of glacier margin, seasonal correction of the volume change, the construction of the presented mass balance records and associated error analysis was carried out by E. Magnússon based on fruitful discussions with J. M. C. Belart and F. Pálsson. All figures in this paper were made by E. Magnússon and J. M. C. Belart as well as tables. P. Crochet and H. Ágústsson contributed to the handling and interpretation of the meteorological data.

Acknowledgements

This work was carried out within SVALI funded by the Nordic Top-level Research Initiative (TRI) and is SVALI publication number 70. It was also financially supported by alpS GmbH. This work is a contribution to the Rannís grant of excellence project, ANATILS. We thank the National Land Survey of Iceland and Loftmyndir ehf. for acquisition and scanning of the aerial photographs. This study used the recent lidar mapping of the glaciers in Iceland that was funded by the Icelandic Research Fund, the Landsvirkjun Research Fund, the Icelandic Road Administration, the Reykjavík Energy Environmental and Energy Research Fund, the Klima- og Luftgruppen (KoL) research fund of the Nordic Council of Ministers, the Vatnajökull National Park, the organization Friends of Vatnajökull, the National Land Survey of Iceland and the Icelandic Meteorological Office. We thank Cristopher Nuth and an anonymous reviewer for constructive reviews, as well as Etienne Berthier for helpful comments on the manuscript. We thank Alexander H. Jarosch for fruitful discussions on the subject of this paper. Auður Agla Óladóttir is thanked for introducing the tools and methodology used in the error analysis of this study to E. Magnússon. Edited by: I. M. Howat