3.1. Data

A complete and unique inventory of New Zealand glaciers was carried out by T. Chinn using aerial photographs from 1978 (Chinn, unpublished). Glacier outlines were photo-interpreted from the oblique aerial photographs recorded in 1978 and were digitized onto topographic base maps at a scale of 1 : 63 360 (imperial system where 1 in = 1 mile). The contour interval of these maps was 30.5 m (100 ft). In addition to mapping all the glaciers, 48 index glaciers, well distributed along the Southern Alps, were chosen for annual snowline survey. Oblique aerial photographs of the index glaciers have been taken every year since 1978, except in 1979, 1990 and 1991.

For updating the New Zealand glacier inventory in the main glacierized area, the central Southern Alps, an image from the ASTER sensor was used. ASTER captures high spatial resolution data in 14 bands, operating in three subsystems: visible and near-infrared (VNIR) consisting of three spectral bands with a spatial resolution of 15 m; SWIR consisting of six spectral bands with a resolution of 30 m; and thermal infrared (TIR) consisting of five spectral bands (Table 1). The best available image was selected considering the following criteria: (1) cloud-free, (2) end of summer, and (3) year with a negative mass balance. The image chosen was recorded on 14 February 2002 and was a processed level-1B image (Fig. 1, image ID: SC:AST_L1B.003:2019534558). The year 2002 was a very suitable year for compiling the glacier inventory since it had the most negative mass balance between 2001 and 2005 (Reference Chinn, Heydenrych and SalingerChinn and others, 2005). This meant that most of the surrounding snowpatches remaining in late summer in other years were gone and could not be misinterpreted as glaciers (Reference ShangguanShangguan and others, 2006). The image contains ∼40% of the total glacier area mapped in the 1978 glacier inventory. Six of the forty-eight index glaciers used for the annual snowline monitoring are also situated within the image (Fig. 1). Landsat images might have permitted us to investigate a longer time series, but were not considered here due to their lower spatial resolution (30 m). ASTER was believed to have the appropriate spatial resolution to map New Zealand glaciers. Changes in glacier area of the small, alpine glaciers, which contribute by far the largest numbers of individual glaciers in New Zealand, would have been more difficult to detect.

Table 1. ASTER instrument characteristics (adapted and modified from Reference Gao and LiuGao and Liu, 2001)

The image was orthorectified (PCI Geomatica OrthoEngine software) to ensure an appropriate mapping accuracy. A total of 31 survey points were collected in the field from late February until late April 2005 using a Trimble Pro-XR GPS instrument. These features were relatively easily detectable on the satellite image and were chosen as ground control points (GCPs). These were, for instance, centers of bridges, sharp intersections between scree slopes and vegetation, or characteristic rock features in the fringe of glaciers. Only 22 of these GCPs were later found to be of high enough quality to be used for further processing, i.e. easy to separate from the surrounding pixels in the image and easy to recognize in the field. The GPS data were post-processed using the carrier phase technique (sub-metre theoretical accuracy) with reference data collected from the Trimble 40 GPS Base Station located at the University of Otago. The field-collected survey points were divided into GCPs and independent check points (CPs). The first orthorectification performed with a digital elevation model (DEM) from the Land Information New Zealand (LINZ) topographic database, with a pixel size of 15 m, showed large shifts in pixel position when overlaying geographical features from LINZ. This was likely due to the planimetric and altimetric extrapolation of the geometric sensor model resulting from the spatially clustered and uneven three-dimensional (3-D) distribution of the field-collected control points (Reference ToutinToutin, 2004). To improve the sensor model and orthorectification result, height points from the New Zealand topographic database were added. Finally 28 GCPs and 15 CPs were used for the orthorectification. The final orthorectification error (root-mean-square (RMS)) for the GCPs (X and Y) was 15.3 and 14.5 m, respectively, and for the CPs (X and Y) was 17.1 and 14.7 m, respectively.

3.2. Extracting land ice area using classification methods

