Winter mass balance of Drangajökull ice cap ( NW Iceland ) derived from satellite sub-meter stereo images

Sub-meter resolution, stereoscopic satellite images allow for the generation of accurate and high-resolution Digital Elevation Models (DEMs) over glaciers and ice caps. Here, repeated stereo images from Pléiades and WorldView2 (WV2) of Drangajökull ice cap (NW-Iceland) are combined with in situ estimates of snow density and densification of firn and fresh snow to provide the first estimates of the glacier-wide geodetic winter mass balance to be obtained from satellite imagery. 15 Statistics in snowand ice-free areas reveal similar vertical relative accuracy (<0.5 m) with and without ground control points, demonstrating the capability for measuring seasonal snow accumulation. The winter (14 October to 22 May) mass balance of Drangajökull was 3.33 ± 0.23 mw.e. with ~60% of the accumulation occurring by February, in good agreement with nearby ground observations. The repeat DEMs yield on average 22% less elevation change than the length of 8 winter snow cores due to (1) the difference in time between in situ and satellite observations, (2) firn densification and (3) elevation 20 changes due to ice dynamics. The contributions of these 3 factors were of similar magnitude. This study demonstrates that seasonal geodetic mass balance can, in many areas, be measured from sub-meter resolution satellite stereo images.


Introduction
Monitoring glacier changes improves understanding of the close connection between glacier mass balance and climate (Vaughan et al., 2013).Glacier monitoring is based on in situ and remote sensing measurements and has confirmed the strong sensitivity of glaciers to climate change.Monitoring has provided evidence for the continuous retreat and mass loss currently taking place in most glaciated regions on Earth (Vaughan et al., 2013;Zemp et al., 2015).
Observations of mass balance provide a valuable shortterm overview of the glacier's mass budget and its implications for water storage, runoff and regional climate (e.g., Huss et al., 2008;Radić and Hock, 2014).In addition, these observations can reveal trends and patterns in glacier mass evolution and are commonly used in glacier modeling (e.g., Huss et al., 2008;Aðalgeirsdóttir et al., 2011).Seasonal records of glacier mass changes, however, are sparse, and many glaciated areas in the world are not currently monitored due to high cost and logistical challenges (Ohmura, 2011).
The most widely used method for measuring winter mass balance is the glaciological method (i.e., snow probing, snow pits and/or shallow cores).With adequate spatial sampling, this method can be used to estimate glacier-wide mass balance with errors of 0.1 to 0.3 m water equivalent (m w.e.; Fountain and Vecchia, 1999;Ohmura, 2011).Remote-sensing-based methods, such as repeated airborne surveys (Machguth et al., 2006;Sold et al., 2013;Helfricht et al., 2014) or unmanned aerial vehicles surveys (Bühler et al., 2016;De Michele et al., 2016), are occasionally used for measuring snow accumulation.These methods allow for the creation of highly accurate and detailed digital elevation models (DEMs) that are compared when measuring changes in elevation and volume due to snow accumulation.
Satellite stereo images with sub-meter resolution (e.g., from WorldView or Pléiades with nearly global coverage) are available for the creation of accurate and detailed DEMs.The high spatial and radiometric resolutions of these images allow for the statistical correlation of features on lowcontrast surfaces, including ice, snow and shadowed terrain (e.g., Berthier et al., 2014;Holzer et al., 2015;Willis et al., 2015;Shean et al., 2016).The DEMs obtained from these sensors have been tested and assessed in numerous studies, reporting relative DEM accuracy ranging from 0.2 to 1 m (Berthier et al., 2014;Lacroix et al., 2015;Noh and Howat, 2015;Willis et al., 2015;Shean et al., 2016).This accuracy indicates high potential for the usage of these sensors in measuring changes over short time intervals for glaciers with relatively high mass balance amplitude (half of the difference between winter and summer mass balance, Cogley et al., 2011).For example, sequential Pléiades DEMs have recently been successfully used for measuring snow thickness in mountainous areas (Marti et al., 2016).
In this paper, we evaluate the capabilities of Pléiades and WorldView2 (WV2) DEMs for measuring winter mass balance over an Icelandic ice cap.A processing chain is developed for constructing co-registered DEMs from submeter resolution optical stereo images.Co-registration is performed without external reference data, enabling application to remote glaciated areas where such data is lacking.Calculation of geodetic winter mass balance is constrained with in situ density measurements as well as simple firn and snow densification models.Finally, we validate our remote sensing results with in situ measurements of snow thickness.
2 Study site and data

Drangajökull ice cap
Approximately 11 000 km 2 of Iceland is covered by glaciers (Björnsson and Pálsson, 2008).Glaciological mass balance studies have been conducted on the three largest ice caps: Vatnajökull (since 1991, Björnsson et al., 2013), Langjökull (since 1997, Pálsson et al., 2012) and Hofsjökull (since 1988, Jóhannesson et al., 2013) (Fig. 1).Field campaigns are carried out twice per year to record the winter and summer mass balance at selected survey sites (Björnsson and Pálsson, 2008;Björnsson et al., 2013), and the measurements reveal typical mass balance amplitude of ∼ 1.5-3 m w.e.(Aðalgeirsdóttir et al., 2011;Pálsson et al., 2012;Björnsson et al., 2013) and even higher in some other glaciated areas such as Mýrdalsjökull and Öraefajökull ice caps (south Iceland) where limited mass balance surveys in the accumulation area have shown winter accumulation of 5-7 m w.e.(Guðmundsson, 2000;Ágústsson et al., 2013).These measurements have improved understanding of the impacts of climate change on glacier mass balance in the North Atlantic and have provided glacial runoff estimates, which are important for water resource management in Iceland.
The study area, Drangajökull ice cap, is located in NW Iceland (Fig. 1) between ∼ 60 and ∼ 900 m a.s.l. and has a total area of 143 km 2 (in 2014).Due to its distance from the Irminger Current, its climate is substantially different from other Icelandic glaciers near the south coast or in the central part of the island (Jóhannesson et al., 2013;Harning et al., 2016a, b).Geodetic observations have revealed that the average glacier-wide mass balance of Drangajökull during the period 1946-2011 was moderately negative (−0.26 ± 0.04 m w.e. a −1 ; Magnússon et al., 2016a).The same observations revealed a striking difference in the mass balance between the western and eastern sides of the ice cap during this period, −0.16 ± 0.05 and −0.41 ± 0.04 m w.e. a −1 , respectively.The spatial distribution of the winter snow accumulation is a likely cause of this difference.
Relatively recent records of in situ mass balance measurements on this ice cap, together with the several meters of expected amount of snow accumulation during the winter, make Drangajökull an appropriate site for developing the described remote sensing methods.Additionally, Drangajökull's relatively small area makes it suitable for testing Pléiades and WV products (DEMs and orthoimages) because the ice cap is covered entirely or nearly entirely within a single stereo pair, eliminating the need for mosaicking and alignment of multiple DEMs from different dates, which would introduce additional complications and errors.

