Updated inventory of glacier ice in New Zealand based on 2016 satellite imagery

Abstract

The only complete inventory of New Zealand glaciers was based on aerial photography starting in 1978. While there have been partial updates using 2002 and 2009 satellite data, most glaciers are still represented by the 1978 outlines in contemporary global glacier databases. The objective of this project is to establish an updated glacier inventory for New Zealand. We have used Landsat 8 OLI satellite imagery from February and March 2016 for delineating clean glaciers using a semi-automatic band ratio method and debris-covered glaciers using a maximum likelihood classification. The outlines have been checked against Sentinel-2 MSI data, which have a higher resolution. Manual post processing was necessary due to misclassifications (e.g. lakes, clouds), mapping in shadowed areas, and combining the clean and debris-covered parts into single glaciers. New Zealand glaciers cover an area of 794 ± 34 km² in 2016 with a debris-covered area of 10%. Of the 2918 glaciers, seven glaciers are >10 km² while 71% is <0.1 km². The debris cover on those largest glaciers is >40%. Only 15 glaciers are located on the North Island. For a selection of glaciers, we were able to calculate the area reduction between the 1978 and 2016 inventories.

1. Introduction

The cryosphere plays a major role in Earth's climate system (e.g. Vaughan and others, 2013). Mountain glaciers are recognized as key indicators of climate change (e.g. Haeberli, 2004) and are retreating and losing mass at historically unprecedented rates (Zemp and others, 2015). Snow and ice cover for the Southern Alps/Kā Tiritiri o te Moana (further called Southern Alps) of New Zealand (NZ) provide an important natural resource that supports power generation, primary production and water resources. Significant glacier shrinkage, including the formation of proglacial lakes, has occurred in recent decades, e.g. at Hooker Glacier (Kirkbride and Warren, 1999; Robertson and others, 2013; Purdie and others, 2014). Mackintosh and others (2017) highlighted the existence of strong climate teleconnections between high latitude climate modes and the mid latitudes of the Southern Hemisphere, emphasizing the global context of NZ glacier variations.

The last complete glacier inventory of NZ from 1978 (South Island/Te Waipounamu) and 1988 (North Island/Te Ika-a-Māui) (further called South Island and North Island, respectively) was constructed manually from oblique aerial photographs, taken at the end of summer, and geodetic maps (Chinn, 2001). It has a positional uncertainty of ±30 m corresponding to the accuracy of the topographic base maps (Gjermundsen and others, 2011). In total, 3144 glaciers larger than 0.01 km² were counted, covering a total area of 1158 km² and representing ~53 km³ of ice (Chinn, 2001). The inventory has been partly updated by Gjermundsen and others (2011) for the year 2002 (~40% of total area) and by Sirguey and More (2010) for the year 2009 (~32% of total area), both using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data. In comparison to 1978, glacier area was reduced by ~17% in 2002 (0.71% a⁻¹; Gjermundsen and others, 2011) and by ~11% in 2009 (0.35% a⁻¹) for these limited areas in the central Southern Alps. The lower area reduction until 2009 probably results from different image (mis-)interpretation and glacier advance between the 1980s and 2008 (Chinn and others, 2005; Purdie and others, 2014). As there are no glacier outlines for the 2002 update and no publication for the 2009 update, they cannot be compared directly and their difference in area reduction remains uncertain.

To help in the compilation of glacier inventories, remote-sensing techniques and automated computer processing have increasingly been applied since the 1980s (Ohmura, 2009). The advantage of these methods is the ability to map a large number of glaciers even in high elevation and remote areas (Kargel and others, 2005; Bhambri and Bolch, 2009). A regular update of glacier inventories is essential at time intervals of a few decades, i.e. equivalent to typical response times of mountain glaciers, preferably by using remote sensing (Haeberli and others, 2002). The response time of glaciers in NZ has a mean value of ~10 years (Hoelzle and others, 2007), with fast reacting glaciers such as Te Moeka o Tūawe/Fox (further called Fox) and Kā Roimata o Hine Hukatere/Franz Josef (further called Franz Josef) Glaciers responding within 5–8 years and large and low-gradient valley glaciers such as Haupapa/Tasman (further called Tasman) Glacier within the order of 100 years (Chinn, 2001). There have been considerable fluctuations in glacier extent since 1978 (advances and re-retreat) (Chinn and others, 2005; Purdie and others, 2014), mostly at glaciers with a short response time, and the current extent is unknown. Hence, the existing inventory from 1978 does not correspond to the contemporary state of glaciers nor to internationally established recommendations for glacier monitoring (GCOS, 2016).