The automatic extraction of debris-free glaciers from multi-spectral satellite imagery relies on the spectral uniqueness of ice and snow compared with surrounding targets (soil, vegetation, rock): high reflectance in the visible (0.4–0.7 μm), and medium and low reflectance in the NIR and SWIR (0.7–2.5 μm), respectively. Glacier boundaries have been mapped successfully in various places in the world using different techniques: (1) photo interpretation and manual delineation on colour composites (Reference Williams, Hall and ChienWilliams and others, 1997; Reference Khromova, Osipova, Tsvetkov, Dyurgerov and BarryKhromova and others, 2006; Reference BeedleBeedle and others, 2008); (2) image segmentation (or thresholding) of principal components (Reference Sidjak and WheateSidjak and Wheate, 1999), single bands (e.g. Thematic Mapper 3 (TM3); Reference Bronge and BrongeBronge and Bronge, 1999), single band ratios (e.g. ASTER 3/4; Reference KääbKääb and others, 2003) or spectral indices (e.g. NDSI; Reference Hall, Bayr, Bindschadler, Shöner, Hardy and FrankensteinHall and others, 2001); and (3) unsupervised and supervised classifications (e.g. ISODATA, maximum likelihood algorithm; Reference Aniya, Sato, Naruse, Skvarca and CasassaAniya and others, 1996). Thresholding of band ratios or spectral indices is generally the preferred approach because of its simplicity, low requirements in terms of time and processing power, and ease of implementation in an operational context. Only one study reports on the extraction of land ice in New Zealand and it used the NDSI approach (Reference Mathieu, Chinn and FitzharrisMathieu and others, 2009). However, the delineation of ice/non-ice targets was not assessed in great detail as the main aim was to assess the suitability of ASTER data to extract ELAs. Advanced image-processing techniques such as spectral unmixing have been used efficiently for retrieving subpixel fraction cover either for snow or ice (Reference Nolin, Dozier and MertesNolin and others, 1993; Reference Klein and IsacksKlein and Isacks, 1999). However, a per pixel approach was considered more appropriate in this work because subpixel fraction cover derived from a 15 m spatial resolution ASTER image would not provide any additional accuracy benefits when comparing with the 1978 national glacier inventory. Indeed, the glacier boundaries in the earlier inventory are believed to be at best ±30 m accurate, corresponding to the accuracy of the topographic base maps (1 : 63 000) used to report the photo-interpreted features.

In this study, we assess the suitability of three common classification techniques to extract land ice in the New Zealand context: (1) band ratio (ASTER 3/4), (2) NDSI (ASTER 1 − ASTER 4)/(ASTER 1 + ASTER 4), and (3) supervised classification (maximum likelihood algorithm). The ratio of TM4/5 (ASTER 3/4 equivalent) was found to be the most appropriate method for the new Swiss glacier inventory (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002). As we were only working with one image, no atmospheric correction was performed. Further, ratios were performed directly on raw digital numbers as this was found to be the most robust approach in an operational context (e.g. Reference AlbertAlbert, 2002; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002, Reference Paul, Huggel and Kääb2004a).

In the case of the ASTER 3/4 ratio, thresholds of 1.8, 2.0, 2.2 and 2.4 were tested to separate glacier and non-glacier. For the NDSI, thresholds of 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65 and 0.7 were tested. Nearly all lakes in the area were classified as glaciers, independent of their turbidity in the NDSI-based glacier products. Turbid water, like the Tasman Lake, was still classified as glacier even with the highest threshold tested. However, this was not a significant problem for the ASTER 3/4 ratio. All lakes were clearly visible in the colour composite and were therefore manually digitized and removed from the NDSI results. The supervised classification was only tested for one small test area of the image. Supervised classification demands substantial work for selecting the best possible training sites.

Significant misclassifications on the glaciers were expected because of the presence of heavily debris-covered glacier tongues in many of the largest glaciers in the area, in addition to the extreme vertical relief causing large shadows. Consequently, it was decided to manually digitize the glacier areas of the entire image. This option was considered the most reliable approach to quantify glacier changes between 1978 and 2002. The manual delineation was done by photo-interpreting the image displayed as colour composite using the VNIR channels (NIR, red and green). It is important to note that individual glacier boundaries were not mapped in this study, but only the land ice area as a whole. Glacier area changes were calculated from glacier area polygons as measured in 1978 and 2002, using simple change overlay techniques. In addition, more detailed maps were produced to analyse how the changes were distributed in space.

3.3. Validation of the classification methods

The classification methods were validated by (1) qualitatively comparing GPS outlines for selected glaciers with the manual and automated classifications and (2) quantitatively comparing the manual classification (i.e. digitized product) with the automated classifications (indices and supervised classification). GPS outlines of selected glaciers were surveyed in the field for the validation of the classification methods. To reduce the potential negative effects of the time difference between the fieldwork (2005) and the image collection (2002), index glaciers were preferred since aerial photographs of these glaciers existed for each year between 2002 and 2005. The index glaciers are monitored annually during the snowline survey, and possible changes could therefore be carefully considered. The surveyed glaciers are well distributed over the study area, covering both eastern and western glaciers (Fig. 1). Four of the six index glaciers located in the area were of convenient size and accessibility for GPS ground mapping: Jalf/Bauman, Chancellor Dome, Langdale, and Ridge glaciers (Fig. 2). In addition, Tasman (which is also an index glacier), Hooker, Fox, Franz Josef and Stocking glaciers were also investigated in the field (Fig. 3). GPS outlines of the individual glaciers were obtained by walking along the boundary between snow and rock. The GPS receiver was set to record a new position every 5 s, calculated from at least five satellites. This criterion was decided in order to ensure a high positional accuracy of the GPS data. The glacier outlines were post-processed using the carrier phase technique with reference data collected from the Trimble 40 GPS Base Station. Both the carrier phase processing technique and the constant lock of the GPS receiver with at least four satellites ensured systematic sub-metre accuracy, even with a limited number of position records for a single point (i.e. dynamic GPS data collection). Higher accuracy is achieved by using the higher frequency of the L1/L2 carrier signal as opposed to the pseudo-random code. Hence, the accuracy of the GPS glacier outlines was much more dependent on the field interpretation rather than on inaccuracies in the GPS recordings. For some glaciers and because of safety issues, the GPS positions could not be recorded at the exact boundaries (e.g. frontal part of Hooker and Tasman Glaciers (Fig. 3)). In these cases, the average distance to the exact boundaries was noted and the positions were later adjusted in Geographical Information System (GIS) software. Once the manually digitized map was assessed with the field-surveyed GPS outlines, this map was used as an efficient way to select the best thresholds for the index-based classifications and to validate all automated classifications.