Satellite stereo images
Two pairs of Pléiades (French Space Agency, CNES) stereo images were acquired over Drangajökull: the first on 14 October 2014 (beginning of the winter) and the second on 22 May 2015 (end of the winter; Table 1 and Fig. 2).An additional dataset of stereo images was acquired from WV2 (Dig-italGlobe Inc via the US National Geospatial Intelligence Agency) on 13 February 2015, covering ∼ 92 % of the ice cap (Table 1 and Fig. 2).
Pléiades and WV2 images have a spatial resolution of 0.7 and 0.5 m at nadir, respectively.The images are encoded in 12 bits (Pléiades) and 11 bits (WV2).The base to height (B / H) ratio from the stereo pairs ranges between 0.4 and 0.5 (Table 1), providing excellent stereo geometry while minimizing occlusions due to steep topography.
The October 2014 Pléiades images were acquired 1 day after the second significant snowfall of the winter (Fig. 2), showing fresh snow covering most of the imaged area.Fine  The images of May 2015 contain areas with clouds at the southern border of Drangajökull, mostly located off-glacier (Fig. 2), and a few thin clouds over the ice cap, though the glacier surface remains visible.The February 2015 orthoimage reveals a similar off-glacier snow extent as the images of May 2015 (Fig. 2).

Lidar
A lidar DEM was produced from an airborne survey in July 2011 (Fig. 1) as part of larger effort to survey all Icelandic glaciers and ice caps from 2008-2012 (Jóhannesson et al., 2013).For Drangajökull, this survey covered an extensive ice-free area outside of the ice cap, up to ∼ 10 km from the ice margin at some locations.The survey was carried out with an Optech ALTM 3100 lidar, with a typical point cloud density of 0.33 pts m −2 .A DEM with 2 m posting was produced from the point cloud (Magnússon et al., 2016a).An uncertainty assessment performed on another lidar dataset from the same sensor acquired in similar conditions revealed an absolute vertical accuracy well within 0.5 m (Jóhannesson et al., 2011).