By mid-2013, the Global Land Ice Measurements from Space (GLIMS) database contained ~58% of the global glacier area outside the ice sheets (Pfeffer and others, 2014) and is now globally complete with multitemporal coverage of nearly 100 000 glaciers (http://glims.colorado.edu/glacierdata/). The Randolph Glacier Inventory (RGI) was initiated to create a globally complete dataset for IPCC AR5 based on existing outlines in GLIMS and a dedicated community effort to provide outlines for the missing regions. The NZ glaciers are until now represented in these global databases as outlines from the year 1978, combined with the 2009 update in GLIMS.

Glacier monitoring data are strongly biased toward the Northern Hemisphere (Chinn and others, 2012), especially concerning multi-year coverage of glacier outlines. Therefore, we created a new and consistent inventory of NZ glaciers based on satellite imagery from 2016. Updated information about glaciers in NZ will constitute a considerable benefit to the global glacier and climate monitoring networks (Chinn and others, 2012).

2. Study area

New Zealand is located approximately between 34°–46°S and 166°–178°E. It mainly consists of two larger islands (the North and South Islands) with an oceanic maritime climate dominated by a strong west-east precipitation gradient (Chinn and others, 2005; WGMS, 2008). The high (up to 3724 m a.s.l.) but narrow (~50 km) mountain ranges of the Southern Alps stretch northeast to southwest and parallel to the west coast of the South Island, acting as a barrier oriented perpendicular to the prevailing westerly winds. Mean annual precipitation is 3 m a⁻¹ along the western coastal plains rising to 15 m a⁻¹ close to the Main Divide (central watershed between east and west) (Chinn, 2000). The glaciers at and immediately east of the Main Divide benefit from frequent spillovers (Chater and Sturman, 1998), but precipitation in the eastern ranges drops to ~1 m a⁻¹ with a steep precipitation gradient (Sturman and Wanner, 2001). The abundant precipitation results in numerous late season and perennial snowfields that are difficult to distinguish from glaciers.

One percent of global mountain glacier area is located in NZ (Chinn, 2001). The glaciers of NZ are mostly distributed within the Southern Alps and adjacent ranges of the South Island (42.0°–45.9°S, 167.3°–173.8°E) (Gjermundsen and others, 2011; Chinn and others, 2012) (Fig. 1). Summits range from 1850 m a.s.l. in Fiordland (at the southern end of the South Island) to over 3000 m a.s.l. in the central Southern Alps, and to 2000 m a.s.l. in the north (Chinn and others, 2005) with the highest peak at 3724 m a.s.l. (Aoraki/Mt Cook) (https://www.openstreetmap.org/). On the North Island, Mt Ruapehu (2797 m a.s.l.), a stratovolcano, has a few small glaciers (~0.4% of the total glacier area).

Fig. 1. Location of the two study regions in New Zealand (red squares in the upper left inset) and glacier distribution on the South Island (main image) and North Island (upper center inset) shown on a shaded relief of the 15 m DEM by Columbus and others (2011). New Zealand outline: ESRI shapefile).

After the first mapping of glaciers in NZ in the 1860s–1890s, the glaciers reduced their area and volume until the mid-1970s (Chinn and others, 2008; Salinger and others, 2008). Since 1978, annual flights by the National Institute of Water and Atmospheric Research in small aircraft, referred to as the end-of-summer-snowline (EOSS) measurements, have been used to photograph and measure the annual snow line as an estimate for the equilibrium line altitude (ELA) at 50 index glaciers, and terminus position at up to 126 glaciers. However, the terminus positions are not photographed every year depending on the flight plan and cloud cover. These data and ground measurements showed a considerable glacier advance that started during the early 1980s and lasted until about 2008 (Chinn and others, 2005; Purdie and others, 2014). This advance was most prominent at the large but steep Fox and Franz Josef Glaciers, but at least 58 of the measured glaciers advanced during this period. The cause of the advance was lower temperatures with a contribution from increased precipitation (Mackintosh and others, 2017). This advance was not recorded at those large debris-covered valley glaciers with proglacial lakes (Winkler and others, 2010) as a result of their long response times and the insulating effect of the debris cover (Chinn and others, 2012). These debris-covered glaciers have generally experienced mass loss since the late-1800s evident by a steady lowering of ice-surface levels (downwasting) accompanied by small or no changes in terminus position for many decades (Chinn and others, 2005, 2012; Winkler and others, 2010). Between the 1970s and the 1990s, lowering ice levels of these protracted response glaciers reached the turning point for the development of proglacial lakes and a thereby induced acceleration of ice loss in their lower parts (Chinn and others, 2012).

Glaciers in New Zealand are situated in a maritime climate with abundant precipitation to the west of the Main Divide and much drier conditions to the east of it (Chinn, 1999; Hoelzle and others, 2007; Gjermundsen and others, 2011). They exhibit a high-alpine character with mostly individual glaciers (Chinn and others, 2005). Many calve over steep slopes, partly resulting in (disconnected) regenerated glaciers or – due to the abundant snowfall – avalanche deposits in valley floors and topographic depressions. Such features are difficult to handle in glacier inventories as their separation gives extreme values for topographic information (i.e. geometric/topographic information from the glaciers is less readily relatable to climate), and the presence of flowing ice under snow and ice deposits is hard to determine. The largest valley glaciers with heavy debris cover and proglacial lakes are situated in the central Southern Alps (Chinn and others, 2005), close to Aoraki/Mt Cook. Most of the NZ glaciers are of small (<0.32 km²) or medium size (0.32−2.56 km²; Chinn, 2001) and show short response time characteristics such as large mass turnover and high terminus ablation rates (Chinn and others, 2012). Many NZ glaciers are debris-covered (WGMS, 2008) making it difficult to delineate glacier extents. Previous estimates indicate debris cover comprises ~8% of the total glacier area in the central Southern Alps (Anderson and Mackintosh, 2012) and 25% of the area of the large valley glaciers (Chinn, 1996).

3. Material and methods

3.1 Input data

The last inventory of 1978 for the South Island and 1988 for the North Island by Chinn (2001) gave a detailed and thorough analysis of NZ glaciers and acts as a valuable database. Twenty basic parameters were mapped for each glacier, including identification (ID) and position information, dimension measurements (e.g. area, length), aspect and type. A particularly important parameter, often not included in other glacier inventories, was the measurement of the debris-covered area. This database is available as an unpublished hard copy, but the digitally available GLIMS outlines from 1978 do not correspond to these original outlines as a different person digitized them. Therefore, some random errors were included in the 1978 inventory that are now inherent in the current GLIMS and RGI databases. Hence, we only used selected features such as GLIMS-ID and drainage divides from this dataset and did not perform a complete comparison or change assessment. Instead, we performed a change analysis based on the hard copy data manually digitized for a small selection of glaciers in the central glacier area around Aoraki/Mt Cook that covers a substantial portion (32% of glacier area) of the total inventory.

Compilation of glacier inventories from automated classifications is an established procedure (Frey and others, 2012). Clean and slightly dirty glacier ice can be mapped automatically due to the different spectral properties of ice and snow in the shortwave infrared (SWIR) compared to other surface types (Paul and others, 2015) and their high reflectance in the visible and near infrared (Racoviteanu and others, 2009; Gjermundsen and others, 2011). Minimum mapping area of glaciers is usually 0.01 km² to exclude most of the often smaller seasonal and perennial snowfields (Paul and others, 2009; Pfeffer and others, 2014). Glaciers in NZ are small (more than 50% of the glaciers in the 1978 inventory are smaller than 0.08 km² (Chinn, 2001)), meaning the satellite images need to have a sufficient spatial resolution to properly identify them (Leigh and others, 2019).

Landsat 8 OLI and Sentinel-2 MSI satellite images from 2016 were chosen to create a new glacier inventory for NZ (Table 1). They are generally captured around 10.30 am local time over NZ (i.e. NZDT), the acquisitions being ~15 min apart. In summer, the sun is sufficiently high in the sky at that time so that minimal shadowing occurs, i.e. images captured in the months from November to February were optimally suited. However, this needed to be balanced against the often abundant seasonal snow. Usually seasonal snow will not have melted until late in the summer, or even until early autumn, i.e. late February to April. Therefore, images are preferably taken at the end of the ablation period with a minimum amount of seasonal snow and no/sparse cloud coverage (Williams, 1986; Paul and others, 2009, 2013; Racoviteanu and others, 2009). Some care is required as early snowfall in autumn can partly remain in higher elevations. To minimize the effect of seasonal snow during image processing, it is useful to select an image from a negative mass-balance year (Gjermundsen and others, 2011). As additional requirements, the new inventory should be as recent as possible and all images should date from one single year.

Table 1. Specifications of ASTER, Landsat 8 OLI and Sentinel-2 MSI satellites (modified after Winsvold and others, 2016)

^a 10 days over New Zealand in 2015/2016.

On the basis of the above criteria, images from February to April 2016 were chosen, as they had sparse seasonal snow, limited shadowing and a nearly total coverage of all glaciers. Nonetheless, cloud cover always provides challenges in alpine areas and especially in the maritime climate of New Zealand. Several overlapping satellite images had to be combined to cover the total area not only to overcome cloud cover but also to merge the images with a minimum of snow in higher elevations (February) and the completely melted early autumn snowfall in lower elevations (April).

Landsat 8 has a 30 m resolution in the multispectral bands, which is comparably coarse but has a high radiometric resolution (Table 1) and is therefore suitable for these mapping purposes. The images chosen for glacier mapping are listed in Table 2 with the location of footprints in Figure 2.

Fig. 2. Footprints of the Landsat 8 OLI (brown) and Sentinel-2 MSI (green) images used for glacier mapping. The two red markers (Region A and B) show the locations used for the uncertainty assessment (see Fig. 4).

Table 2. Selected parameters of the applied Landsat 8 OLI/TIRS and Sentinel-2 MSI images (source Landsat parameters: https://landsatonaws.com, source Sentinel parameters: https://remotepixel.ca/) and a comment on snow and cloud conditions for each scene

We encountered no issues due to geolocation mismatch between the Landsat 8 and Sentinel-2 images. One Sentinel-2 image (Number 11 in Table 2 from 16.03.2016, 58GGR) was slightly shifted and, therefore, spatially corrected.

In addition to suitable satellite scenes, an appropriate DEM is required. We used the freely available NZSoSDEM v1.0 with a spatial resolution of 15 m from the School of Surveying, University of Otago, based on the Land Information New Zealand (LINZ) Topo250 topographic map series (Columbus and others, 2011).

3.2 Mapping glaciers

Clean and debris-covered ice were mapped semi-automatically using Landsat 8 images as Sentinel-2 data resulted in incomplete patchy glacier outlines with the thresholds we selected for the band ratio method. The band ratio approach was used for clean ice (ratio: red/SWIR), i.e. ratio band 4/band 6 with a threshold of 1.9 (mean value, varied for each Landsat image) and an additional threshold in the QA band (cloud confidence). For comparison, we also mapped the glacier area of one satellite image (central glacier area) with the Normalized Difference Snow Index. Both methods showed more or less the same result with a relative difference of 1.6%. Hence, we judged the band ratio as sufficient for mapping glaciers in this region. The results were 3 × 3 median filtered and the outlines were additionally smoothed using the PAEK method in ArcGIS.

Automated mapping of debris-covered glaciers is still a challenge. We mapped debris-covered ice semi-automatically using a maximum likelihood classification. This method is a supervised classification in which the user draws training sites to identify pixels within arbitrary classes. The whole image is then classified based on the likelihood that each pixel falls within a particular class. In this analysis, three inputs were used: (1) The multispectral Landsat 8 bands, (2) the thermal Landsat 8 bands (3) and a slope layer, which had been created from the DEM. The first step was to identify ~20 training sites for each of the four classes ‘Water’, ‘Bare Ice’, ‘Glacier Debris’ and ‘Miscellaneous Debris’. The next step was to carry out a principal components analysis (PCA) for analyzing correlation between bands and rotating them so that data obsolescence is highlighted. The PCA eigenvalues indicated that only the first three bands of the PCA contained significant information, which were stacked with the slope layer as input to the maximum likelihood classifier. The resulting glacier debris class was clipped to the GLIMS glacier extent (including a buffer area), as it was assumed that the debris-covered extent of the glaciers would not have grown or shifted significantly outside the previous extent since 1978. Upglacier migration of debris cover was allowed within the GLIMS extent. This reduced the number of misclassified pixels.

For the North Island, we used an existing glacier outline dataset created by manual digitization from Google Earth (Image source: Maxar Technologies; Acquisition date: 14 April 2016), which was guided by field observations of glacier margins made between 2012 and 2015 (Eaves, 2015). This dataset was edited using the methods below.

The result of this mapping was two datasets: one set with the clean-ice glacier parts and another set with the debris-covered parts.

3.3 Correction of raw glacier outlines

Manual post processing was necessary due to misclassifications (e.g. lakes, clouds), mapping in shadowed areas, and the connection of bare and debris-covered sections (Fig. 3). Additional input data for the post processing were mainly the Sentinel-2 MSI images (Fig. 2). As the band ratio mapping with Landsat and the maximum likelihood classification of debris cover was not satisfying, especially compared to the extent of the glaciers on the Sentinel images, the majority of the outlines was corrected manually. Therefore, the initial input served only as a basis for the resulting mapping. Hence, the final inventory is more or less based on manual editing and not on automated mapping. Over cloudy areas in the Landsat 8 as well as in the Sentinel-2 scenes (max. 1% of glacier area), orthophotos from LINZ (resolution: 0.75 m, date varying: 2004–2016), Google Earth images and the 1978/88 outlines from GLIMS (Bolch and Chinn, 2013) were used. These additional images served as a guide, where missing parts were completed by creating straight lines or following the curvature of the glacier outline. Because of the wide variety and mix of data, we do not set a precise date to our outlines but assign the date range February to April 2016 which covers almost all (>99%) of the outlines.

Fig. 3. Raw outlines from automatic mapping (green = clean ice, brown = debris cover) with final outlines (blue) for McKerrow (left) and Mueller Glaciers (right) (Background: Landsat 8 OLI from 13.02.2016).

3.4 Application of drainage divides and finalization of dataset

As a last step, the connected ice units were split by the glacier ice-divides derived from the 1978 inventory (RGI Consortium, 2017), and the GLIMS IDs were assigned to each ice area. We noted that several ice-divides in the 1978 inventory are in the wrong position but have opted to keep them for consistency with the earlier dataset. Several debris-covered areas within one glacier were summed up and treated as one area for the subsequent analysis (e.g. only one minimum elevation was calculated for statistics and not several). The ice-divides are sensible from a water resource perspective and have to be considered when changes of glacier extent are interpreted in climatic terms. Differences between perennial snow and glaciers are not visible from space as the former is ‘snow that persists on the ground year after year’ (https://nsidc.org/cryosphere/glossary/term/perennial-snow). Despite careful glacier mapping, perennial snowfields may be included in our glacier dataset, especially if they are connected to glacier areas. We follow the common approach of assuming that any identified ice body with an area >0.01 km² is classified as a glacier (Leigh and others, 2019). We assigned each glacier as ‘debris-covered’ if it exhibited any debris cover independent of its amount. The final dataset will be submitted to the GLIMS database. For this submission, we plan to adapt the incorrect ice-divides. This correction will be based on the DEM used in this study.

3.5 Assessment of topographic information

All topographic information, such as elevation, aspect and slope, was derived from the NZSoSDEM v1.0 and added as additional information to each glacier and to each debris-covered area.

3.6 Uncertainty assessment

For assessing the uncertainty due to mapping, we undertook a digitizing experiment. Two regions (Fig. 4, location in Fig. 2) were chosen to map the glaciers manually. These images contained clean ice glaciers, debris-covered glacier parts and glacier area in cast shadow. Three experts digitized the glacier area independently multiple times with a maximum total overlay of nine outlines for one glacier. No limit on glacier size was set. There were 55 glaciers on both images, and 39 were mapped more than once. In the following analysis, we included only these 39 glaciers, but disregarded the changing number of glacier outlines.

Fig. 4. (a) Region A and (b) region B (location see Fig. 2) of the digitizing experiment for the uncertainty assessment. Overlay of multiple glacier mappings by three different GIS operators who are experts in glacier mapping. (c) Statistics from the digitizing experiment concerning glacier area and glacier (image) properties.

The largest glacier had a mean area of 4.32 km² and the smallest a mean area of 0.003 km². As can be seen in Figure 4c, the largest digitizing variability (up to 44%) was recorded for glaciers <0.01 km² (that are excluded from the inventory). This high percentage is a consequence of a relatively consistent pixel digitizing error for features with very small areas. The std dev. for the smallest glaciers (<0.1 km²) included in the inventory was between 5 and 35%. Relative digitizing variability decreased with mean glacier size. Glaciers in cast shadow (black asterisks in Fig. 4c) had a similar mapping precision to unshadowed glaciers. The glacier with cast shadow and debris cover (red square in Fig. 4c) is the large glacier complex in Figure 4a with two debris-covered tongues. As this complex was partly mapped in several pieces and partly as one glacier, the evaluation was assessed as one glacier. The mapping deviation (~2%) is also low. Generally, the std dev. was 7.5% for the debris-covered glaciers (not area weighted). The area weighted std dev. for all glaciers >0.01 km² was 4.3%. This was added as the uncertainty value to the total inventory.

4. Results

In total, we map 2918 single glacier units with an area of 794 ± 34 km² in 2016 (Table 3). There are 141 debris-covered glaciers with a total 82 ± 4 km² of debris area (10% of total area). The elevation of all glaciers ranges from 420 to 3724 m a.s.l. with a mean of 1953 m a.s.l.

Table 3. Inventory parameters

^a Number of debris cover is the number of glaciers that have debris cover (not the number of debris-covered areas as there can be multiple debris-covered areas per glacier).

There are 15 glaciers in the inventory on the North Island (all on Mt Ruapehu) with a total area of 3.0 km². Seven of these glaciers are debris-covered with a debris area of 1.3 km² (47%). All glaciers are smaller than 1 km² with a mean of 0.2 km². Forty percent of the glaciers on the North Island are very small (<0.1 km²) and only one glacier is larger than 0.5 km². The minimum elevation is 2196 m a.s.l., maximum 2742 m a.s.l. and mean elevation is 2494 m a.s.l.

Most of the glaciers in NZ are small (Fig. 5). Of the 2918 glaciers mapped, 2064 (71%) are in the smallest size class <0.1 km², covering a total area of 68 km² (9%). Many glaciers (610) are also found in the next size class, 0.1–0.5 km² covering a total area of 132 km². Both classes together comprise 92% of the total glacier number, but only cover 25% of the glacier area. Ninety-six percent of the glaciers are <1 km² with 36% of the glacier area. Only a small amount (12 km²; 4%) of these glaciers is debris-covered. On the other hand, only seven glaciers are >10 km², and together cover an area of 206 km² (26%). The two largest size classes with 0.8% of the total glacier number cover 41% of the total glacier area and nearly all are debris-covered with a coverage of 16%.

Fig. 5. Area and number of glaciers according to size classes. White = clean glacier area, grey = debris cover. Note log y-scale for number of glaciers.

The minimum area of debris cover per glacier is 0.01 km² and the maximum 22.92 km² with a mean of 0.58 km² and a relative std dev. of 35%. There are 27 glaciers with more than one debris-covered part. Minimum elevation of the debris cover is 420 m a.s.l., maximum elevation 2702 m a.s.l. and mean elevation 1562 m a.s.l. There are 127 glaciers that have substantial debris cover, with at least 5% debris-covered area (Table 4).

Table 4. Inventory parameters of debris-covered glaciers with a threshold debris-covered area proportion of 5%

Prominent glaciers in NZ are Tasman, Fox and Franz Josef Glaciers as they are frequently visited by tourists. Tasman is, at 82.8 km², the largest and highest glacier in NZ (Fig. 1) and lies at the east side of Aoraki/Mt Cook. About 28% of its area is debris-covered, and it stretches from 3724 to 800 m a.s.l. This maximum elevation of Tasman is associated with a tributary flowing off the eastern flanks of Aoraki (Linda Glacier). The proglacial lake of the calving Tasman Glacier (known as Tasman Lake) has an area of 6.9 km². Fox and Franz Josef Glaciers are west of the Main Divide immediately adjacent to Tasman Glacier, i.e. connected in the highest part of their accumulation areas. They have areas of 32.5 and 32.0 km², respectively, and only a small amount of debris cover on their tongues (Fox: 0.43 km², Franz Josef: 0.20 km²). The minimum elevation of Fox Glacier (420 m a.s.l.) is the lowest of all NZ glaciers and this glacier reaches up to 3062 m a.s.l. Franz Josef Glacier stretches from 2870 to 593 m a.s.l. Neither glacier has a proglacial lake.

Figure 6 shows the minimum and maximum elevation for all glaciers. Within each size class, we calculated the mean for all glaciers and for glaciers with any debris cover. The mean values of the maximum elevation rise with increasing size class and vice versa for the minimum elevation, which implies a larger altitudinal range with increasing size class (see Table 3). The mean minimum and maximum elevation of the smallest glaciers with debris cover are nearly 500 m lower than those of all glaciers. The mean minimum elevation is always lower for the debris-covered glaciers while the mean maximum elevation gradually coincides with the mean maximum elevation of all glaciers. All glaciers >10 km² are partly debris-covered and, therefore, the mean values are identical. The mean elevation H_mean is in a small range (1911–1993 m a.s.l.) for all size classes of all glaciers (Table 3) and shows a large variability in case of the debris-covered glaciers with a threshold debris-covered area of 5% (Table 4).

Fig. 6. Minimum (blue) and maximum (green) elevation and their mean (solid line) for all glaciers according to size class. The mean of debris-covered glaciers is indicated by a dashed line. Note log scale on x-axis.

The mean glacier elevation distribution is depicted in Figure 7. All glaciers on the North Island belong to the highest categories with H _mean >2300 m a.s.l. The large Tasman (H _mean = 1884 m a.s.l.) and Murchison (H _mean = 1799 m a.s.l.) Glaciers dominate with their light reddish colors east of the Main Divide, and Fox (H _mean = 2119 m a.s.l.) and Franz Josef (H _mean = 2088 m a.s.l.) Glaciers with their strong reddish colors on the west of the Main Divide (Fig. 7c). Nevertheless, glaciers with higher and lower mean glacier elevations appear on both sides of the Main Divide (Fig. 8). The mean elevation of the western glaciers is mostly higher compared to the eastern ones and these glaciers do not appear in the northern and northwestern aspect sector.

Fig. 7. Mean glacier elevation for selected glacier areas (mean value at ~1950 m a.s.l.). Dotted-area represents debris cover.

Fig. 8. Mean aspect vs mean elevation for glaciers >0.1 km² with mean elevation of aspect sector in respect to location (east or west of the Main Divide).

The distribution of glacier area according to elevation shows a small areal amount in the lowest (<800 m a.s.l.) and highest elevation bins (>2800 m a.s.l.; Fig. 9). About 60% of the glacier area is located between 1800 and 2200 m a.s.l. and 81% between 1600 and 2400 m a.s.l. The maximum debris-covered area (~60%) is located between 1000 and 1400 m a.s.l. and 90% between 800 and 1600 m a.s.l., with none above 2800 m a.s.l.

Fig. 9. Distribution of area with elevation of the total glacier area (white) and the debris-covered area (dark grey) in 100 m elevation bins.

Most of the glacier area is located in the southern aspects, which includes both total glacier area and the debris-covered parts (Fig. 10). The debris-covered area is less frequently located in the northern aspects than the total glacier area (Fig. 10b).

Fig. 10. (a) Absolute (km²) and (b) relative (%) total glacier area (blue) and debris cover (brown) according to aspect. Note different scales and values that are given on the E-axis.

While Figure 10 shows the distribution of glacier area vs aspect based on the individual 15 × 15 m DEM pixels, we also compare the distribution of the total glacier area and debris cover vs mean aspect of the whole glacier areas. The debris-covered parts are also evenly distributed over all aspects as the total glaciers, especially for the areas <1 km². The larger areas are mainly located in the SW- to E-sections, the largest ones mainly in the south. In addition, we only compare the mean aspect vs the amount of debris cover in the ablation area (approximated by the mean glacier elevation) for glaciers >0.1 km² according to the location (east or west of the Main Divide). No ablation area with any amount of debris cover is located in the north or northwest and only sparsely in the west and northeast. The ablation areas are distributed over all other aspects with an accumulation around SW for the eastern glaciers and between S and SE for the western ones. No correlation was found between the amount of debris cover and mean aspect in the ablation area.

As the glacier data from the 1978 inventory from GLIMS/RGI are not completely consistent with the original, paper-based inventory, no area change assessment of the total area was possible. Nonetheless, we extracted 14 glaciers from the original data from 1978 (Chinn, 2001), and compared these with the respective glaciers from the 2016 inventory (Fig. 11). In 2016, all of these glaciers are partly debris-covered, with eight of them calving into proglacial lakes. This sample is not representative of all NZ glaciers and contains only the largest glaciers.

Fig. 11. Tasman, Fox and Franz Josef Glaciers with surrounding glacier area (for location see Fig. 1). Green outline = glacier area from 1978, blue area = clean ice glacier area from 2016, brown area = debris-covered glacier area from 2016 (background image: Sentinel-2 MSI from 12 February 2016).

All of these 14 large glaciers split into several parts with a minimum of four single glacier units and a maximum at Tasman of 26 single glacier units. Nevertheless, Tasman is a special case as it is a compound valley glacier with a number of tributary glaciers, many of which have their own names. Some of these tributary glaciers have become detached and now no longer directly feed into Tasman, resulting in the high number of glacier splits. The second highest number of split glaciers is Whymper Glacier with 13 single glacier units, and next are Murchison and Hooker Glaciers with 12 single glacier units, respectively.

The area of the largest glacier (Tasman) reduced from 101.6 to 82.8 km² (−19%; Fig. 12). Fox and Franz Josef, both west of the Main Divide and without proglacial lakes, reduced only slightly with −3 and −6% change, respectively. Murchison Glacier reduced from 33.9 to 24.4 km² (−28%). It is, as Tasman, located east of the Main Divide, has a large debris-covered tongue (total glacier nearly 30% debris-covered) and a large proglacial lake. In this (unrepresentative) example, no general difference was visible between glaciers east and west of the Main Divide and with or without proglacial lakes concerning area change. Fox and Franz Josef only had a small amount of debris cover (~1%) and showed the lowest area reduction. All other glaciers have debris-covered areas of >10% up to more than 40%. The mean area reduction for these 14 glaciers was 21% (−0.55% a⁻¹), and the reduction without Fox and Franz Josef was 27% (−0.71% a⁻¹).

Fig. 12. Glacier area change of selected debris-covered glaciers (see Fig. 11). Green asterisks = glaciers without proglacial lakes, blue asterisks = glaciers with proglacial lakes. Note log scale on x-axis.

5. Discussion

5.1 Uncertainties

Many challenges can occur in compiling glacier inventories using remote sensing. The main requirement for their compilation is the availability of suitable satellite images. This can be a problem in the maritime climate in NZ, which causes frequent cloud cover in some mountainous regions. Using Landsat 8 and Sentinel-2 satellites, this problem mostly could be resolved using just one satellite image. In some rare cases (max. 1% of the glacier area), the outline of a glacier was not totally visible. In these cases, we connected the existing outlines by more or less straight lines or followed the curvature of the visible glacier borders by checking the glacier outlines visible in the additional sources (see Materials and method). As we mapped semi-automatically on the Landsat 8 images and did the manual post processing on the Sentinel-2 MSI images, there was no conflict due to the different spatial resolutions of both sources.

One of the largest limitations involved in this research is found in the debris cover outlining process. We tested a number of other classification schemes such as Minimum Distance, Mahalanobis Distance, Parallelepiped and Spectral Angle mapper classification. We found that the maximum likelihood classifier applied to the PCA bands yielded the overall most robust result as a basis for manual correction supported by visual inspection using Sentinel-2 imagery. Because there were a large number of pixels highlighted as ‘debris’ in the maximum likelihood classification that were outside of the glacier extent (e.g. river gravels, miscellaneous debris throughout the landscape, see Fig. 3), the glacier outlines from 1978 were used as a mask. This introduced an uncertain assumption that the debris-covered extent would not have grown, except in an up-glacier direction, since 1978. As we completed the mapping by manual post processing, this potential error was greatly reduced or even eliminated.

Paul and others (2013) set the mean error of mapping glacier outlines to 3% using high- and medium-resolution images, with higher values for smaller and lower for larger glaciers. The required measurement uncertainty for glacier area from GCOS (2016) is 5%. Based on our uncertainty assessment, the evaluated mapping error was 4.3%. This value is in the range of the GCOS required uncertainty but is higher than the result of Paul and others (2013). Nevertheless, in their study they also found a large variation in mapping results. We therefore consider our result to be reliable and sufficiently accurate.

In addition, some mapping errors may be included in the new glacier inventory due to remaining perennial snowfields. The mass-balance year 2015/16, derived from EOSS measurements, was a negative year for the NZ glaciers (Willsman and others, 2017) but there were some late spring snowfalls leading to late snow, even if the snow line on the glaciers was quite high at the end of April. To avoid these mapping errors, we used images from February and April 2016 for Sentinel-2. The February images had the most snow-free areas in higher elevations (no new snow) and the April images the most snow-free areas in the lower elevations (the new snow has melted in lower levels but remained in the high). By applying these time steps, we have minimized the mapping error due to seasonal snow. Besides perennial snowfields connected to glaciers, there are perennial snowfields without any glacier contact. To avoid mis-mapping of these areas, we applied an area threshold of 0.01 km², assuming that these snowfields have smaller sizes. Nonetheless, it is possible that we mapped some of the larger snowfields erroneously.

Finally, elevations are based on a DEM (Columbus and others, 2011) which, due to the rapid retreat of some glaciers, no longer represents the surface elevation of all glacier areas in 2016. Hence, some of the lower elevation results may be overestimated.

5.2 Comparison with earlier glacier outlines

During the comparison of our inventory with the inventory from 1978 and the partial update from 2009, we found several errors in these older datasets.

The glacier outlines from 2009 were mainly based on ASTER satellite images. These images have a coarser spatial and especially radiometric resolution compared to the Sentinel-2 images (Table 1). Therefore, they did not provide any information in shadow, as revealed by the snow in shadow visible in the example shown in Figure 13. The green lines in this figure indicate the glacier outline from 2009 and run exactly at the border of the cast shadow (orange circle). ASTER passes NZ ~1 h later than Sentinel-2 MSI (11:30 NZDT), therefore the shadows are mainly consistent between 2009 and 2016 for late summer imagery in Sentinel-2 MSI/Landsat 8 OLI /ASTER. The shadowed area appeared probably black in the ASTER image used but glacier ice becomes visible in the Sentinel-2 image from 2016. This phenomenon happened often. Glacier recession and re-advance are ruled out due to the high elevation (upper part of accumulation areas) where generally little change occurred and that no advance/general growth happened between 2009 and 2016, based on knowledge from annual field-work and EOSS flights. Hence, the glacier outlines from 2009 excluded many glacier areas lying in cast shadow and identified these areas as glacier reduction. Using the Sentinel-2 images, we could see that this is incorrect and there is still glacier ice in these areas. We therefore did not compare the new inventory with the partial update from 2009.

Fig. 13. Mapping error in the 2009 glacier area update. The orange arrows point to glacier area in cast shadow (Glacier area 1978 and 2009: GLIMS, background image: Sentinel-2 from 12 February 2016).

The uncertainties included in the 1978 inventory were more complex and not as straightforward to distinguish. The mapping of the inventory took several years and was thoroughly performed using the support of oblique aerial images. This mapping was done manually on paper sheets. Afterwards, these paper data were digitized, but unfortunately this was done by a different person who was less familiar with glaciers. This means that while the overall calculation of areal ice coverage would be close to reality, random errors for specific glaciers were included in the 1978 inventory that are now inherent in the current GLIMS and RGI databases. Therefore, we did not compare our inventory with the 1978 inventory but only chose selected glaciers from the original 1978 paper sheet data to make a non-representative comparison of several large glaciers in the central glacier area. Before a detailed analysis of 1978–2016 changes takes place, the paper inventory needs to be re-digitized. If we only take the total area of 1158 km² from Chinn (2001), the reduction in glacier area is 31 or 0.8% per year.

5.3. Characteristics and change of New Zealand glaciers

Area change of a sample of 14 glaciers in the central Southern Alps showed a reduction of 21% (0.55% a⁻¹), since 1978. Results of this study were compared with two other glacier inventories from the Southern Hemisphere with slightly different inventory dates. At a lower latitude than NZ, Drenkhan and others (2018) mapped the glacier area change in the tropical Peruvian Andes (glacier area is ~1/6 of NZ) 1988–2016, for 2016 also using Sentinel-2 MSI. The area change (−37%) and change rate (−1.33% a⁻¹) are higher compared to the New Zealand sample.

On a similar latitude as NZ, Meier and others (2018) examined the glacier area in Patagonia (~30 times the NZ glacier area). The pattern of the distribution of number and area of glaciers is comparable and also the distribution of glacier area vs elevation. Additionally, they also found a positive relationship between maximum elevation and glacier area. The Andean main ridge has a clearer effect on the glacier mean elevation than the Main Divide in NZ. The area change (−8.8 ± 5.2%) and change rate (−0.29 ± 0.17% a⁻¹) 1986–2016 are much lower compared to our sample data. In an also quite maritime but small glacier area (~96 km²) in northern Norway on the Northern Hemisphere, the relative distribution of glacier number and area in 2014 showed a different pattern compared to NZ, with a lower number of smaller glaciers and a smaller area covered by the largest ones (Stokes and others, 2018). The area reduced 1988–2014 by 9% with a change rate of 0.32% a⁻¹; much less compared to NZ. The glaciers in this region also advanced in the 1980s until the beginning of 2000 and showed nearly no debris cover.

Approximately 10% of the glacier area in NZ is debris-covered. This proportion of debris cover is larger than the global average of 4.4% (Scherler and others, 2018), similar to the Karakoram and Pamir area (10%; Mölg and others, 2018), and less than the Greater Caucasus (13%; Tielidze and others, 2020). While Scherler and others (2018) find that the NZ region has the second-highest proportion of debris-covered area globally (after Caucasus and the Middle East) this result should be treated with caution given that the 10% coverage found here is significantly less than the 18% coverage of debris found by Scherler and others (2018). Another study (Herreid and Pellicciotti, 2020) found a similar debris value of 20%. A minor part of this difference resulted from the use of Landsat images from different time periods and, hence, potentially snow-covered debris areas. The main part of this discrepancy, related to both other studies, was attributed to the use of old glacier outlines in RGI (15% false-positive glacier area (Herreid and Pellicciotti, 2020)), which this study resolves by providing updated glacier outlines.

Almost all (89%) ice <1000 m a.s.l. is debris-covered. Some of these low elevation debris-covered glacier trunks are completely (e.g. Douglas Glacier, Fig. 11) or almost completely (e.g. Balfour Glacier, Fig. 11) separated from their accumulation areas by cliff faces. These regenerated glaciers are fed by ice avalanching and are active glaciers, 3.1 km (Douglas) and 6.7 km (Balfour) long.

Mean elevation of glaciers is higher to the west of the Main Divide compared to the east of the Main Divide. This result holds for overall average values, but is also true for all aspects except NW (Fig. 8). For glaciers with similar hypsometry, there is a direct relationship between mean glacier elevation and ELA (Furbish and Andrews, 1984), and hence this result seems to contradict another work that shows a clear increase in ELA from west to east across the main divide (Chinn and Whitehouse, 1980). The most likely reason for that is a systematic difference in morphology of glaciers east and west of the Main Divide. For example, if we consider the 10 largest glaciers in the Southern Alps the mean elevation of the four glaciers east of the Main Divide (Tasman, Murchison, Hooker and Lyell Glaciers) is low because they have extensive debris-covered tongues in the valley floor. Conversely, the six largest glaciers west of the main divide (Fox, Franz Josef, Bonar, Olivine Ice Plateau, Lambert and Volta Glaciers) are all largely clean glaciers, which terminate at relatively high elevations, or have steep, narrow tongues. While these are only 10 glaciers, collectively they cover 233 km² (29% of all ice area).

6. Conclusions and summary

In this study, we mapped the glacier area of New Zealand in 2016 first semi-automatically by means of optical remote sensing (Landsat 8 OLI). Manual post processing, mainly by using Sentinel-2 MSI images, was necessary due to misclassifications (e.g. lakes, clouds), mapping in shadowed areas, and the connection of bare and debris-covered sections, mostly based on Sentinel-2 MSI images. The glacier area of the 2918 mapped glacier entities is 794 km² with a debris cover of 10%. Of these glaciers, 71% are very small (<0.1 km²) and only seven glaciers >10 km². Tasman Glacier as a glacier complex with 26 tributary glaciers is the largest glacier with 82.8 km². Fifteen glaciers are located on the North Island with a total area of 3.0 km² and a debris cover of 47%. Minimum elevation of the glacier area is 420 m a.s.l. and maximum elevation 3724 m a.s.l. with a mean elevation of 1953 m a.s.l. Glaciers west of the Main Divide have mostly a higher mean elevation compared to glaciers in the east. There is no correlation between the amount of debris cover and the mean aspect in the ablation area.

Uncertainties of this inventory are mainly based on snow patches erroneously mapped as glacier area, the interpretation of debris-covered glaciers in shadow and/or potentially covering already disconnected dead ice. The value of this uncertainty is difficult to quantify, but based on our uncertainty mapping, we found a value of 4.3% resulting in a total glacier area of 794 ± 34 km².