Fig. 2. Comparison of the 2005 GPS-surveyed glacier outlines with the manual digitizing based on the 2002 ASTER image for the five smallest glaciers studied in the field.

Fig. 3. Comparison of the 2005 GPS-surveyed glacier outlines with the manual digitizing based on the 2002 ASTER image for Tasman (left) and Hooker Glacier tongues (right). The part with no data in the GPS line (indicated by the black arrow) is where this pond was located at the time of the fieldwork (i.e. no longer a pond, but part of the large lake).

4.1. Performance of manual interpretation

In this section we first present the assessment of the manual interpretation because it served as a validation set for the automated classification methods. The manually digitized glacier outlines were compared with the GPS-collected glacier outlines in order to assess the accuracy of the interpretation. For the smaller glaciers, there is a good correlation between the manual digitizing and the field outlines (except from Stocking Glacier, which has been retreating between 2002 and 2005). This indicates that the manual digitizing was fairly accurate for these small glaciers, and that they have changed little during the 3 years (Fig. 2). The validation data as well as the operator experience suggest that the manual digitizing is accurate within 2 pixels (30 m) (±15 m). This would ignore the fact that glacier/nonglacier boundaries are not hard but tend to be fuzzy in nature, and that small discrepancies might still be attributable to time difference between the imagery and the fieldwork. Although difficult to evaluate, this error is likely to increase in the case of steeper glaciers, in shadow where the snow/ice and rock contrast will not be as good, and in the case of debris-covered areas. The termini of the larger glaciers had changed too much in the same period to be used for validation purposes. Nevertheless, the lateral margins of these glaciers also seemed to match relatively well (within 2–3 pixels), despite the difficulty in separating rock and debris-covered ice (Fig. 3). Considering the previous observations the overall accuracy of the digitizing of the 2002 glacier outlines is estimated to be within 3 pixels (±22.5 m).

4.2. Performance of automated classification techniques

The automated classifications were compared visually with the GPS outlines and the manually digitized outlines. Then confusion matrices were computed between the manually digitized glacier maps (acting as a reference dataset) and the automated classifications (Tables 2, 4 and 6). Producer’s, user’s and overall accuracy were calculated from the confusion matrices (Tables 3, 5 and 7). The overall accuracy is the total proportion of pixels classified correctly. User’s and producer’s accuracy are equivalent to type I and type II errors (for details see Skidmore, 1999). For simplicity, accuracy measures are only presented for the threshold values selected as optimum (2.0 for ASTER 3/4 ratio and 0.5 for NDSI).

Table 2. Confusion matrix results for the band ratio image ASTER 3/4. Values are given in % relative to the total number of pixels within the image (N = 1.7 × 10⁷ pixels)

Table 3. Producer’s, user’s and overall accuracy for the two classes (glacier and non-glacier areas) for the band ratio image ASTER 3/4, with and without 3 × 3 median filter. For these results the threshold is set at 2.0

Table 4. Confusion matrix results for the NDSI image ((ASTER 1 − ASTER 4)/(ASTER 1 + ASTER 4)). Values are given in % relative to the total number of pixels within the image (N = 1.7 × 10⁷ pixels)

Table 5. Producer’s, user’s and overall accuracy for the two classes (glacier and non-glacier areas) for the NDSI image. For these results the threshold is set at 0.5

Table 6. Confusion matrix results for the supervised classification using the maximum likelihood algorithm. Values are given in % relative to the total number of pixels within the image subset of Franz Josef Glacier area (Fig. 1; N = 2.8 × 10⁵ pixels)

Table 7. Producer’s, user’s and overall accuracy for the two classes (glacier and non-glacier areas) derived from supervised classification using the maximum likelihood algorithm