In situ and meteorological measurements
In situ mass balance measurements have been carried out by the Icelandic Meteorological Office (IMO) and the National Energy Authority on Drangajökull annually since 2005, typically at the end of May (winter mass balance) and again at the end of September (summer mass balance).Snow cores are drilled at six to eight locations at the end of each winter, except for the 2013 campaign (no measurements collected due to bad weather) and the extensive 2014 campaign, where 12 survey sites were measured (Fig. 1).For winter mass balance, the length, volume and weight of each segment of the core drilled were measured, allowing retrieving bulk snow density, snow thickness and the winter mass balance at each location (Fig. 1).Similar procedures for drilling are described in many previous studies (e.g.Guðmundsson, 2000;Thorsteinsson et al., 2002;Ágústsson et al., 2013).The position was measured using a handheld GPS at each core location.We used the in situ data collected at eight of these locations in spring 2015 for data calibration and validation.These measurements were carried on 19 June 2015, which is 1 month later than usual due to an unusually cold spring.All available in situ records of snow density from 2005-2014 were also included in this study.
Additionally, a manually interpolated map of in situ net mass balance for the glaciological year 2013-2014 was obtained using measurements at the 12 mass balance survey sites and a 110 km profile of snow depth from ground penetrating radar (GPR) traversing through all the survey sites (unpublished data, IMO and IES).The locations of survey sites and the GPR profiles were chosen to represent the spatial variation and elevation dependence of the snow cover.The interpolation method is described for a similar dataset by Pálsson et al. (2012).A map of the Drangajökull bedrock topography (Magnússon et al., 2016b) was also used in this study, and daily precipitation and temperature measurements for 2014-2015 from the meteorological station Litla Ávík (LÁ, station #293, 40 km SE of Drangajökull, 15 m a.s.l., Fig. 1) were obtained from IMO (public data, www.vedur.is).

Methods
This section is organized as follows: in Sect.3.1 we describe the processing of remote sensing data to obtain co-registered DEMs, in Sect.3.2 we explain how we derive glacier-wide geodetic winter mass balance from the remote sensing observations and in situ calibration data, and in Sect.3.3 we evaluate the results obtained from remote sensing by comparing them with in situ snow thickness measurements.

Processing of satellite data
Two different schemes (Fig. 3) were used to obtain the DEMs and the difference of DEMs (dDEM), spatially co-registered (e.g., Nuth and Kääb, 2011).Spatial calculations are done in the conformal conic Lambert projection, ISN93 (details at www.lmi.is).Scheme A involves lidar-derived ground control points (GCPs) as a reference, whereas scheme B involves common snow-and ice-free areas in the datasets.From each scheme, statistics of elevation difference in snow-and icefree areas were calculated to verify that the dDEM is unbiased and to quantify its relative accuracy.

Scheme A: processing of Pléiades images using lidar-derived GCPs
The shaded relief lidar DEM was used as a reference for extracting GCPs (Berthier et al., 2014).The GCPs were typically large boulders surrounding the ice cap and on two of the nunataks exposed within the ice cap.These boulders were chosen as GCPs because they are easily recognized in both the lidar hillshade and the stereo images and because they adequately spread horizontally and vertically throughout the study area (e.g., Nuth and Kääb, 2011).Each pair of Pléiades stereo images was processed separately using the ERDAS Imagine (©Intergraph) software as follows: 40 tie points (TPs) were automatically measured on each stereo pair, and an additional 10 GCPs were manually digitized, five of which were common in the October 2014 and May 2015 Pléiades images.The original image's rational polynomial coefficients (RPCs) were thus refined by including the GCPs and TPs in the bundle adjustment.
After RPC refinement, a DEM was produced from each stereo pair by pixel-based stereo-matching with the routine enhanced automatic terrain extraction (eATE).Images were resampled to twice the native pixel size (i.e., to ∼ 1.4 m), which balances the speed of processing and DEM quality.A triangulated irregular network (TIN) was produced from the point cloud and used for sampling a DEM in regular grid spacing of 4 × 4 m.An orthoimage (0.5 × 0.5 m pixel size) was also produced from the image closest to nadir of each pair.
Lidar-derived GCPs from ice-free areas have often been used in photogrammetric studies on glaciers (e.g., James et al., 2006;Berthier et al., 2014;Magnússon et al., 2016a).In the case of Pléiades and WorldView, a few GCPs are sufficient to remove most of the horizontal and vertical biases in the resulting DEMs (Berthier et al., 2014;Shean et al., 2016).

Scheme B: processing of Pléiades images with DEM co-registration
In this approach, the DEMs were produced from the pair of stereo images with the original RPCs.This work was carried out with the open source software Ames Stereo Pipeline (ASP, version 2.5.3)developed by NASA (Shean et al., 2016).The processing chain uses the routine stereo, producing a point cloud from each pair of stereo images, followed by the routine point2dem, which produces a gridded DEM (4 × 4 m grid size) and an orthoimage (0.5 × 0.5 m pixel size) for each pair of stereo images.
Areas with thin semitransparent clouds covering the ice cap in the May 2015 Pléiades images (Fig. 2) produced data gaps in the DEM.These image fragments were processed separately and then mosaicked and superimposed over the initial May 2015 Pléiades DEM and orthoimage.The correlation performed in these areas was based directly on the fullresolution images, instead of a pyramidal correlation from subsample images.This improved the correlation (Shean et al., 2016), resulting in full coverage of these areas (Fig. 2).
The snow-and ice-free areas were delineated from the May 2015 Pléiades orthoimage using a binary mask obtained by setting up a cutoff value of < 0.2 for the top of atmo-sphere absolute reflectance.These images show clear contrast between snow and bare ground (Fig. 2), making image segmentation an efficient approach for the identification of bare ground.
The DEMs were co-registered using the routine pc_align in ASP software, based on the iterative closest point (ICP) algorithm for co-registration of two point clouds (Shean et al., 2016).The ICP was performed in two steps: (1) the snowand ice-free areas of the May 2015 Pléiades DEM were used as a slave DEM, and the entire October 2014 Pléiades DEM was used as a master DEM.A transformation matrix with six parameters (three translations and three rotations) was calculated between the master and slave DEMs.(2) The transformation matrix was applied to the entire May 2015 Pléiades DEM.The applied transformation is quantified by the vector joining the centroids of the May 2015 Pléiades DEM before and after co-registration; this vector has a north component of 8.28 m, a west component of 7.57 m and a vertical component of 12.85 m.A slight planar tilt of 0.002 • was corrected between the DEMs.

February 2015 WV2 DEM
The WV2 data was collected and processed as part of the ongoing US National Science Foundation ArcticDEM project.A gridded DEM with 2 × 2 m grid size was produced with the Surface Extraction with TIN-based Search-space Minimization (SETSM) software (Noh and Howat, 2015), using the RPC sensor model and no GCPs.The 13 February 2015 orthoimage acquired from WV2 was also provided in 2 m pixel size.Since the raw WV2 images were not available in this study, the February 2015 WV2 DEM was co-registered to the October 2014 Pléiades DEM using the ICP algorithm as described in the previous section (scheme B).First, the WV2 DEM, originally in polar stereographic projection, was reprojected and bilinearly resampled to 4 × 4 m.Then, the ICP algorithm was applied to the ice-free areas from the May 2015 Pléiades orthoimage after manually aligning it to the February 2015 WV2 orthoimage and verifying a similar distribution of snow-free areas between the orthoimages of February and May 2015.The vector joining the centroids of the WV2 DEM before and after co-registration has components 10.32 m to the north, 4.63 m to the east and an 8.81 m shift in the vertical.A slight planar tilt of 0.002 • was corrected between the DEMs.

Statistics of elevation differences in snow-and ice-free areas
Statistical indicators of bias and data dispersion were calculated from the dDEM in snow-and ice-free areas using the October 2014 Pléiades DEM as a reference.This included number of cells, median, mean, standard deviation (SD) and normalized median absolute deviation (NMAD, Höhle and Höhle, 2009) over snow-and ice-free terrain.The bare ground areas in the May 2015 images (Fig. 2) were selected for the uncertainty analysis of the dDEM.In the October 2014 Pléiades images, off-glacier snow was on average less than 20 cm thick and therefore negligible in the error analysis (further described in Sect.4.1).
Since the terrain of the ice cap is substantially different (i.e., much smoother) than its ice-free surroundings, statistics were also calculated after filtering snow-and ice-free areas based on (1) a high slope exclusion filter in which snowand ice-free areas with slopes > 20 • were masked out, as performed in previous similar studies (Magnússon et al., 2016a) acknowledging that only 1 % of the ice cap area exhibits slopes larger than 20 • ; and (2) a shadow filter in which shadows were masked out from analytical hillshading (Tarini et al., 2006) using the sun position at the time of acquisition for the respective images.Shadows of the October 2014 Pléiades DEM and February 2015 WV2 DEM revealed much higher levels of noise than sun-exposed areas, and were mostly localized on snow-and ice-free areas, covering < 4 % of the ice cap in the February 2015 WV2 DEM.
DEM uncertainty based on SD or NMAD conservatively assumes totally correlated errors in the dDEM (Rolstad et al., 2009).However, the spatial autocorrelation inherent in the DEM may produce substantially lower uncertainty estimates than calculated by simple statistics (Rolstad et al., 2009;Magnússon et al., 2016a).A sequential Gaussian simulation (SGSim) was performed over the masked snow-and ice-free areas (Magnússon et al., 2016a) in order to calculate a likely bias-corrected mean elevation difference on the ice cap.

Computation of glacier-wide mass balance
Three dDEMs were produced from the different combinations: dDEM t2 t1 , dDEM t3 t2 and dDEM t3 t1 , where t 1 = 14 October 2014, t 2 = 13 February 2015 and t 3 = 22 May 2015 (Fig. 4).The glacier-wide geodetic winter mass balance was calculated as Bw where t f denotes the date of the last DEMs used and h dDEM is the average elevation change over the ice cap observed from the remote sensing data (dDEMs).The term ρ Snowt f is the bulk snow density at the time of the latter DEM, and C t f t 1 represents the spatially averaged densification of the firn layer, h Firn , and the fresh snow, h Snowt 1 , existing on the glacier surface at t 1 .The density and firn densification terms are quantified from field measurements (Sect.3.2.2,3.2.3 and 3.2.4).The mass balance Bw t 3 t 2 is calculated as the difference between Bw t 2 t 1 and Bw t 3 t 1 .Alternatively, the glacier-wide geodetic winter mass balance can be obtained relative to the summer surface, covered by fresh snow at t 1 , assuming that firn or ice does not reap- pear on the glacier surface after t 1 .This approach results in Bw (2) In this case the date of the summer surface is not fixed, and it can vary over the ice cap (Cogley et al., 2011).This surface is, however, typically used as the reference when obtaining the winter balance from in situ mass balance measurements.

Average elevation change
The average elevation change over the ice cap, h dDEM , is extracted from the dDEMs.The extent of the ice cap was digitized from the October 2014 Pléiades orthoimage, following the criteria defined in previous studies (Jóhannesson et al., 2013;Magnússon et al., 2016a) for glacier digitation, which excludes snowfields located at the western and southern sides of the ice cap.We assume that uncertainties in geodetic mass balance caused by digitization of the ice cap outlines are negligible due to the high image resolution.
The data gaps in the dDEMs within the ice cap occur in large shadows north of nunataks in October 2014 and February 2015 and in the south-easternmost part of the ice cap in February 2015 (Fig. 2).These shadows led to < 1 % data gaps for dDEM t3 t1 and ∼ 8 % gaps for dDEM t 2 t 1 and dDEM .

Bulk snow density
The average bulk snow density on Drangajökull at the end of the winter 2014-2015 was ρ Snowt 3 = 554 kg m −3 (SD = 14 kg m −3 ), as deduced from eight snow cores at elevations ranging from 300 to 920 m a.s.l.This density value is used for conversion of volume to water equivalent for the geodetic winter mass balance calculations based on dDEM The midwinter (i.e., 13 February) density of snow is expected to be lower than the bulk snow density measured at the end of the winter.The value ρ Snowt 2 = 500 ± 50 kg m −3 was adopted for the mass balance calculations based on dDEM t 2 t 1 .This lower value of the snow density was observed in a few occasions on Drangajökull during early spring measurements (i.e., 2014 field campaign at the end of March, Fig. 7), and its uncertainty is accordingly large due to the lack of measurements.
The bulk density of snow accumulated for the period 3-14 October, ρ Snowt 1 , is estimated as 400 kg m −3 , which is typical for newly fallen snow on ice caps in Iceland (unpublished data, IES).The bulk density of snow fallen after the May Pléiades images is ρ Snowt 3 −t 4 = 515 kg m −3 , where t 4 = 19 June 2015 (date of the in situ measurements).This is estimated as an average value of snow density in the uppermost segment of each core measured in the field.

Firn densification
Densification of the firn layer leads to a continuous lowering of the bottom of the annual snow pack and an underestimate of snow volume changes estimated from the dDEM (Sold et al., 2013).The total area covered by firn at the end of the 2014 ablation season was 91 km 2 , or about 64 % of the ice cap, based on the extent of snow in a Landsat 8 image acquired on 16 September 2014 (data available from the US Geological Survey, http://earthexplorer.usgs.gov/).Similar spatial distribution of firn areas are inferred from the map of net annual mass balance of the year 2013-2014, showing 58 % of the ice cap with positive mass balance at the end of the summer.
The 2013-2014 net mass balance distribution was used to correct for firn densification, assuming this was a typical year of mass balance for Drangajökull.The net annual surface elevation change due to firn densification vertically integrated over the entire firn column should correspond to the average annual accumulation layer transformed from end-of-the-year snow density to ice (Sold et al., 2013), as where b n+ is the mass balance of 2013-2014 (in units of kg m −2 ) over the accumulation area (positive, by definition), and ρ Firn u and ρ Firn l are the upper and lower values of density of the firn profile, estimated as ρ Firn u = 600 kg m −3 and ρ Firn l = 900 kg m −3 .These values of density in the firn layer are consistent with the literature (Cuffey and Paterson, 2010) and with a measured deep density profile obtained on Hofsjökull ice cap in central Iceland (Thorsteinsson et al., 2002).
For simplicity, the firn densification was distributed linearly over the time span covered (0.603 year for t 3 1 and 0.334 year for t 2 1 ), under the assumption that the firn densification does not vary seasonally.Slight variations can occur in the firn densification process through time, due to accumulation variability and rain -and meltwater percolation (Ligtenberg et al., 2011).The mean values of the firn densification maps, 0.41 and 0.23 m for t 3 1 and t 2 1 respectively, were scaled by the firn area within the ice cap in order to calculate the glacier- The above quantification of the firn densification is based on the mass balance measured extensively during a single year (2013-2014) and assumes equal net accumulation between years as well as a constant densification rate within the glaciological year.An uncertainty of 50 % in the firn correction was used for the error budget of the mass balance (Table 3), due to the assumptions and approximations involved in this method.

Fresh snow densification in the reference DEM
The October 2014 Pléiades DEM, used as a reference for the winter mass balance calculations, contains the first two snowfalls of the winter (Fig. 2), starting on 3 October.This thin snow layer densifies over time from settling, rainfall and compression (e.g., Ligtenberg et al., 2011), causing a lowering of the reference surface and leading to an underestimation of the total winter snow.The snow densification correction was calculated as where W t 1 is the average thickness of the fresh snow (in m w.e.) at t 1 and ρ c is the bulk density of same snow layer at time t f , assuming that the entire fresh snow layer at t 1 is preserved during the period t 1 to t f .ρ c is estimated as 600 kg m −3 for both Bw t 2 t 1 and Bw t 3 t 1 .The first term on the right hand side of Eq. ( 4) corresponds to the h Snowt 1 , which is spatially averaged in Eq. ( 2).The value of W t 1 at a given The Cryosphere, 11, 1501-1517, 2017 www.the-cryosphere.net/11/1501/2017/location was estimated as where P is daily precipitation (in m) and T + is average daily temperature for days when it is above 0 • C, but otherwise T + = 0 • C. α is a snow fall switch, taking the value 1 only if average daily temperature is below 1 • C, otherwise it is 0. β (t * ) takes the value 1 if W t * −1 is positive but is 0 otherwise to avoid accumulation of negative new snow.ddf is a simple degree-day melt factor for snow assumed to be 0.0055 m w.e.• C −1 , as obtained for snow on Langjökull ice cap, central Iceland (Guðmundsson et al., 2009).
The daily precipitation values P were obtained by scaling the daily precipitation values from LÁ for each in situ location by comparison of the net precipitation at LÁ through the entire winter (P LÁ = 0.684 m, Fig. 5) and the measured accumulation at each in situ location, resulting in a scaling factor between ∼ 2 (V1, bw 2014-2015 = 1.54 m w.e.) and ∼ 7 (V6, bw 2014-2015 = 4.93 m w.e.).This assumes that all precipitation that falls on the ice cap through the winter remains in the snowpack, including rain, which is assumed to percolate into the cold snow pack where it refreezes as internal ice layers.The daily temperature values, T , were obtained for each in situ location by projecting temperature records from LÁ, using an elevation lapse rate of −0.006 • C m −1 , as has been measured for Langjökull ice cap (Guðmundsson et al., 2009).
The values of W t 1 and consequently h Snowt 1 were obtained at each in situ site and averaged to obtain the glacier-wide C t f t 1 h Snowt 1 and h Snowt 1 for Eqs.(1) and ( 2) respectively.The in situ locations are fairly evenly distributed over the elevation range of the ice cap and are therefore considered to be representative of the glacier-wide calculations.Based on the observed temporal and spatial variability, we conservatively estimate the uncertainties of h Snowt 1 and C t f t 1 h Snowt 1 to be 50 and 75 %, respectively.

Error propagation
Assuming that the variables in Eq. ( 1 where ρ Snow is the uncertainty in bulk snow density, h dDEM is the uncertainty in average elevation change obtained from dDEM, C {h Firn } is the uncertainty in firn correction and C h Snowt 1 is the uncertainty in snow correction for the reference DEM.Table 2 summarizes the values and uncertainties of each variable affecting the calculation of the geodetic winter mass balance.The uncertainty of Bw t 3 t 2 is calculated as the quadratic sum of uncertainties of Bw t 2 t 1 and Bw t 3 t 1 .The error equation for Eq. ( 2) is analogous to Eq. ( 5), replacing the term C h Snowt 1 by h Snowt 1 .

Comparison of Pléiades-based elevation changes and in situ measurements
For validation of results, the elevation difference at the in situ locations was extracted using bilinear interpolation from dDEM t3 t1 from scheme A, since this scheme is fixed to the same reference frame as the in situ GPS coordinates (lidar frame, Fig. 3).The resulting elevation difference, h dDEM t 3 t 1 was compared with the snow thickness, h Snow in situ , measured at the in situ locations in the 2015 campaign.
Three main factors cause differences in results between the remote sensing and the glaciological method (Sold et al., 2013): (1) the time difference between the DEMs and in situ surveys, (2) firn densification and (3) surface emergence or submergence due to ice dynamics.The corrected satellitebased elevation difference cdDEM t3 t1 for comparison to in situ data is cdDEM where C {h Firn } is the correction due to firn densification (Sect.3.2.3)and h Snow t 1 is the correction due to snow accumulated before t 1 (Sect.3.2.4).h Snowt 3 −t 4 is the correction for snow accumulation and ablation between t 3 (the 22 May Pléiades DEM) and the in situ snow thickness measurements, calculated in the same way as h Snowt 1 , using ρ Snowt 3 −t 4 and allowing for net negative values (i.e., the switch β in Eq. 5 is omitted).dh dyn is the surface emergence and submergence due to ice dynamics (Sect.3.3.1).The magnitude/sign of these corrections differ between the accumulation and ablation areas (Fig. 6).

Ice dynamics
We compare two methods for estimating the effect of ice dynamics on local surface elevation change, dh dyn , during the study period (e.g., Jarosch, 2008;Sold et al., 2013): The emergence and submergence velocities dh dyn icetools were calculated using a full-Stokes ice flow model with the icetools library (Jarosch, 2008) and the finite element package, Fenics.The model calculates a 3-D velocity field resulting from the ice deformation, given the glacier geometry.The bedrock DEM (Magnússon et al., 2016b) and the October 2014 Pléiades DEM were used as inputs.The 2-D horizontal velocities measured with GPS in the 2013-2014 field campaigns were used to calibrate the ice flow rate factor, A. The annual emergence and submergence velocities across the ice cap were computed on a 200 m regular grid and scaled by 0.603, a factor to represent the time span t 1 − t 3 (14 October to 22 May), assuming constant velocities through the glaciological year.
Assuming that the glacier is in a steady state, the longterm average surface net balance (divided by the density of ice) equals in magnitude to the emergence and submergence velocities across the glacier (Sold et al., 2013).Acknowledging that there is significant year-to-year variability in surface net mass balance, the net mass balance measurements from the year 2013-2014, scaled by the water (1000 kg m −3 ) to ice (900 kg m −3 ) conversion factor, were assumed to be representative of local annual emergence and submergence velocities.The obtained values at the in situ locations were then scaled to represent dh dyn bn2013-2014 over the time span t 1 −t 3 .

Uncertainty on elevation difference derived from satellite data
The statistics obtained from the dDEMs in snow-and icefree areas (Table 2) allow for a quantitative comparison of the different methods and datasets used in the study.The statistics show smaller SDs and NMADs outside of the areas of high slopes and shadows due to the dependency of the DEM accuracy on the steepness of the terrain (Toutin, 2002;Müller et al., 2014;Lacroix, 2016;Shean et al., 2016) and the presence of shadows (Shean et al., 2016; Table 2).The vertical bias obtained after DEM co-registration ranges from 0 to 0.1 m based on the median, and the NMAD reveals random errors < 0.5 m in both schemes A and B as well as in the co-registered WV2 DEM.Both schemes yield a similar result for elevation difference, h dDEM , on the ice cap.Details on the distribution of errors in the snow-and ice-free areas, as well as histograms of the distribution, are presented in the Supplement.
The thin layer of snow in the October 2014 Pléiades images (Fig. 2) could slightly skew the statistics.The snow thickness is expected to be less than 20 cm outside the ice cap based on snowfall observations on 13 October at locations V1, V2 and V5 (the closest in situ locations to the icefree areas, Fig. 1), ranging from 0.13 m at V1 (291 m a.s.l.) to 0.27 m at V2 (668 m a.s.l.).The snow line was observed at an elevation of ∼ 50 m a.s.l. in the October 2014 Pléiades images, and the majority (> 60 %) of the cells used for the statistics are at a lower elevation than V1.The results obtained from SGSim provide an uncertainty estimate of 95 % for the dDEM on the ice cap.The SGSim results from both schemes agree well and are within the uncertainty obtained from NMAD in the snow-and ice-free areas, which further supports the robustness of the two methods of DEM processing.All proxies used show almost no bias in the dDEMs (Table 2).The NMAD was kept as a conservative metric for dDEM uncertainty, since the presence of snow in the October 2014 Pléiades images may have affected the results from the SGSim in presumed snow-and ice-free areas, especially in close vicinity of the ice cap, leading to an erroneous bias estimate on the ice cap.

Maps of elevation differences and glacier-wide mass balance
Schemes A and B lead to similar elevation differences and uncertainty based on statistical analyses (Table 2).Since it contains fewer data gaps, scheme B was preferred for pro-ducing elevation difference maps (Fig. 4) and for the study of volume changes and the geodetic mass balance.The firn and fresh snow densification lead to a minor addition (∼ 8 %) to the elevation difference, h dDEM (Table 3).Hence, the maps of dDEMs themselves reveal useful and realistic information about the pattern of snow accumulated in Drangajökull and surroundings (Fig. 4).The western half of the ice cap received more snow than the eastern half, with an average elevation difference h dDEM = 5.91 m between October 2014 and May 2015, in comparison with the eastern half, h dDEM = 5.03 m, during the same period, as suggested in Magnússon et al. (2016a).Significant snow accumulation was also observed in several snowfields outside the ice cap between October 2014 and May 2015.
The glacier-wide geodetic winter mass balance is Bw = 3.33 ± 0.23 m w.e. for the period 14 October 2014-22 May 2015, calculated from Eqs. ( 1) and ( 5).The mass balances obtained for the two periods of the same winter are Bw = 2.08 ± 0.28 m w.e.(14 October 2014 to 13 February 2015) and Bw = 1.26 ± 0.37 m w.e.(13 February to 22 May 2015).The glacier-wide geodetic winter mass balances from the start of the glaciological year, obtained from Eqs. ( 2) and ( 5), are Bw = 3.55 ± 0.27 m w.e. for the period 3 October 2014-22 May 2015, and Bw = 2.27 ± 0.31 m w.e. was calculated from 3 October 2014 to 13 February 2015.We quantify the error of each calculated mass balance and determine the weight of each variable from Eq. ( 5) in the total error budget (Table 3).

Pléiades vs. in situ data
As expected, the in situ measurements of snow thickness yield substantially higher values than the uncorrected difference in elevation measured from dDEM t 3 t 1 (May 2015 Pléiades DEM minus October 2014 Pléiades DEM) in the accumulation area (Fig. 6), with an average difference of 2.56 m for points V3, V6, V7 and J2.Conversely, at Point V1 in the ablation area, the in situ measurements of snow thickness are lower (difference of −0.98 m) than the difference in elevation from dDEM t 3 t 1 .The areas closer to the ELA (points V2, V4 and V5, Fig. 1) show better agreement between glaciological and remote sensing methods before applying corrections (Table 4).
The estimated corrections applied for calculating dDEM t 1 are summarized in Table 4.Each correction has a different impact on the overall comparison, depending on the location of the in situ measurement.The highest corrections were estimated from ice dynamics deduced from the records of mass balance, dh dyn bn , reaching up to 1.69 m of emergence at location V1 in the lower part of the ablation area.Corrections typically span from 0 to 1 m (Table 4).
The estimated correction for the snowfall and ablation in the time difference between the beginning of winter (3 October) and the first satellite acquisition (14 October), h Snowt 1 , assumes the start of winter with the first snowfall, on 3 Ocwww.the-cryosphere.net/11/1501/2017/The Cryosphere, 11, 1501-1517, 2017 Table 3. Glacier-wide geodetic winter mass balance and associated error, calculated from Eq. ( 1).The elevation difference, h dDEM , is observed from remote sensing data, while the bulk snow density (ρ Snow ) and densification of firn (C{h Firn }) and fresh snow (C{h Snowt 1 }) are inferred values from field measurements.For each variable its value and the associated error are shown, and in the row below its conversion into mass balance is shown.Bw ρ Snow shows the contribution of the bulk snow density into the uncertainty in the mass balance.The total uncertainty of Bw is computed as the quadratic sum of the uncertainty (in m w.e.) of the elevation difference, firn and fresh snow densification, and bulk snow density.tober 2014.However, imagery from Landsat and MODIS reveal ice on the low glacier areas in the days before the snowfall on 13 October 2014.At this location it was therefore assumed that the later snowfall marked the beginning of the winter (Table 4).
The mean difference between the in situ measurements and the difference in elevation from dDEM t 3 t 1 is 1.34 m (SD = 1.43,N = 8).The mean difference and its standard deviation are significantly reduced after applying the corrections, obtaining a mean difference of 0.52 m (SD = 0.46) when calculating dDEM t 3 t 1 using dh dyn icetools and a mean of 0.34 m (SD = 0.64) when calculating dDEM t 3 t 1 using dh dyn bn2013-2014 .

Pléiades and WorldView DEMs for measuring snow accumulation
We measure the glacier-wide geodetic mass balance during the winter of 2014-2015, as well as two sub-periods of the same winter, by differencing DEMs obtained from satellite data.In our calculations, we incorporate corrections for snow density and densification of firn and fresh snow, based on in situ measurements.This technique can be applied in small and medium size glaciers (typically ∼ 1000 km 2 can be stereoscopically covered at once based on the capabilities of Pléiades and WorldView), with sufficiently high mass balance amplitude (∼ 0.5-1 m w.e. or higher).The main advantages of using stereoscopic satellite images are repeatability and coverage of remote glaciated areas.The use of external reference data for bundle adjustment prior to stereo correlation, such as lidar-based or GPS-based GCPs, does not improve the relative accuracy of the Pléiades and WorldView DEMs used here (Table 2).
Combining data from Pléiades and WorldView allows for high spatial resolution within a short (3-4 month) interval.The availability of these data and the presented processing strategy allow, to our knowledge, for the first optical satellite-based measurement of winter accumulation on a glacier.Both sensors result in a similar level of accuracy (Table 2) and their combination enables more detailed studies of glacier changes.The ArcticDEM project (data available at http://arcticdemapp. s3-website-us-west-2.amazonaws.com/explorer/)freely offers multi-temporal DEMs of the Arctic region collected since ∼ 2010 with dense temporal repetition (more than 30 DEMs during the last 6 years in certain regions of Greenland, e.g., Willis et al., 2015), therefore providing a high potential for similar studies of geodetic mass balance on seasonal timescales.
The two DEM processing schemes have advantages and disadvantages.Scheme A provides DEMs, orthoimages and dDEMs in an absolute reference system, based on a geodetic network where the lidar DEM is fixed (or similar if GPSbased GCPs are used).This scheme is appropriate when limited unchanged areas are available or if there are identifiable features for extraction of GCPs.This approach, however, requires external spatial information and tedious manual GCP selection.Scheme B uses a highly automated workflow and is independent of spatial information other than the satellite images and camera model information.Co-registration based on scheme B, while ideally requiring well-distributed static control surface, can be applied with an adequate distribution of slope and aspect over limited control surfaces (Shean et al., 2016).The three different processing software (ERDAS Imagine, ASP and SETSM) provided satisfactory results for obtained dDEMs.

Correction of physical glacier phenomena for calculating geodetic winter mass balance
In addition to the remote sensing data, the in situ measurements of the bulk snow density and the densification of the firn layer and fresh snow are needed to retrieve the glacier-wide geodetic winter mass balance (Eqs. 1 and 2).Ice dynamics do not affect the glacier-wide geodetic winter mass balance due to mass conservation (Cuffey and Paterson, 2010).
The Cryosphere, 11, 1501-1517, 2017 www.the-cryosphere.net/11/1501/2017/ .The table lists all corrections applied pointwise to the Pléiades elevation differences dDEM t 3 t 1 to make them comparable to the in situ measurements (see text for details).The table also compares two approaches carried out for correction of surface emergence and submergence velocities: (1) dh dyn icetools using a glacier ice flow model (Jarosch, 2008) and ( 2) dh dyn bn2013-2014 using records of mass balance (Sold et al., 2013).c 1 dDEM The sensitivity of the mass balance calculation was tested with different snow densities measured during the 2005-2014 field campaigns in Drangajökull (Fig. 7).The glacierwide geodetic winter mass balance is reduced by 1 % when the average of all previous density records is used instead of the mean 2015 bulk snow density.The minimum average bulk snow density recorded (511 kg m −3 in 2011) results in 8 % lower mass balance, and the maximum average bulk snow density recorded (583 kg m −3 in 2008) results in a 5 % higher mass balance.We obtained similar discrepancies by using snow density records from other Icelandic ice caps.Bulk snow density measured on Mýrdalsjökull ice cap in 2010 (Ágústsson et al., 2013) and on Langjökull ice cap in 2015 produced a 3 and 10 % overestimation and underestimation of mass balance, respectively.
Bulk snow density can vary substantially between different glaciers or between different years in the same area.Individual years, however, show relatively low scatter of bulk snow density distribution over the different in situ locations on Drangajökull (Fig. 7).The low scatter indicates that bulk snow density measurements taken at one or many points on a date close to that of the satellite acquisitions, if adequately selected for the whole ice cap, should give reasonable results for glacier-wide geodetic winter mass balance calculations.
The firn densification model assumes a temporally constant annual mass balance in the accumulation area, which is a significant source of uncertainty due to high inter-annual climate variability.Other methods can be used for a more accurate correction for firn densification, such as deep core drilling (Thorsteinsson et al., 2002), or robust firn layer observations and modeling (e.g., Sold et al., 2015).For large areas, such as catchments of the Greenland Ice Sheet, a firn densification model such as IMAU-FDM (Ligtenberg et al., 2011), forced by a surface mass balance model such as the RACMO2.3(Noël et al., 2015) can also be applied.However, the resolution (typically 11 km) of these models may be too coarse to resolve a relatively small Icelandic ice cap such as Drangajökull.
The densification caused by fresh snow potentially present at the time of acquisitions of the reference (initial) DEM needs to be studied differently for each case and will depend www.the-cryosphere.net/11/1501/2017/The Cryosphere, 11, 1501-1517, 2017 on the amount of snow falling between the beginning of the glaciological winter and the satellite acquisition.If satellite images are acquired prior to the start of the winter, this effect disappears, and a correction due to surface melt should be assessed (e.g., by using a degree-day model as in Eq. 5).Densification of fresh snow corrected by Eq. ( 1) leads to smaller uncertainty than shifting the mass balance to the beginning of the season using Eq. ( 2), and the uncertainty associated with Eq. (2) will increase with the length of the time period from the start of the winter to t 1 .Firn and fresh snow densification have little effect on the geodetic winter mass balance, increasing it by 8 % (Table 3), indicating that even if these variables remain unknown (i.e., in remote areas), adequate calculations of geodetic mass balance can be performed with moderately increased uncertainties, ranging between 5 and 10 % for glaciers with mass balance amplitude similar to Drangajökull.The error in geodetic mass balance is primary controlled by our knowledge of physical glacier phenomena (bulk snow density and densification of firn and fresh snow) and, to a lesser degree, by the accuracy of the derived maps of elevation differences from the satellite data (Table 3).

Validation of results: remote sensing vs. in situ
The glacier-wide geodetic mass balances suggest that ∼ 60 % of the winter accumulation occurred during the first 4 months of the winter (14 October 2014-13 February 2015, Table 3).Precipitation records at a weather station ∼ 40 km from the ice cap indicate the same ratio of accumulation for the two time periods: 342 mm (62 % of total) between 14 October 2014 and 13 February 2015 and 218 mm (38 % of total) between 13 February and 22 May 2015 (Fig. 5).The consistency of the ratio of accumulation in the two sub-periods observed at the weather station and calculated from the satellite images is encouraging and also supports the applicability of the corrections applied due to differences in time between in situ and geodetic mass balance observations.The temporal offset between the glaciological and the geodetic measurements results in some ambiguity in the definition of the beginning and end of the mass balance season.Glaciological measurements generally use the previous summer layer as reference, which ensures a well-defined starting point of the mass balance year, despite the fact that the date chosen for the spring campaign (i.e., the winter balance end date) is not objectively defined.For example, two snow events occurred in late May and early June, which can either be considered part of the winter or summer balance seasons.The timing of remote sensing surveys are further dependent on sensor tasking and favorable weather (cloud-free) conditions, and, as a consequence, a temporal offset between glaciological and geodetic observations is likely to occur.
The points V1-V4 are located at Leirufjarðarjökull (Fig. 1), a surge-type glacier (Björnsson et al., 2003;Brynjólfsson et al., 2016).The dynamics of this glacier outlet are, by nature, not in balance with the rate of accumulation or ablation, and thus the calculation of emergence and submergence velocities from the net annual mass balance is inappropriate at these locations.On the other hand, an underestimation of submergence velocities is observed over the southern areas using the full-Stokes ice flow model, possibly explained by the lack of basal sliding in the ice flow model.Only minor elevation changes were detected in this part of the glacier in the past decades (Magnússon et al., 2016a), and it is not known to surge; hence, the net annual mass balance approach may be more suitable in this area.

Conclusions
This study shows the capabilities of sub-meter satellite stereo images for measuring winter mass balance.The DEMs created from Pléiades and WV2 satellite stereo images reveal relative accuracy of 0.2-0.3m (for slopes < 20 • ), which allows measuring the evolution of snow accumulation in two periods of the winter on Drangajökull ice cap.Two methodologies used for the processing of DEMs yield similar accuracy and elevation changes with and without using GCPs, showing that the processing of modern sub-meter satellite stereo images for measuring glacier elevation change can be performed without external reference data, such as lidar or GPS data, as long as areas of stable (snow-and ice-free) terrain are present in the imagery to serve as relative control.
The glacier-wide geodetic winter mass balance was 3.33± 0.23 m w.e. for 14 October 2014-22 May 2015, with ∼ 60 % of the accumulation occurring between 14 October 2014 and 13 February 2015.Besides the remote sensing observations, the glacier-wide geodetic winter mass balance calculation requires knowledge of the bulk snow density for volume to water equivalent conversion and a correction for firn and fresh snow densification, which are estimated in this study from in situ measurements.The uncertainty in the bulk snow density is the largest contributor to the uncertainty in glacier-wide geodetic winter mass balance and is significantly larger than the uncertainty in the average elevation change and the firn and fresh snow densification.
Densification of firn and fresh snow produces a systematic but minor (8 %) increase to the mass balance obtained from the geodetic method.This contribution may vary for individual cases depending on the climatic conditions and the timing of snowfall events relative to reference (i.e., start of winter) image acquisition.Uncertainties in geodetic winter mass balance can be minimized with records of bulk snow density and previous years' mass balance.Extrapolation of snow density from other glaciers with different characteristics can, however, lead to slightly larger errors (up to 10 %).
The satellite-derived map of elevation change and eight in situ measurements of snow thickness are in agreement after correcting for three phenomena of sub-meter to meterlevel elevation change: (1) the difference in time between The Cryosphere, 11, 1501Cryosphere, 11, -1517Cryosphere, 11, , 2017 www.the-cryosphere.net/11/1501/2017/ in situ campaigns and satellite acquisitions, (2) the effect of firn densification in the accumulation area and (3) the vertical component of the ice flow motion.While glacier winter mass balance measurements have been sparse due to the difficulty in obtaining field measurements and the low contrast of snow-covered terrain preventing photogrammetric surveying, we demonstrate that sub-meter satellite imagery may offer a powerful new tool for glacier mass balance monitoring on sub-annual timescale.The potential for this approach is enhanced by the rapid increase and availability of optical satellites collecting stereo images in glaciated regions with dense temporal resolution.Due to the relative accuracy of the DEMs and uncertainties in snow density and firn and fresh snow densification, repeated DEMs are capable of obtaining useful estimates of the glacier-wide seasonal mass balance in areas where expected mean thickness of winter snow exceeds 1 m.The accuracy is improved significantly when satellite data and in situ information are combined.

Figure 1 .
Figure 1.Area of study and data collected.(Left) Mosaic of Iceland from Landsat 8 images, mosaicked by the National Land Survey of Iceland.The blue rectangle locates the Drangajökull ice cap, and a blue dot indicates the location of the meteorological station "Litla Ávík" (LÁ).L, M, V and H represent the locations of Langjökull, Mýrdalsjökull, Vatnajökull and Hofsjökull ice caps, respectively.(Right) A shaded relief representation of a lidar DEM covering Drangajökull and vicinity in the summer 2011 (Jóhannesson et al., 2013).Margins of the ice cap are shown as a black polygon, and the equilibrium line altitude (ELA) obtained from the mass balance measurements over 2013-2014 is shown with a green dashed line.Blue dots indicate location of the in situ measurements.Locations labeled V1-7 have been measured since 2005, whereas locations labeled J1-5 were only measured in 2014 except J2, which was also measured in 2015.Black rectangles show the footprints of the Pléiades images, and a green rectangle shows the footprint of the WV2 DEM.

Figure 3 .
Figure 3. Flowchart of the different schemes studied for obtaining unbiased DEMs.Rectangles indicate processing steps and parallelograms indicate products.Orange squares indicate processing with ERDAS software, and green squares indicate processing with ASP software.

Figure 4 .
Figure 4. Elevation difference based on Pléiades and WV2 data.(a) Elevation difference from October 2014 Pléiades DEM to February 2015 WV2 DEM.(b) Elevation difference from February 2015 WV2 DEM to May 2015 Pléiades DEM.(c) Elevation difference from October 2014 Pléiades DEM to May 2015 Pléiades DEM.A black polygon indicates the glacier margin in October 2014.The yellow dashed line shows the boundary between the eastern and western halves of the ice cap.Contours on the ice cap were smoothed with a Gaussian filter of 9 × 9 window size.(d) Longitudinal profile A-A with surface elevation (black line, in meters above ellipsoid, m a.e.) and snow thickness (blue) over the glacier and ice-free areas.The red dashed lines indicate the location of the glacier margins.

t 3 t 1 .
The estimated uncertainty in bulk snow density is ±27 kg m −3 , obtained from the SD from all available bulk snow density measurements in Drangajökull since the first field campaign in 2005.This error includes the uncertainty in density caused by (1) errors in measurements and (2) likely snow densification between the May 2015 Pléiades images and the June 2015 field campaign.

Figure 5 .
Figure 5. Cumulative precipitation (clear blue) and temperature (red line) for the winter 2014-2015 (1 October 2014 to 19 June 2015) from the station Litla Ávík.Blue dashed lines show the time of acquisition of satellite stereo images.
) are not correlated to one another, the error in the mass balance calculation is ob-

Figure 6 .
Figure 6.Sketch of the different factors, marked in red and indicated with red arrows, affecting the comparison between the glaciological (3 October 2014-19 June 2015) and geodetic (14 October 2014-22 May 2015) methods.Light blue represents snow fallen in winter, and dark blue represents preexisting ice and firn.

Figure 7 .
Figure 7.The density values obtained at each in situ location for field campaigns 2005-2015.Each circle represents the average density of the shallow core at each in situ location.Blue filled circles show the average density measurements.Black "+" shows the averaged density measured on Langjökull, and black "×" shows the averaged density measured on Mýrdalsjökull ice cap in year 2010(Ágústsson et al., 2013).The 2013 campaign was not carried out due to bad weather conditions.

Table 1 .
Dates, type of data (split between remote sensing and in situ data), sampling and specifications of datasets used in this study.

Table 2 .
Magnússon et al., 2016a)the dDEMs in snow-and ice-free areas, and mean elevation difference on the ice cap, h dDEM .N represents number of data points.The three bottom rows indicate the statistics after masking slopes > 20 • and shadows.Bias-corrected SGSim represents the mean elevation bias from 1000 simulations and the standard deviation of the simulations (details inMagnússon et al., 2016a).

Table 4 .
Comparison of values of snow thickness, h snow insitu , measured in the field, and elevation difference obtained from Pléiades DEMs, h 1 , using the two different approaches, and Res 1 and Res 2 are the residuals between the glaciological and geodetic methods after applying the corrections.C{h Firn } h Snow h Snow dh dyn c 1 dDEM t