Fusion of multi-source satellite data and DEMs to create a new glacier inventory for Novaya Zemlya

: Monitoring glacier changes in remote Arctic regions are strongly facilitated by satellite data. This is especially true for the Russian Arctic where recently increased optical and SAR satellite imagery (Landsat 8 OLI, Sentinel 1/2), and digital elevation models (TanDEM-X, ArcticDEM) are becoming available. These datasets offer new possibilities to create high-quality glacier inventories. Here, we present a new glacier inventory derived from a fusion of multi-source satellite data for Novaya Zemlya in the Russian Arctic. We mainly used Landsat 8 OLI data to automatically map glaciers with the band ratio method. Missing debris-covered glacier parts and misclassified lakes were manually corrected. Whereas perennial snow fields were a major obstacle in glacier identification, seasonal snow was identified and removed using Landsat 5 TM scenes from the year 1998. Drainage basins were derived semi-automatically using the ArcticDEM (gap-filled by the ASTER GDEM V2) and manually corrected using fringes from ALOS PALSAR. The new glacier inventory gives a glacierized area of 22,379 ± 246.16 km2 with 1474 glacier entities >0.05 km2. The region is dominated by large glaciers, as 909 glaciers <0.5 km2 (62% by number) cover only 156 ± 1.7 km2 or 0.7% of the area, whereas 49 glaciers >100 km2 (3.3% by number) cover 18,724 ± 205.9 km2 or 84%. In total, 41 glaciers are marine terminating covering an area of 16,063.7 ± 118.8 km2. The mean elevation is 596 m for all glaciers in the study region (528 m in the northern part, 641 in the southern part). South-east (north-west) facing glaciers cover >35% (20%) of the area. For the smaller glaciers in the southern part we calculated an area loss of 5% (52.5 ± 4.5 Abstract: Monitoring glacier changes in remote Arctic regions are strongly facilitated by satellite data. This is especially true for the Russian Arctic where recently increased optical and SAR satellite imagery (Landsat 8 OLI, Sentinel 1/2), and digital elevation models (TanDEM-X, ArcticDEM) are becoming available. These datasets offer new possibilities to create high-quality glacier inventories. Here, we present a new glacier inventory derived from a fusion of multi-source satellite data for Novaya Zemlya in the Russian Arctic. We mainly used Landsat 8 OLI data to automatically map glaciers with the band ratio method. Missing debris-covered glacier parts and misclassiﬁed lakes were manually corrected. Whereas perennial snow ﬁelds were a major obstacle in glacier identiﬁcation, seasonal snow was identiﬁed and removed using Landsat 5 TM scenes from the year 1998. Drainage basins were derived semi-automatically using the ArcticDEM (gap-ﬁlled by the ASTER GDEM V2) and manually corrected using fringes from ALOS PALSAR. The new glacier inventory gives a glacierized area of 22,379 ± 246.16 km 2 with 1474 glacier entities >0.05 km 2 . The region is dominated by large glaciers, as 909 glaciers <0.5 km 2 (62% by number) cover only 156 ± 1.7 km 2 or 0.7% of the area, whereas 49 glaciers >100 km 2 (3.3% by number) cover 18,724 ± 205.9 km 2 or 84%. In total, 41 glaciers are marine terminating covering an area of 16,063.7 ± 118.8 km 2 . The mean elevation is 596 m for all glaciers in the study region (528 m in the northern part, 641 in the southern part). South-east (north-west) facing glaciers cover >35% (20%) of the area. For the smaller glaciers in the southern part we calculated an area loss of ~5% (52.5 ± 4.5 km 2 ) from 2001 to 2016.


Introduction
Glacier outlines are a key dataset for analysing glacier fluctuations (changes in length [1], area [2,3] and volume/mass [4]), determining their characteristics (e.g., [1,5]), run-off contribution [6] or future behaviour [7]. In addition, they allow spatially constrained calculations of glacier specific products such as flow velocities and stable terrain for accuracy assessment of glacier volume changes derived from two digital elevation models (DEMs) [8]. Due to the strong glacier shrinkage observed globally [9,10], it is necessary to update their extent at about decadal time-scales [11]. This is also true for the high Arctic, which has warmed almost twice as fast as the global average over the last few decades, likely due to feedback mechanisms from diminishing sea ice cover [12]. As glaciers in the Russian Arctic contain about 10% of the total glacier mass [9], they could potentially make a large contribution to sea level rise. However, the entire region lacks field-based mass balance observations using the direct glaciological method or length change measurements [10] and our current knowledge about their past changes is mostly derived from remote sensing data. From the islands in the Russian Arctic, reference [13] reported the strongest mass loss for Novaya Zemlya (NVZ) using ICESat laser altimetry and GRACE gravimetry and [14] also observed mostly retreating outlet glaciers over the 1992 to 2010 period for its northern ice cap, in particular for marine terminating glaciers. More recently, reference [15] reported high flow velocities and substantial thinning near the glacier fronts for NVZ.
The comparison of flow velocity fields derived from radar interferometry also revealed that the drainage divides provided in the Randolph Glacier Inventory 5.0 (RGI) for the separation of glacier entities had partly poor quality as they crossed the flow lines [16]. Moreover, the only available dataset of glacier outlines in the Global Land Ice Measurements from Space (GLIMS) database and the RGI is partly generalized (roughly digitized), has no clear size threshold, and most glaciers have only a date range (of about 10 years according to [13]) rather than a clear acquisition date. For modelling purposes, it is desirable to have both a clear date for each glacier polygon and only a few years difference for a larger region. This is in especially true for NVZ as several of its glaciers have been identified as surge type [17].
This study presents a new glacier inventory (outlines and topographic characteristics) for Novaya Zemlya. Due to varying image quality and availability, Landsat 8 OLI scenes of 2013 and 2015 were used for the northern region (dominated by one large ice cap) and a scene from 2016 for the smaller mountain and valley glaciers in the south. The main purpose of this study is to provide (a) updated and (b) temporally clearly defined glacier extents covering (c) a short time span along with (d) improved drainage divides and (e) complete topographic attributes for each glacier. The latter two items are derived from the recently published very high-resolution ArcticDEM [18] and interferometric fringes from ALOS PALSAR SAR images. Additionally, we used the ASTER GDEM V2 to fill data voids in the ArcticDEM and a Landsat 5 TM scene from 1998 to remove seasonal snow. Whereas the large ice cap and the fluctuations of its outlet glaciers have already been studied in detail, little is known for the local glaciers and ice caps (GIC) in the south. We have thus also derived glacier area changes between 2001 and 2016 for this region.

Study Region
The study region covers the glaciers and ice caps on the two islands of NVZ in the Russian Arctic stretching from about 73 to 77 • N ( Figure 1). Small valley and mountain glaciers dominate the southern island and are present up to 74.1 • N on the northern island. Further north, an isolated mountain range that is not connected to the large northern ice cap can be found. Another small mountain chain is attached to the large ice cap, but its glaciers are glaciologically independent as they are only connected in the ablation area. The large main ice cap on the northern island is the most prominent feature of NVZ. It has a length of about~400 km along the main north-south ridge and a mean east-west extension of~80 km. Its maximum elevation is 1600 m a.s.l., whereas the southern part reaches heights of about 1340 m a.s.l. [19]. Different types of glaciers are found on NVZ but most of them are ice caps and their outlet glaciers or simple cirque and valley glaciers. A special type is the niche glacier that can be found in the abundant topographic depressions where drifting snow can accumulate. As the possible ice in these locations is seldom visible and might not flow, these features can be named perennial snowfields that are not counted in our glacier inventory. Whereas most of the glaciers that are separated from the main ice cap are land terminating, the large outlet glaciers of the main ice cap are mostly marine terminating, so that calving plays an important role in their mass balance [13][14][15]. In addition, several glaciers have been identified as surge type [17]. Debris cover on glaciers is rare. It is assumed that glaciers in NVZ are predominantly cold based [20].
The climate in NVZ is influenced by the North Atlantic-derived water of the Barents Sea where the island provides an orographic barrier for eastwards moving Atlantic cyclonic systems [21]. The warmest and wettest conditions (~400 mm precipitation) are found in the southwest (~−5 • C) and become colder (~−10 • C) and dryer (~250 mm precipitation) to the northeast [21]. Precipitation exhibits large seasonal and interannual variation with peak precipitation rates in August and driest conditions in April.

Satellite Data
For the glacier mapping we selected seven of the most suitable (with minimum seasonal snow and largely cloud free) Landsat 8 Operational Land Imager (OLI) standard terrain corrected (Level 1T) scenes available from the glovis.usgs.gov archive for 2013/2015 (each 3 scenes) and 2016 (1 scene, Table 1). In the northern part of the island, clouds prohibited glacier mapping either on the eastern Figure 1. The islands of Novaya Zemlya, with outlines from the new NVZ glacier inventory and the RGI 5.0. Landsat footprints are shown in black and stable terrain for DEM co-registration is in pink.

Satellite Data
For the glacier mapping we selected seven of the most suitable (with minimum seasonal snow and largely cloud free) Landsat 8 Operational Land Imager (OLI) standard terrain corrected (Level 1T) scenes available from the glovis.usgs.gov archive for 2013/2015 (each 3 scenes) and 2016 (1 scene, Table 1). In the northern part of the island, clouds prohibited glacier mapping either on the eastern or western side using only one Landsat 8 scene, so here a mosaic of two Landsat 8 scenes was used (B, C Remote Sens. 2017, 9, 1122 4 of 19 in Table 1). Moreover, seasonal snow was a significant problem in the central-western part of the island. In this case, we used two Landsat 5 scenes from 1998 (E, H) with near perfect mapping conditions to map glaciers without seasonal snow and the Landsat 8 scenes from 2015 (D, G) for mapping the front position of the glaciers. In this study, we additionally processed L-band Synthetic Aperture Radar (SAR) data from the first Advanced Land Observation System (ALOS) Phased-array L-band Synthetic-Aperture Radar (PALSAR) developed by the Japan Aerospace Exploration Agency (JAXA) and the Japanese Resources Observation System Organization (JAROS) ( Table 2). The radar operated at a frequency of 1.27 GHz (24 cm wavelength) from 2006 to 2011. The data analyzed herein were acquired from 2008 to 2010 along ascending (evening) orbits of the satellite, at horizontal receive and transmit polarization (HH), in Fine Beam Single (FBS) mode, with a 35 • look angle off nadir, a 70 km swath width, and 4.7 m × 3.4 m resolution in slant-range (across-track) and azimuth (along-track) direction, respectively. The data were obtained as raw individual frames from JAXA via the European Space Agency (ESA). Frames were concatenated and processed into single look complex tracks using the Gamma Remote Sensing processor (http://www.gamma-rs.ch).

Digital Elevation Models
For the new glacier inventory we decided to use the new ArcticDEM from the National Geospatial-Intelligence Agency of the National Science Foundation released in September 2016 (Release 3: v2.0) ( Table 3). The ArcticDEM is derived from Worldview and GeoEye-1 optical stereo imagery in an automated processing chain [18] and is provided in a stripe (with 2 m resolution) or mosaic (5 m) format. The projection is polar stereographic north as this is the most suitable metric projection that can be applied over the entire Arctic without the need for zones. We used the 5 m mosaicked version acquired between 2013-2015, as it has been feathered to reduce gaps and edge-matching artifacts. The datasets (IDs: 43_58-43_59, 44_55-44_59, 45_54-45_57, 46_54-46_55) were obtained from the data repository of the Polar Geospatial Center of the University of Minnesota.
Remote Sens. 2017, 9, 1122 5 of 19 The data gaps in the ArcticDEM (mostly small gaps in the margin of the island and sometimes large gaps probably due to interpolation difficulties) were filled with ASTER GDEM V2 data (Table 3) available from the NASA Reverb homepage released in October 2011 by NASA and Japan's Ministry of Economy, Trade and Industry (METI) [22,23]. The ASTER GDEM V2 was compiled photogrammetrically by automated digital matching of stereo pairs from all suitable ASTER images acquired between 2000 and 2012 at a horizontal resolution of 30 m and georeferenced to the WGS84 ellipsoid with geographic coordinates (latitude, longitude) [24,25]. The root mean square error (RMSE) for the ASTER GDEM V2 is given as 8.7 m, and the absolute vertical accuracy, expressed as a linear error at the 95% confidence level, is 17.0 m [26].
The intermediate DEM (IDEM) from TanDEM-X is based on NVZ on acquisitions from 12 December 2010 to 26 March 2012 (Table 3). It is provided in geographic coordinates with a horizontal resolution of 3 , corresponding to approximately 90 m. The indicated absolute horizontal and vertical accuracies are <10 m [27]. The TanDEM-X IDEM has been used during the SAR processing for topographic-phase removal and geocoding. We also used outlines from the RGI version 5 as a guide and for change assessment in the southern part of NVZ [28,29]. The glacier outlines were derived by manual digitalization from SPOT5 and Landsat imagery acquired between 2000 and 2010 [13]. The RGI 5.0 gives a glacierized area of 22,128.2 km 2 (southern part: 1604.7 km 2 ) with 480 (238) glacier entities [13]. A distinctive size threshold is not particularly indicated (we estimate the minimum size threshold is 0.25 km 2 ) and thus smaller glaciers are often not mapped. There is further some generalization due to manual digitation of all outlines. Drainage divides were derived from an unknown DEM and glacier entities lack an exact acquisition date, i.e., only the date range from 2000 to 2010 is given.

Glacier Mapping, Change Assessment, Quality Assessment and RGI Acquisition Date Determination
Glaciers were mapped with the well-established semi-automated band ratio method (e.g., [30]) using raw digital numbers of Landsat 8 OLI bands 6 (red) and 8 (shortwave infrared, SWIR) and bands 3 (red) and 5 (SWIR) for Landsat 5. The optimal threshold for the ratio image (red/SWIR) was chosen interactively for each scene with pixels being classified as ice for a ratio >2 for OLI and >1.6 for TM (scene dependent). The resulting binary map was smoothed with a median filter (3 × 3 kernel) to reduce noise and converted into a vector shapefile format with polygon topology. Manual editing was applied to remove sea ice, correct outlines of some glaciers in shadow and with debris cover (only few), and to remove lakes (mainly in the south). More challenging was the manual editing of calving termini, which were sometimes difficult to separate from sea ice. Clouds were on either side of the island for the footprint 176/006, so that we used two scenes for this region to map glaciers. Most of the seasonal and perennial snow was removed by applying a size threshold of 0.05 km 2 . For the region that was impacted by large amounts of seasonal snow (scene E/D and H/G in Figure 1), the Landsat images from 1998 were used for the initial glacier mapping and the scenes from 2015 were used to update all terminus positions ( Figure 2).

Figure 2.
Raw glacier mapping from 2015 (yellow outlines) overlaid with a Landsat imagery from 1998 (band composition: SWIR, NIR, RED) used for glacier mapping due to better snow conditions. The glacier with magenta outlines was used for DEM comparison over glacierized terrain.
Accuracy assessment of glacier outlines is difficult as appropriate reference data (i.e., with a higher spatial resolution and the same snow conditions) are often missing. Moreover, the uncertainty for glacier outlines is higher when manually edited (shadow, debris cover, calving termini) and for smaller glaciers than for larger ones. For these reasons, we calculated three uncertainty values: for the southern part, the northern part and the entire region. For clean ice the variability in the location of outlines was found to be about one pixel (depending on the operator's threshold used) [31]. This value was thus used to determine the precision of the outlines with the buffer method, i.e., expanding/reducing their horizontal extent by ±15 m outward/inwards. Before the calculation, all glacier entities of the northern part were merged into one polygon (to avoid double counting of buffers on the ice cap), whereas in the southern part glacier entities were not merged. The buffering increased the area for the northern part to 152.7 km 2 (~0.73%) and decreased it by 154.57 km 2 (~0.75%) and for the southern part to 102.38 km 2 (~5.99%) and 89.96 km 2 (~5.93%), respectively. Subsequently, we calculated the mean of the area differences and its percentage from the total for both regions (north: 153.63 km 2 or 0.74%, south: 96.17 km 2 or 5.96%) and used the percentages as a measure of uncertainty. The uncertainty for all glacier outlines in NVZ is given by the sum of both area uncertainties of the southern and the northern part (153.6 km 2 + 96.1 km 2 = 249.8 km 2 ) or 1.1% of the total area.
As the acquisition date of the RGI 5.0 outlines of NVZ is not known, we compared outlines visually to Landsat quick looks to determine their possible acquisition date. A scene from 8.8.2001 was identified as the possible source for the outlines in the southern part of NVZ. Accuracy assessment of glacier outlines is difficult as appropriate reference data (i.e., with a higher spatial resolution and the same snow conditions) are often missing. Moreover, the uncertainty for glacier outlines is higher when manually edited (shadow, debris cover, calving termini) and for smaller glaciers than for larger ones. For these reasons, we calculated three uncertainty values: for the southern part, the northern part and the entire region. For clean ice the variability in the location of outlines was found to be about one pixel (depending on the operator's threshold used) [31]. This value was thus used to determine the precision of the outlines with the buffer method, i.e., expanding/reducing their horizontal extent by ±15 m outward/inwards. Before the calculation, all glacier entities of the northern part were merged into one polygon (to avoid double counting of buffers on the ice cap), whereas in the southern part glacier entities were not merged. The buffering increased the area for the northern part to 152.7 km 2 (~0.73%) and decreased it by 154.57 km 2 (~0.75%) and for the southern part to 102.38 km 2 (~5.99%) and 89.96 km 2 (~5.93%), respectively. Subsequently, we calculated the mean of the area differences and its percentage from the total for both regions (north: 153.63 km 2 or 0.74%, south: 96.17 km 2 or 5.96%) and used the percentages as a measure of uncertainty. The uncertainty for all glacier outlines in NVZ is given by the sum of both area uncertainties of the southern and the northern part (153.6 km 2 + 96.1 km 2 = 249.8 km 2 ) or 1.1% of the total area.
As the acquisition date of the RGI 5.0 outlines of NVZ is not known, we compared outlines visually to Landsat quick looks to determine their possible acquisition date. A scene from 8.8.2001 was identified as the possible source for the outlines in the southern part of NVZ.
In this study, we additionally calculated glacier-specific area changes for the southern part of the study region ( Figure 1). Deriving this change is more difficult as drainage basins need to be the same and each glacier still should be represented as one polygon in both time periods. We have therefore merged the RGI 5.0 outlines and intersected them with the drainage basins obtained in this study. We manually selected 222 glaciers in the southern part (represented by one polygon) and derived area changes for each of them.

SAR Processing (Fringes)
We computed a total of eight interferometric phase images from the ALOS-1 PALSAR-1 with baselines shorter than 1000 m at 100 m postings ( Table 2) [32]. Differential SAR interferograms were processed in the two-pass approach using the TanDEM-X IDEM for topographic-phase removal and geocoding [27]. Standard processing steps included co-registration of the SAR acquisitions, multi-looked interferogram generation using a 4 × 8 window, estimation and removal of the topographic phase using the IDEM, and terrain corrected geocoding. The acquisition date of the IDEM matches that of the ALOS-1 PALSAR-1 dataset close enough to not expect any major topographic signal left on the differential interferograms. In addition, because the perpendicular baseline of the ALOS-1 PALSAR-1 image pairs is below 1000 m (Table 2), a topographic height change associated with a 2P cycle is about 60 m whereas the TanDEM-X vertical accuracy is below 10 m. Hence, phase signals can be interpreted as ice surface displacement in the satellite line-of-sight direction with possible atmospheric disturbances. Whereas coherence is low for the lower sections of the outlet glaciers, it was sufficiently high over the slower moving interior parts of the ice cap to obtain clearly visible fringes.

DEM Fusion
The 42 tiles of the ArcticDEM were first mosaicked and then resampled bilinearly (for smoothing of artefacts) to a resolution of 15 m. The 34 granules of the ASTER GDEM V2 were also first mosaicked, then resampled bilinearly to 15 m and finally projected to a polar stereographic projection. Before merging both DEMs, they were co-registered according to [33]. In this case, we defined the ArcticDEM as the reference (master), the ASTER GDEM V2 as the slave DEM, and digitized four areas around the island where we assumed terrain is stable ( Figure 1). The co-registration revealed shifts of the ASTER GDEM V2 in x (10 m), y (−16 m) and z (11 m) directions. After co-registration, the mean elevation difference of the stable terrain decreased to −1 m (before 9.3 m) and the standard deviation to 9 m (before: 9.09). The resulting dataset was subsequently layer-stacked to have the same x and y extensions, and no-data values in the ArcticDEM were replaced with values from the ASTER GDEM V2 (Figure 3). This procedure allowed us to close most of the data gaps in the ArcticDEM (in the following referred as the fused DEM).

Drainage Divides
To split contiguous ice masses into individual glacier entities, digital drainage divides are required.
Here we use the semi-automated method described in [34], which is a combination of the methods by [35,36] and uses standard GIS functions for watershed analysis. As we use a different DEM (fused DEM) for the divides as it was used for orthorectification of the Landsat scenes by USGS, a large number of sliver polygons resulted after digital intersection along mountain crests.
The resulting drainage basins were different in comparison to the divides in RGI 5.0 (cf. Figure 1). However, a closer visual comparison with the flow-velocities derived from SAR data revealed further inconsistencies (Figure 4a) [37]. If the outlines of the drainage basins derived from the fused DEM Remote Sens. 2017, 9,1122 8 of 19 did not follow the fringe patterns, they were manually edited to the extent possible. For consistency, the divides in the previous inventory guided the final selection of basins. In the final product, divides no longer cross surface flow but contain the correct glacier entities (Figure 4b).

Topographic Parameters
After glacier complexes (the outcome of the mapping, basically one large polygon) were digitally intersected with the drainage divides (cf. Figure 4b), zonal statistics were calculated based on the glacier outlines and the fused DEM or products derived from it (e.g., slope and aspect grids) as proposed by [38]. The calculation of the topographic parameters (minimum, mean, median and maximum elevation, mean slope and aspect) for each glacier were performed according to [39]. Furthermore, glacier hypsometry (area-elevation distribution), and size/number distribution statistics are derived both for the entire region and the southern part.

Topographic Parameters
After glacier complexes (the outcome of the mapping, basically one large polygon) were digitally intersected with the drainage divides (cf. Figure 4b), zonal statistics were calculated based on the glacier outlines and the fused DEM or products derived from it (e.g., slope and aspect grids) as proposed by [38]. The calculation of the topographic parameters (minimum, mean, median and maximum elevation, mean slope and aspect) for each glacier were performed according to [39]. Furthermore, glacier hypsometry (area-elevation distribution), and size/number distribution statistics are derived both for the entire region and the southern part.

Glacier Characteristics
The new glacier inventory for NVZ has an area of 22,379.0 ± 246.2 km 2 with 1474 glacier entities when considering only glaciers larger than 0.05 km 2 . The large ice cap in the north together with surrounding local glaciers consists of 589 glaciers with an area of 20,784.4 ± 153.8 km 2 whereas the southern part has an area of 1612.6 ± 96 km 2 with 885 entities (Figure 1). From the entire sample, 41 glaciers are identified as marine terminating with an area of 16,063.7 ± 118.8 km 2 . Plotting the area and number of glaciers per size class reveals a distribution similar to other regions, like Greenland or Patagonia [40,41], but more extreme (Figure 5a). A high number of glaciers < 0.5 km 2 (909 or 62%) cover only 156 ± 1.7 km 2 (or 0.7%) of the total area and 49 glaciers (3.3%) > 100 km 2 have an area of 18,724 ± 2059 km 2 covering 84% of the total area. Small glaciers are mostly found in the south, whereas most of the large glaciers are part of the dominant ice cap in the north. The mean glacier size for the entire region is 15.2 km 2 whereas for the southern part it is only 1.8 km 2 .
The area-elevation distribution is depicted in Figure 5b for the whole region and the southern part of NVZ in per-cent. Glaciers are found from sea level to an elevation of more than 1100 m. The mean elevation is 596 m for all glacier entities (528 m for the northern part) and 641 m for those in the southern part (571 m for the overall glacierized surface). Most of the ice coverage can be found at an elevation interval of 600 to 700 m (whole island) and 500 m (southern part), respectively.
The mean, mid-point or median elevation of a land terminating glacier can be used as a proxy for its balanced-budget equilibrium line altitude ELA0 [42]. Its values are thus reflecting climatic conditions in a region, in particular precipitation amounts (e.g., [43]). As many glaciers in the northern part are marine or lake terminating and glacier separation is to some extent subjective for flat ice caps, we only display mean elevation of glaciers in the southern part in Figure 6. The spatial analysis of this sample reveals no latitudinal but some longitudinal dependencies, i.e., an increase in mean elevation with distance from the western coast. Low mean elevations are visible for glaciers close to the west coast and these are steadily increasing inwards. The highest elevations can be found on the main mountain crest.

Glacier Characteristics
The new glacier inventory for NVZ has an area of 22,379.0 ± 246.2 km 2 with 1474 glacier entities when considering only glaciers larger than 0.05 km 2 . The large ice cap in the north together with surrounding local glaciers consists of 589 glaciers with an area of 20,784.4 ± 153.8 km 2 whereas the southern part has an area of 1612.6 ± 96 km 2 with 885 entities (Figure 1). From the entire sample, 41 glaciers are identified as marine terminating with an area of 16,063.7 ± 118.8 km 2 . Plotting the area and number of glaciers per size class reveals a distribution similar to other regions, like Greenland or Patagonia [40,41], but more extreme (Figure 5a). A high number of glaciers <0.5 km 2 (909 or 62%) cover only 156 ± 1.7 km 2 (or 0.7%) of the total area and 49 glaciers (3.3%) >100 km 2 have an area of 18,724 ± 2059 km 2 covering 84% of the total area. Small glaciers are mostly found in the south, whereas most of the large glaciers are part of the dominant ice cap in the north. The mean glacier size for the entire region is 15.2 km 2 whereas for the southern part it is only 1.8 km 2 .
The area-elevation distribution is depicted in Figure 5b for the whole region and the southern part of NVZ in per-cent. Glaciers are found from sea level to an elevation of more than 1100 m. The mean elevation is 596 m for all glacier entities (528 m for the northern part) and 641 m for those in the southern part (571 m for the overall glacierized surface). Most of the ice coverage can be found at an elevation interval of 600 to 700 m (whole island) and 500 m (southern part), respectively.
The mean, mid-point or median elevation of a land terminating glacier can be used as a proxy for its balanced-budget equilibrium line altitude ELA 0 [42]. Its values are thus reflecting climatic conditions in a region, in particular precipitation amounts (e.g., [43]). As many glaciers in the northern part are marine or lake terminating and glacier separation is to some extent subjective for flat ice caps, we only display mean elevation of glaciers in the southern part in Figure 6. The spatial analysis of this sample reveals no latitudinal but some longitudinal dependencies, i.e., an increase in mean elevation with distance from the western coast. Low mean elevations are visible for glaciers close to the west coast and these are steadily increasing inwards. The highest elevations can be found on the main mountain crest.  By area, most of the glaciers are oriented towards the south-east (>35%), north-west (20%) and north (~15%) (Figures 7a and b). Mostly small glaciers are found in the other sectors. When focusing only on the southern part, most of the glacier area and numbers of individual glaciers are exposed to the south-east (>20%), east (>19%) and north-east (~13%).   By area, most of the glaciers are oriented towards the south-east (>35%), north-west (20%) and north (~15%) (Figures 7a and b). Mostly small glaciers are found in the other sectors. When focusing only on the southern part, most of the glacier area and numbers of individual glaciers are exposed to the south-east (>20%), east (>19%) and north-east (~13%).  By area, most of the glaciers are oriented towards the south-east (>35%), north-west (20%) and north (~15%) (Figure 7a,b ). Mostly small glaciers are found in the other sectors. When focusing only on the southern part, most of the glacier area and numbers of individual glaciers are exposed to the south-east (>20%), east (>19%) and north-east (~13%).  By area, most of the glaciers are oriented towards the south-east (>35%), north-west (20%) and north (~15%) (Figures 7a and b). Mostly small glaciers are found in the other sectors. When focusing only on the southern part, most of the glacier area and numbers of individual glaciers are exposed to the south-east (>20%), east (>19%) and north-east (~13%).  Plotting mean elevation versus glacier aspect (Figure 7b) reveals about 100 m lower elevations for glaciers facing east than facing southwest for glaciers on the southern part of NVZ. Glaciers facing north to east have lower mean elevations than glaciers facing other aspect directions.
A scatter plot displaying minimum and maximum elevation versus glacier size for individual glaciers is depicted in Figure 8a. The plot indicates for the whole island and the southern part a high variability of values for small glaciers, whereas for glaciers larger than~5 km 2 a dependence of size on maximum elevation and minimum elevation is visible. Most of the largest glaciers, however, are marine terminating so that no dependence on size exists. A comparison of mean slope vs. glacier size is depicted in Figure 8b. It shows that large glaciers have rather low slope values, whereas small glaciers tend to cover the full range of slopes.
Remote Sens. 2017, 9, 11 11 of 18 Plotting mean elevation versus glacier aspect (Figure 7b) reveals about 100 m lower elevations for glaciers facing east than facing southwest for glaciers on the southern part of NVZ. Glaciers facing north to east have lower mean elevations than glaciers facing other aspect directions.
A scatter plot displaying minimum and maximum elevation versus glacier size for individual glaciers is depicted in Figure 8a. The plot indicates for the whole island and the southern part a high variability of values for small glaciers, whereas for glaciers larger than ~5 km 2 a dependence of size on maximum elevation and minimum elevation is visible. Most of the largest glaciers, however, are marine terminating so that no dependence on size exists. A comparison of mean slope vs. glacier size is depicted in Figure 8b. It shows that large glaciers have rather low slope values, whereas small glaciers tend to cover the full range of slopes.

Glacier Changes
The glacier area changes for the whole island and glaciers larger than 0.05 km 2 between the RGI (22,128.2; 1604 km 2 ) and the new inventory (22,379; 1612 km 2 ) cannot be determined due to differences in the interpretation of glacier areas and unknown acquisition dates. We have thus only calculated the total area in the RGI 5.0 for the southern part of the island for glaciers larger 0.5 km 2 (1543.8 km 2 ) as here we know the acquisition year (~2001). Our value for the 2016 NVZ inventory is 1491.3 ± 87.9 km 2 , giving an overall area decrease of −52.5 ± 4.5 km 2 or −3.4% in 15 years (−0.22% a −1 ). As the RGI glacier outlines were rather conservative (debris cover and ice in steep terrain is often not included) and our new inventory is the opposite, not only losses but also glacier area increases are visible in Figure 9. Analyzing the changes glacier by glacier gives area increases of about 25 km 2 and area losses of −126.6 km 2 or −7.9% in 15 years (−0.52% a −1 ). Figure 9 also shows that values of area loss and area gain increase towards smaller glaciers and become very similar when larger glaciers are considered. The area gains are likely due to the different interpretation of glacier entities (see example in Figure 10) rather than to glacier growth.

Glacier Changes
The glacier area changes for the whole island and glaciers larger than 0.05 km 2 between the RGI (22,128.2; 1604 km 2 ) and the new inventory (22,379; 1612 km 2 ) cannot be determined due to differences in the interpretation of glacier areas and unknown acquisition dates. We have thus only calculated the total area in the RGI 5.0 for the southern part of the island for glaciers larger 0.5 km 2 (1543.8 km 2 ) as here we know the acquisition year (~2001). Our value for the 2016 NVZ inventory is 1491.3 ± 87.9 km 2 , giving an overall area decrease of −52.5 ± 4.5 km 2 or −3.4% in 15 years (−0.22% a −1 ). As the RGI glacier outlines were rather conservative (debris cover and ice in steep terrain is often not included) and our new inventory is the opposite, not only losses but also glacier area increases are visible in Figure 9. Analyzing the changes glacier by glacier gives area increases of about 25 km 2 and area losses of −126.6 km 2 or −7.9% in 15 years (−0.52% a −1 ). Figure 9 also shows that values of area loss and area gain increase towards smaller glaciers and become very similar when larger glaciers are considered. The area gains are likely due to the different interpretation of glacier entities (see example in Figure 10) rather than to glacier growth.
included) and our new inventory is the opposite, not only losses but also glacier area increases are visible in Figure 9. Analyzing the changes glacier by glacier gives area increases of about 25 km 2 and area losses of −126.6 km 2 or −7.9% in 15 years (−0.52% a −1 ). Figure 9 also shows that values of area loss and area gain increase towards smaller glaciers and become very similar when larger glaciers are considered. The area gains are likely due to the different interpretation of glacier entities (see example in Figure 10) rather than to glacier growth.  Remote Sens. 2017, 9,11 12 of 18 Figure 10. The problem of area increase and decrease exemplified for a small region in the southern part of Novaya Zemlya. Positive area changes are due to mostly the semi-automatic mapping method used in this study. Negative area changes are due to glacier shrinkage through time.

The New Inventory
The glacier inventory presented here has several improvements compared to the version available in the RGI 5.0. Firstly it covers the entire region, with glacier outlines, the DEM and derived topographic parameters all referring to a consistent period of four years. Secondly, all glaciers now have a time stamp, which is important for modelling their past and/or future evolution [44]. For two scenes we needed, however, to consult two Landsat 5 scenes from 1998 to get coverage free of seasonal snow. The latter is mandatory for glacier mapping, in particular in Arctic regions with their abundant perennial snowfields (e.g., [40]). New satellite missions, such as Sentinel 2A/B, will likely improve the situation in the future, as images are acquired much more frequently and the higher spatial resolution will improve identification of seasonal snow. Thirdly, we have improved the location of drainage divides using the new fused DEM in combination with fringes from SAR sensors. Further improvements are expected from more complete versions of the ArcticDEM or by filling of the still existing data voids with the TanDEM-X DEM or IDEM instead of the noisy ASTER GDEM V2. The latter was not possible yet as the already available IDEM was obtained for a different purpose. Figure 10. The problem of area increase and decrease exemplified for a small region in the southern part of Novaya Zemlya. Positive area changes are due to mostly the semi-automatic mapping method used in this study. Negative area changes are due to glacier shrinkage through time.

The New Inventory
The glacier inventory presented here has several improvements compared to the version available in the RGI 5.0. Firstly it covers the entire region, with glacier outlines, the DEM and derived topographic parameters all referring to a consistent period of four years. Secondly, all glaciers now have a time stamp, which is important for modelling their past and/or future evolution [44]. For two scenes we needed, however, to consult two Landsat 5 scenes from 1998 to get coverage free of seasonal snow. The latter is mandatory for glacier mapping, in particular in Arctic regions with their abundant perennial snowfields (e.g., [40]). New satellite missions, such as Sentinel 2A/B, will likely improve the situation in the future, as images are acquired much more frequently and the higher spatial resolution will improve identification of seasonal snow. Thirdly, we have improved the location of drainage divides using the new fused DEM in combination with fringes from SAR sensors. Further improvements are expected from more complete versions of the ArcticDEM or by filling of the still existing data voids with the TanDEM-X DEM or IDEM instead of the noisy ASTER GDEM V2. The latter was not possible yet as the already available IDEM was obtained for a different purpose.

Interpretation of the Results
The glacier statistics presented above are largely dependent on the topography of the islands of NVZ. This is particularly visible for the area-aspect distribution of the glaciers. More than 35% of all glaciers are oriented towards the southeast. This is different from other regions, because one would expect the larger glaciers to the North or on the windward rather than on the leeward side of a mountain range. To a certain degree this distribution is also visible in the southern part, but here we also found a gradient in mean elevation from the Barents Sea to the Kara Sea ( Figure 6), indicating that precipitation is decreasing from the coasts inwards which was also detected in other maritime regions such Greenland or Alaska [21,40,45]. For the marine terminating glaciers, the concept of mean or median elevation as a proxy for a balanced-budget ELA does not really apply so the given mean elevation should not be interpreted in this way. Glaciers that reach down to the lowest elevations are typically also the largest glaciers (Figure 8a) with the smallest mean slopes (Figure 8b). However, local topographic conditions (shading), excess nourishment or marine terminating tongues can cause a deviation from this general scheme. Despite these factors, the lower mean elevation of the land-terminating glaciers in the south at least indicates that precipitation is higher there.

Comparison with Other Studies
The calculated ice coverage of NVZ (22,379 ± 335.6 km 2 ) can be compared quantitatively with three other studies. Grant et al., 2009 [17] presented a number of~22,984 km 2 which is an estimate based on [46] and thus representative for the mid 1990s. Moholdt et al., 2012 [13] gave~22,100 km 2 using manual digitalization of SPOT5 and Landsat imagery (2000-2010), whereas [19] reports 24,400 km 2 of ice cover based on aerial photography for the years around 1950. Regarding the number of glaciers, 480 glacier entities are found in the RGI dataset [17], report 690 (based on the study of [19]) and we mapped 1474, largely due to the lower size threshold.
A qualitative and quantitative comparison can only be made with [13] due to its digital availability. In this case, we obtain for the whole island a slightly larger glacier area. Our larger value might also result from considering glaciers with sizes between 0.05 and 0.25 km 2 and possibly due to a larger area fraction mapped at high altitudes. This is also visible in the derived area changes for the southern part. As shown in Figure 9, area losses and gains were derived in particular for the smaller glaciers, whereas for the larger ones area losses are more similar. This demonstrates the principle difficulties when calculating area changes from datasets with a different origin and that possible differences in interpretation can result in larger area changes than in reality (cf. also [47]).

DEM Quality
Whereas optical satellite imagery (Landsat and Sentinel 2A/B) are increasingly available for Arctic regions, precise DEMs were lacking until very recently. The SRTM DEM is only available up to a latitude of 60 • N [48] and the (optical) ASTER GDEM V2 has a limited quality over the featureless terrain of an ice cap and is missing an exact time stamp. This situation has recently improved considerably with the high-resolution ArcticDEM now being available [18] and it will further improve with the release of the 12 m resolution TanDEM-X DEM [49].
As we have used a fused product (Arctic DEM void filled with ASTER GDEM V2), to determine topographic attributes of glaciers, we co-registered both DEMs. Additionally, we performed a cross-comparison of both DEMs on glacierized terrain with the raw DEMs (no co-registration). For this purpose, we selected two randomly selected glaciers (glacier a = southern part, glacier b = northern part; cf. Figure 1) and calculated zonal statistics of elevation values (minimum, maximum, mean and standard deviation) from both DEMs ( Table 4). The differences in these parameters for glacierized terrain are comparably small for both DEMs. The differences in maximum elevation of glacier (a) is 15 m; whereas for glacier (b) it is only about −2 m. The opposite is detectable for minimum elevation. This impacts on the range value derived for both glaciers having a +12 and −17 m difference. As both DEMs are acquired at different points in time, the differences in the mean (16.5 m and 3.6 m) is partly due to the missing co-registration but also due to real glacier elevation change.
In the latter case this would indicate that glaciers in the southern part are loosing mass, whereas in the northern part the situation might be more stable. This assumption is in contrast to other studies [13,50] that have shown that glaciers of NZV are losing mass and retreating. In general, the standard deviation is always slightly higher in the ArcticDEM than in the ASTER GDEM V2. Overall, we think that the combination of these two DEMs is acceptable for the application in this study.
Both the ArcticDEM and the ASTER GDEM V2 have gaps in NVZ (we here only consider gaps within the NVZ glacier inventory) because both DEMs are derived from optical imagery. These gaps may appear where cloud cover or fog obscured the ground or where unfrozen water bodies, or other regions of low radiometric contrast (shadow) exist, so that pixel correlation cannot be resolved by the software [51,52]. A statistical analysis of the gaps reveals that the area of the gaps in both DEMs is nearly the same: The ArcticDEM (ASTER GDEM V2) has an area of 1134.9 km 2 (1100.0 km 2 ). Differences are apparent in mean elevations: The gaps of the Arctic DEM have a mean elevation of 509 m a.s.l., whereas the gaps of the ASTER GDEM V2 are located at around 803 m a.s.l. We attribute the better performance at higher elevations to the better radiometric resolution of the panchromatic bands of the Worldview-1, 2, and 3 and GeoEye-1 satellite sensors of which most of the ArcticDEM data was generated from [18]. The differences in mean elevation of the DEM gaps also indicate that drainage divides on the large ice cap mostly stem from the better quality ArcticDEM.

Uncertainties and Limitations
The accuracy of the NVZ glacier inventory is similar to or even better than for other studies where the semi-automated mapping with the help of a ratio image in Arctic regions has been used (e.g., [40]). The small uncertainty value for glacier outlines of about 1.1% in the mean over all glaciers (obtained with the buffer method, [53]) is a consequence of the large mean glacier size (in absolute terms). We further observed that the uncertainty increases towards smaller glaciers as a more detailed study has recently shown [31] (cf. Figure 10). The uncertainty would be even higher if more manual corrections (e.g., for debris cover) would have to be performed.
Perennial snowfields, however, impact on the result, as they should be excluded in a glacier inventory but are difficult to distinguish from niche glaciers in topographic depressions that can be included. But this is a matter of definition, as in most cases also perennial snowfields might include ice that is not flowing. Moreover, when abundant seasonal snow is present at the same time such as in footprints D and G from 2015 (Table 1), the issue becomes difficult. By applying a size threshold of 0.05 km 2 we likely excluded many of them, but some of the larger entities are included and it can be a matter of preference (and intended use) if these need to be included or not. For scenes D and G the scenes from 1998 clearly helped in identifying and removing as the 2015 scenes have only been used for updates of terminus positions. Finally, by providing our dataset to GLIMS, we hope that a consensus on the interpretation can be reached at some point. As already mentioned in Section 5.2, the largest uncertainty in our dataset is possibly related to wrongly classified perennial snowfields.
Using fringes from SAR interferometry for correction of drainage divides in flat topography was very helpful to refine locations. However, a good temporal match is required as small changes in surface elevation can change the position of drainage divides in flat terrain considerably (e.g., [45]). As the TanDEM-X IDEM used here for geocoding of the SAR images is temporarily close to the ArcticDEM, the corrections applied should be consistent. However, corrections can only be applied when coherence is present which only applies to the slow flowing upper regions of the ice cap. The correction of the divides with the help of fringe data is, however, sometimes also a matter of discussion due to different interpretation of the operators. A detailed check by a second operator is highly recommended for a more robust result.
Splitting glacier complexes into entities is supported by automated watershed analysis but is challenging during selection and manual editing. It also depends on the purpose of the inventory, resulting in a likely smaller number of glaciers for glaciological versus hydrological applications. The combination of the most recent DEM together with fringes derived from radar interferometry likely provided the best possible dataset, but our divides refer to hydrological applications. For glaciological applications many of the divides separating ice caps have to be merged. The provided glacier size distribution would be very different in the latter case. As this dataset will be freely available from the GLIMS database, each user can modify the divides for the intended purpose.
Finally, we also analysed the area difference due to re-projection between different UTM zones, the area-preserving projection (North Pole Lambert Azimuthal Equal Area projection with D WGS 1984 datum) and the polar stereographic projection. Calculating the area for the entire region gives with the Lambert projection an area of 22,390.13 ± 246.1 km 2 and with the polar stereographic projection 21,747.1 ± 239.21 km 2 . This gives an area difference of~11 ± 0.12 km 2 (0.05%) with the former and 32 ± 6.9 km 2 (or about 3%) with the latter projection. This confirms that calculating areas for small regions like NVZ with maximal 3 UTM zones (40,41,42) only introduces uncertainties smaller than other uncertainties described in this paper. In contrast, the value derived from the Polar Stereographic projection is twice as high as the calculated uncertainty. We therefore suggest not calculating area values with UTM projections for regions stretching over more than 3 UTM zones (e.g., Greenland) and not calculating area from a dataset in Polar Stereographic projection.

Conclusions
This study has exemplified how recently available datasets can improve knowledge about glacier distribution in the high Arctic. Combining multi-temporal Landsat data with the new ArcticDEM resulted in a high-quality inventory for the glaciers in Novaya Zemlya within a temporal range of 4 years. Glaciers were mapped semi-automatically using the well-established band ratio method and drainage basins were computed with a semi-automated watershed analysis from the ArcticDEM and data voids filled by the ASTER GDEM V2. The location of the drainage basins was further improved using interferometric SAR data (fringes). We found a total glacier area of 22,379.0 ± 246.2 km 2 with 1474 glacier entities larger than 0.05 km 2 . The region is dominated by the large ice cap in the north and surrounding local glaciers consisting of 589 glaciers with an area of 20,784.4 ± 153.8 km 2 , whereas the southern part has 885 glaciers covering 1612.6 ± 96 km 2 with 885 entities. From the entire sample, 41 glaciers are identified as marine terminating and most of the ice is located at elevations from 600 to 700 m a.s.l. (southern part 500 m a.s.l.).
On the methodological side, we found that glacier definition (the way glaciers are mapped, e.g., including or excluding perennial snow) is crucial for deriving glacier changes when using external datasets. In the case different source data are combined, a careful evaluation is required and maybe a restriction to a common area, using the drainage basins of one dataset as a mask (the older one in the case of dominant retreat). Our interpretation of perennial snowfields is certainly open for discussion; a more strict interpretation is possible. We could only obtain the four-year temporal frame by integrating two scenes free of seasonal snow from 1998 and updating of terminus positions with the 2015 scenes. This was possible here as only few glacier tongues had to be updated, but it might not work elsewhere. Finally, despite the obvious quality differences between the ArcticDEM and the ASTER GDEM V2 we only found minor differences in the statistics of topographic parameters on and off-glaciers. This indicates that differences tend to average out over larger regions and that the ASTER GDEM V2 can be well used after co-registration for this purpose (DEM gap filling) despite its quality issues.