In the confusion matrices the ideal threshold is where both the omission error (non-glacier detection) and commission error (false glacier inclusion) are as low as possible. However, if choosing the lowest value of one of the errors, the other will be very high. Therefore a choice has to be made when choosing thresholds based on what one wants as an outcome. If glacier mapping over a large and heterogeneous area is the main target (as for this case), the ideal threshold should not exclude too much (what is glacier) nor include too much (what is not glacier). If only the small and very steep glaciers were of interest it could be wise to choose a threshold that included more of the glaciers in shaded areas (and a higher commission error would have to be accepted). In the New Zealand context, with a high degree of debris-covered ice on many long valley glaciers, manual correction was required, as most of these areas were left out of the classifications. Thus, for the ASTER 3/4 ratio and the NDSI, the task was to find the threshold where a minimum manual adjusting of land ice boundaries was needed. By inspection of the ASTER 3/4 ratio image, the best threshold was found to be ∼2.0 (Fig. 4), this result being also confirmed with the confusion matrix (Table 2). However, major changes in glacier areas will be detected even though the threshold used is not optimal.

Fig. 4. Comparison of the various glacier extraction techniques for the Ridge Glacier area based on the 2002 ASTER image: (1) manual digitizing, (2) band ratio (ASTER 3/4, threshold 2.0) and NDSI (threshold 0.5). Bottom left panel is a zoom-in of Ridge Glacier.

In Tables 2 and 4 (ASTER 3/4 and NDSI results for the entire image) we note that >20% of the glacier area defined in the photo interpretation is not detected in the classification (Omission error). Most of this error by far can be attributed to debris-cover glacierized areas of long valley glaciers and others which are not detected as glaciers. We also note that applying a 3 × 3 median filter to this classification slightly increases the accuracy (Table 3). This filtering was efficient for removing noise in areas with dense vegetation, but outside these areas it did not influence the classification results significantly. Areas defined as non-glacier in the photo interpretation, but classified as glaciers in both the ASTER 3/4 and the NDSI maps (Commission error), were typically very bright rocks or relatively small snow/icy patches, often located on terrain features where it was not likely to be a glacier (e.g. small hills separated from the surrounding topography). The initial test of supervised classification did not give very promising results. We compared the results obtained for the supervised classification over a limited area around Franz Josef Glacier (Tables 6 and 7; Fig. 5) with the results obtained over the same area for the best thresholds applied to the ASTER 3/4 ratio and the NDSI (results not shown). The supervised classification both excluded more area classified as glacier in the digitizing and included more area not classified as glacier in the other two methods. Table 6 shows this clearly, where both the omission and commission error are very high. The overall accuracy given in Table 7 is also very low. The strips clearly visible in Figure 5a are due to residual errors in the radiometric correction of the SWIR band, and are in the order of 1 digital number.

Fig. 5. Supervised classification (using maximum likelihood algorithm) of Franz Josef Glacier based on the 2002 ASTER image: (a) classification into seven classes; (b) classification merged into two classes (glacier and no glacier); (c) comparison of the classification with the manual digitizing.

4.3. Glacier area change from 1978 to 2002

In the 2002 ASTER image the area of the manually digitized glacier outlines totalled 428 km². This represented almost 11% of the entire image. Of this 428 km², 242 km² were located on the eastern side of the Main Divide and 186 km² on the western side of the Main Divide (Table 8). The glacier polygons of the 1978 inventory were clipped to match the ASTER image. Within the limits of the ASTER image, a total of 513 km² of glacier area was calculated, which accounted for ∼41% of the glacierized area in New Zealand in 1978. Of this total, 296 km² were located east of the Main Divide and 217 km² west of the Main Divide (Table 8).

Table 8. Changes of glacier areas in the central Southern Alps between 1978 and 2002

Changes during the 24 year period in the central Southern Alps indicate that the total glacier area decreased from 513 km² to 428 km² (a loss of 16.6%). We estimated the uncertainty for the calculated areas to be ±1.1 % (or ±5 km²) for the 1978 glacier outlines and ±0.6% (or ±2 km²) for the 2002 glacier outlines. The uncertainty in the total glacier area change was ±6.8% (or ±6 km²) (see section 5.2 for details of uncertainty estimations). The glacier retreat was smaller on the western side of the Main Divide (14.3%) than on the eastern side (18.3%). The two main reasons for this are that (1) the two largest glaciers on the western side (Fox and Franz Josef Glaciers) advanced during the study period, and (2) the development or enlargement of proglacial lakes on the large eastern glacier rapidly shrank the glaciers by lake calving. On the western side, Fox and Franz Josef Glaciers advanced 300–650 m in the time period studied, while La Perouse and Victoria Glaciers both retreated 1900 m (Fig. 6). On the eastern side all the largest glaciers retreated 1–2.1 km, most significantly Tasman Glacier (Tasman Glacier 2100 m, Murchison Glacier 1000 m, Hooker and Classen Glaciers 1700 m; Fig. 7). The enormous loss of ice mass due to the recent growth of these glaciers’ proglacial lakes contributes to such large area losses on the eastern side that the large retreat on most of the other valley glaciers on the western side of the divide (e.g. Victoria and La Perouse Glaciers (Fig. 6)) is small by comparison. If the eastern glaciers with proglacial lakes were discounted, the difference in relative decrease of glacier areas between the eastern and the western side of the Main Divide would be much smaller or even might be reversed.

Fig. 6. Area changes measured for the glaciers west of the Main Divide in the central Southern Alps from 1978 to 2002. Fox and Franz Josef Glaciers advanced several hundred meters whereas the debris-covered glaciers (e.g. La Perouse and Victoria Glaciers) retreated dramatically.

Fig. 7. Area changes measured for the glaciers east of the Main Divide between 1978 and 2002, showing the enormous retreat in glacier area due to the rapid expansion of proglacial lakes with glacier calving.

In contrast, the smaller, high-elevation glaciers showed only minimal retreats compared with the large valley glaciers. Figure 8 illustrates the 1978 and 2002 glacier outlines for the five smallest glaciers studied in the field. It is not easy to detect significant glacier area changes. Some of these glaciers show small retreats, some have approximately the same size, and some even seem to have increased in some parts between 1978 and 2002. The latter observation is consistent with a trend to positive mass balances over the past two decades (Reference Chinn, Heydenrych and SalingerChinn and others, 2005).

Fig. 8. Comparison of the 1978 and 2002 glacier outlines for the smallest glaciers studied in the field.

5.1. Mapping accuracy of automated classification techniques

Mapping debris-covered glacier ice in a multispectral image is a recurring challenge because the spectral signature of the debris often prevails over the ice or snow signature if a certain percentage of the glacier pixels are debris-covered. Often confused with surrounding rocks, debris-covered glacier tongues are commonly manually added to classifications as a post-classification step (e.g. Reference Williams, Hall and ChienWilliams and others, 1997). Reference Paul, Huggel and KääbPaul and others (2004a) describe a semi-automatic method for delineating debris-covered glaciers which combines the advantages of automated multispectral classification for extracting clean glacier ice with slope information derived from a DEM. One important question raised by our research is how debris-covered glacier tongues can be mapped accurately from satellite imagery, even by using manual delineation, when it is nearly impossible to recognize their exact boundaries in the field. However, even in the worst cases misinterpretation in the field will not deviate more than a few pixels from the true boundary, as long as we assume that the ice beneath the debris is part of the glacier.

Snow in shadow tends to have reflectance properties similar to those of bright rock, and is another challenge for the automatic extraction of glaciers. Manual detection of these surfaces was relatively straightforward as long as the interpretation was careful, and topographic maps and other terrain information (e.g. DEMs and aerial photographs) were used. However, these areas were not always easy to extract with automatic methods, depending on the type of snow or ice (clean snow versus dirty ice) and steepness of the slope (snow in steeper slopes tended to be harder to extract). Most of the shaded clean snow could be accurately delineated using the lowest threshold tested for the ASTER 3/4 ratio. However, this threshold would also include large areas which obviously were not glaciers, so a trade-off was required to balance omission and commission errors. Many dirty snow/ice patches in shadow were left out. The NDSI was also shown to detect shaded clean snow relatively well, but not as well as the ASTER 3/4 ratio.

When we applied the NDSI, nearly all lakes in the area were classified as glaciers (Table 4). Turbid water, like the Tasman Lake, was still classified as glacier ice even with the highest threshold tested. The reason for this is the similar reflectance of very turbid lakes (with high load of rock flour in the water) and glacier ice in the green spectral domain. Increasing the NDSI threshold increased the omission error significantly (Table 4). However, this was not found to be a problem with the ASTER 3/4 ratio. In New Zealand the lakes are clearly visible in the VNIR composite image, and the water bodies were therefore manually digitized and removed from the NDSI results. After the water mask was applied to the NDSI results, this method revealed results similar to the ASTER 3/4 ratio (Fig. 4). In some areas the NDSI could be even more efficient. Similar problems of misclassifications of lakes as glaciers were reported by Reference Paul and KääbPaul and Kääb (2005) in Arctic Canada when using the Landsat Enhanced Thematic Mapper (ETM) 3/4 method. The best threshold for the NDSI in this study was ∼0.5 (Table 4), in accordance with other studies (e.g. Reference Mathieu, Chinn and FitzharrisMathieu and others, 2009). Band ratios have been reported to be an especially simple, robust and fast method for glacier mapping over large areas (e.g. in Switzerland (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002) and Peru (Reference AlbertAlbert, 2002)), and were confirmed as good performers in this study in New Zealand. The ASTER 3/4 ratio with a threshold around 2.0 was the most efficient (Table 2). This is also similar to other studies (e.g. Reference Paul, Kääb and HaeberliPaul and others, 2007).

5.2. Accuracy of glacier change detection between 1978 and 2002

The accuracy of glacier change detection between 1978 and 2002 was constrained by the difficulty of photo-interpreting glacier boundaries in both datasets. Figure 8 shows a comparison of the 1978 glacier outlines drawn from oblique aerial photographs and the manual digitizing of glaciers from the 2002 ASTER image for a selected number of glaciers, also delineated in the field. Drawing glacier outlines from oblique photographs is time-consuming and the accuracy depends on the local knowledge of the interpreter as well as the availability of good cartographic information for the area, like topographic maps and DEMs. This is apparent in Figure 8 where more details are evident in the digitizing of the satellite image compared with the 1978 interpretation, even though the ASTER image has a lower spatial resolution than the oblique photographs. The 1978 mapping was done at a coarser scale, mostly because details were practically impossible to transfer to the paper topographic maps. The 1978 digitizing was done on topographic base maps at a scale of 1 : 63 360 (imperial system where 1 in = 1 mile), where 1 mm corresponds to ∼63 m, a distance that includes four ASTER pixels. For the 1978 inventory we believe that in most areas the accuracy of the digitizing onto the base maps is within 1 mm. However, some corners are slightly cut off and some of the outlines are slightly generalized (e.g. Fig. 8). We conservatively estimated the accuracy of the 1978 glacier boundaries to be within 2 mm (±1 mm), giving us a ground-equivalent distance of ±63 m.

The accuracy of the 2002 manual digitizing firstly depends on the quality of the orthorectification. In the final orthorectified product, we achieved a root-mean-square error (RMSE) of 17.1 and 14.7 m for the CPs, in the X and Y directions respectively. In practical terms, these results are slightly more than a pixel (15 m), a standard level of accuracy when orthorectifying satellite imagery. The consistency of these results was verified for the entire image by visual inspection of the orthorectification overlaid by rivers and other vector information from the New Zealand topographic database. The other source of errors is related to the digitizing of the glacier boundaries. We estimated the uncertainties of the 2002 digitized glacier outlines for the most part to be within two pixels (30 m) (i.e. ±15 m). No corner was cut and every pixel was accounted for (Fig. 2). In areas of doubt (e.g. steep areas with glaciers in shadows), aerial photographs, 3-D model and field observations were used to support our interpretation. However, in order to accommodate difficulties of interpreting debris-covered tongues and the problems of interpreting steep mountain faces (e.g. east face of Mount Cook; Figs 9 and 10), we estimated the total uncertainty in the manual digitizing to be within 3 pixels (±22.5 m). By combining the uncertainty in the manual digitizing with the uncertainty of the orthorectification, we obtained an overall uncertainty estimation for the 2002 glacier outlines of ±28.5 m.

Fig. 9. Example of comparison of the 1978 and 2002 glacier inventories: (a) area changes measured in the Mount Cook area and (b) corresponding ASTER orthoimage displayed as it was interpreted. 1 marks the high peak of Mount Cook (3754 m a.s.l.) and 2 marks the middle peak (3717 m a.s.l.). The arrow points to a spot on the east face of Mount Cook (shown in Fig. 10) apparently misinterpreted in the 1978 inventory (in fact rock avalanche debris), while some of the red areas could be misinterpreted in the 2002 inventory.

Fig. 10. Photograph of east face of Mount Cook taken from Langdale Glacier in April 2005, showing the face at the time of year with least snow and ice. 1 is the high peak (3754 m a.s.l) and 2 is the middle peak (3717 m a.s.l.) shown in Figure 9.

In order to calculate the resulting uncertainty in area, we took each single glacier polygon (for both inventories) and generalized them into circles (>600 digitized polygons within the study area in 2002, many containing several glaciers). This is a simplification, since the polygons certainly have more complex shapes; however, it should give us a relatively good error estimate. The uncertainties in area were calculated for each circle, both for the 1978 and the 2002 data. The uncertainty in the total glacier area change was ±6.8% (±6 km²) (Table 8).

Snow pixels show up as bright white in the visible imagery because of their high reflectance. However, whether this snow is part of a glacier or belongs to an occasional snowpatch is impossible to determine from the spectral signature itself. The interpretation is subjective, depending on what the interpreter sees as a glacier or not. In positive mass-balance years, with ample snow and a cold ablation season, many snowpatches will survive the summer and be extracted as glaciers from a satellite image by automated classification. If these areas were separated from the glacier, most could easily be removed by setting a minimum size for glacier polygons (e.g. 0.01 km²). However, if these snowpatches are connected to the glaciers, this will cause a larger area to be classified as glaciers compared with the true glacier area. To minimize these potential confusions, it is important to select an image from a negative mass-balance year.

In addition to the difficulty of distinguishing debris-covered glacier parts from rocks, it was also problematic to photo-interpret the upper extent of glaciers near steep rock faces. For the 1978 glacier inventory, the definition was that if these steep faces contributed ice by flow, they were part of the glacier, but if they were disconnected and only contributed by winter snow avalanches, the faces were not part of the glacier (Chinn, unpublished). However, this distinction was almost impossible to make only from the ASTER image because of the coarse resolution of the image relative to the features. Steep rock faces were not excluded from the total glacier area unless it clearly seemed to be bare rocks. The comparison of the 1978 and 2002 inventory maps shows many red areas around the highest peaks along the Main Divide (see example around Mount Cook in Fig. 9a), suggesting that glaciers would have grown in these areas during the study period, something which is not likely. Most of these red areas were classified as glaciers in 2002 because they were snow-covered at the time of image acquisition, whereas in reality they were just snow-covered rocks (cf. Fig. 9). The arrow in Figure 9 points to an area that was most likely wrongly classified as a glacier in the 1978 inventory. After closer investigation we found that this area is in fact a corridor of active rock avalanche. This is another illustration of the complexity of tracing glacier outlines from aerial photographs or satellite images. However, considering the two inventories as a whole, circumstances like active rock avalanches (leading to classifying part of glaciers as non-glaciers) and snow-covered mountain faces (leading to classifying rocks as glaciers) are minor and hardly influence the total outcome of the two inventories. When estimating the uncertainties in the two inventories, we have accounted for all possible error sources, and our estimates were based on conservative assessments. These estimates point to small uncertainties in the glacier area calculations and we can argue that the overall pattern and magnitude of glacier retreat would not change significantly after a further refinement of these outlines.

5.3. Change of glacier area in New Zealand from 1978 to 2002

The results show a measured retreat of the glaciers in the central Southern Alps of ∼17% during the period 1978–2002. Similar results have been reported in other parts of the world. The new Swiss glacier inventory (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2004b) measured the relative loss in glacier size from 1973 to 1998/99 as ∼20%. However, in Switzerland, glaciers smaller than 1 km² (100 ha) contributed to ∼40% of the total loss, while they cover only 15% of the area. About 18% of this change occurred from 1985 to 1999, illustrating a recent acceleration of the glacier retreat. This strongly contrasts with the New Zealand situation where the largest glaciers have contributed to the largest ice loss. Reference Quincey and GlasserQuincey and Glasser (2009) reported a retreat of Tasman Glacier by 3.5 km between 1990 and 2007. They also found that the area of the proglacial lake doubled between 2000 and 2007 alone. In Norway, the maritime glaciers along the western coast with a large annual mass turnover had a large mass surplus between 1962 and 2000. Conversely, the continental glaciers with small summer and winter balances had a mass deficit over the same period (Reference Andreassen, Elverøy, Kjøllmoen, Engeset and HaakensenAndreassen and others, 2005). Since 2001, all monitored glaciers in Norway have shown a considerable mass deficit. Overall the Norwegian glaciers have retreated during the 20th century. Continental glaciers have, with few exceptions, retreated throughout the century, while many maritime glaciers have been through periods of advance and recession, although recession has been the main process (Reference Andreassen, Elverøy, Kjøllmoen, Engeset and HaakensenAndreassen and others, 2005). Reference Andreassen, Paul, Kääb and HausbergAndreassen and others (2008) reported glacier area shrinkage at investigated glaciers in Jotunheimen National Park, Norway, by 23% since the 1930s for 38 glaciers. The largest reduction occurred in the first part of the period, between 1930 and 1965. Since 1965 the area reduction has been 12% for 164 glaciers. Reference Paul and KääbPaul and Kääb (2005) calculated a decrease in area for 225 glaciers in Cumberland Peninsula, Baffin Island, Arctic Canada. They compared Landsat ETM+ and ASTER scenes from August 2000 with manual delineation of LIA trimlines and moraines. They found an average area loss of 11% over the period. The relative loss of area increased towards smaller glaciers (<10 km²) and was more or less constant for larger glaciers. Reference Khromova, Osipova, Tsvetkov, Dyurgerov and BarryKhromova and others (2006) assessed glacier recession in the east Pamirs, central Asia, using historical data and space images from Russian satellites and ASTER. They found a 7.8% decrease in glacier area during 1978–90, and 11.6% during 1990–2001. In addition to the decrease in glacier area and retreat of glacier fronts, they also found increased debris-covered area and the development of new lakes. All these reports of glacier retreat are similar to what is found in this study, but the reasons behind the recessions might be very different.

Qualitative analysis of our results suggests that the smallest glaciers in New Zealand do not contribute to the largest ice loss (Figs 6–8). The large heavily debris-covered valley glaciers have lost enormous areas over the past few decades, whereas some of the smaller glaciers may even have shown a small advance. This last observation was corroborated by Reference Allen, Owens and SirgueyAllen and others (2008) who reported that Stocking and Eugenie Glaciers (east of the Main Divide) advanced up to 100–140 m between 2002 and 2006.

The clustering of these large valley glaciers in the central Southern Alps covered by the 2002 ASTER image indicates that the area change for the study period is likely to be larger here than for the rest of the glacierized areas in New Zealand. This suggests that if the 2002 glacier inventory had been made for the whole of New Zealand, the total area changes from 1978 would probably have been lower. The response time for alpine glaciers is ∼30 years (Reference Cuffey and PatersonCuffey and Paterson, 2010), while response times for high-velocity glaciers like Franz, Fox and Stocking are in the order of 5–8 years (Reference Fitzharris, Lawson and OwensFitzharris and others, 1999). Some of the large valley glaciers maintained their LIA areas (1850–90) up until the 1970s, giving a response time of ∼100 years. The relatively thick ice of the lower reaches of these glaciers is well insulated beneath thick protective mantles of debris and became largely relict. Ice loss was shown by downwasting (surface lowering) alone, with little to no terminus retreat. These glaciers are in a state of extreme disequilibrium with the present climate. However, during the past few decades, many have reached irreversible tipping points that have triggered proglacial lake expansion into terminus areas (Reference Quincey and GlasserQuincey and Glasser, 2009). These lakes, with their rapid expansion by both melt and glacier calving, are destroying the lower trunks of the large valley glaciers, making the glacier behaviour not directly linked to the climate (Reference Hochstein, Watson, Malengrenau, Nobes and OwensHochstein and others, 1998; Reference Kirkbride and WarrenKirkbride and Warren, 1999; Reference Purdie and FitzharrisPurdie and Fitzharris, 1999; Reference RöhlRöhl, 2006; Reference Winkler, Chinn, Gärtner-Roer, Nussbaumer, Zemp and ZumbühlWinkler and others, 2010). Reference ChinnChinn (1999) reported a reversal of the past century glacier recession trend. Recent interannual fluctuations of ELAs in New Zealand indicate that alternate positive and negative mass-balance years (1998–2008) have followed two decades of positive mass balance (1977–98) (Reference Willsman, Salinger and ChinnWillsman and others, 2008). Figures 6 and 7 show the changes in glacier area for the western and eastern glaciers in the 24 year period. Only the two largest of the west coast glaciers have advanced, whereas all the other large glaciers have retreated. It may be inferred that the western climate is more favourable for glacier advance, compared with east of the Main Divide. However, Reference ChinnChinn (2001) argues that the notion that the ‘western’ glaciers behave differently from the ‘eastern’ glaciers, which is frequently mentioned in the literature, is unfounded. This false impression arises from the many simplistic comparisons made between only two groups of glaciers that are of different types and have very different responses to climate: ‘Frontal changes at the steep and exceptionally responsive Franz Josef and Fox Glaciers, which provide enhanced responses to climate within 5 to 8 years (Fizharris and others, 1999), have repeatedly been compared with the large, low-gradient valley glaciers of the Mount Cook area, which have slow, subdued response times in the order of 100 years’ (Reference ChinnChinn, 2001). Fox and Franz Josef Glaciers have very steep mass-balance gradients and enormous accumulation basins compared with the other valley glaciers presented in Figures 6 and 7. An illustration of their enormous accumulation basins is given in Figure 11 and a comparison is made with the other main valley glaciers west of the Main Divide. The two glaciers with the largest accumulation basins (Fox and Franz Josef Glaciers) have advanced, those with reasonably large basins (e.g. Spencer Glacier, Fritz Glacier, Salisbury Snowfield) seem to be stable, while the valley glaciers with the smallest accumulation basins (e.g. La Perouse and Victoria Glaciers) have shown a marked retreat. It is important to highlight in this context that the hypsometry (i.e. the distribution of glacier area over its altitudinal range) for Fox and Franz Josef Glaciers is also very different from that of the other western valley glaciers. The main part of the accumulation basins for these two glaciers is located above 2000 m a.s.l., while only very small sections of the other valley glaciers are above this elevation. This would make Fox and Franz Josef Glaciers very responsive to small changes in climate, much more so than the other western valley glaciers. Because of their fast reaction time, scientists have tried to correlate the retreat and advance of these two glaciers with climatic parameters. Reference Hooker and FitzharrisHooker and Fitzharris (1999) found a strong link between atmospheric circulation changes, climate variables and the behaviour of Fox and Franz Josef Glaciers. However, causality effects between climate and glacier behaviour are complex, and hypsometry is only one factor in this interplay (Reference Dyurgerov, Meier and BahrDyurgerov and others, 2009). The processes leading to glacier changes were not the subject of this study, but the reader can consult Reference Chinn, Heydenrych and SalingerChinn and others (2005) for a comprehensive summary of annual mass balances in New Zealand and their possible causes during the past few decades.

Fig. 11. Comparison of the area changes (1978–2002) measured for the glaciers west of the Main Divide with delineation of the accumulation basins.