2.1. Overview

Like previous automated tools (Reference Racoviteanu, Paul, Raup, Khalsa and ArmstrongRacoviteanu and others, 2009), the core of our algorithm is based on hydrological modeling tools readily implemented in many established Geographic Information Systems (GIS). These tools map the watershed area that contributes surface flow to a common drainage outlet point (‘pour point’ ; e.g. the site of a river gauging station). By relying on watershed algorithms we implicitly assume that ice divides, topographic divides and hydrological divides are identical, which holds true for most glaciers in confined valleys (Reference Manley, Williams and FerrignoManley, 2008; Reference Schiefer, Menounos and WheateSchiefer and others, 2008; Reference Bolch, Menounos and WheateBolch and others, 2010).

For the identification of pour points, there are differences between glacierized and unglacierized basins that must be considered. In unglacierized river basins, pour points generally coincide with the lowest point of the valley cross section, which makes the calculation of watersheds straightforward (Fig. 1a). However, due to the convex surface of glacier ablation areas, not all ice of the same glacier converges to a single point. Thus, taking the lowest glacier point as the pour point does not necessarily capture the entire surface area of a glacier (Fig. 1b). Hence, multiple pour points may be required to capture the entire glacier surface (Fig. 1c). The identified ‘flowsheds’, one flowshed per pour point, also need to be allocated to individual glaciers in order to obtain the correct outline of each individual glacier within a glacier complex (Fig. 1d). Our algorithm automatically identifies pour points and merges flowsheds associated with individual glaciers.

Fig. 1. Principle of basin separation using hydrological modeling tools. (a) Watersheds (colored areas) calculated in unglacierized terrain using the lowest point of the valley cross section as a single pour point (black circles). (b) Flowsheds in glacierized terrain (colored areas) obtained if only the lowest-lying glacier gridcell is used as pour point. (c) Location of pour points and corresponding flowsheds necessary to capture the entire surface area of the glacier complex. (d) The same situation as in (c) but dashed boxes indicate which pour points belong to the same glacier. Black arrows mark the direction of ice flow.

2.2. Glacier separation

The glacier separation workflow consists of five main steps: (1) preprocessing, (2) identification of pour points, (3) flowshed calculation, (4) allocation of flowsheds to individual glaciers and (5) identification of sliver polygons. Figure 2 illustrates the conceptual workflow. Some of the required empirical parameters have fixed values, others are obtained through calibration.

Fig. 2. Diagram illustrating the conceptual workflow of the separation algorithm. Ovals represent input/output data while rectangles show the main operations. Solid arrows indicate the main progression, and dashed arrows show auxiliary steps.

2.3. Error assessment

The error assessment consists of two main steps: (1) preparation of reference outlines and (2) error allocation. We quantify errors if glaciers are split or merged incorrectly. We do not assess errors that are only due to DEM inaccuracies. Inaccurate DEMs lead to topographic divides that are different from the true topographic divides, which ultimately results in shifted glacier divides. Because such errors do not directly reflect algorithm failure, they are not assessed here.

2.4. Implementation

The glacier separation algorithm is written in the Python™ programming language and uses functions of the Environmental Systems Research Institute ESRI^® ArcGIS 10.1 software package. This algorithm can be run from an integrated development environment or from inside ArcGIS. For larger regions we use a parent script that splits the domain into several subregions and launches the actual script for each of these subregions. The error assessment workflow is also implemented in a Python™ script.

3.1. Calibration

In total, ten parameters (Table 2) were derived for the Juneau Icefield area (Fig. 5a), applying trial-and-error experiments as well as automatized calibration. We used a DEM resampled to a spatial resolution of 40 m. Our experiments showed that this cell spacing retains large-scale relief features and allows for deriving reasonably accurate and smooth flowsheds. For buffer 1, we used 160 m (four cells), and for buffer 2 we chose 80 m (two cells) in order to make sure that none of the glacierized terrain is affected by the creation of the gutter. Trial-and-error procedures showed that an L _gutter of 100 m was sufficient to significantly reduce the number of pour points and to constrain them to lower glacier portions. The parameters A _sliver1, A _sliver2 and C _sliver were determined through visual inspection of raw tool output. In the case of A _sliver1, we consulted reference values from previous studies. Our A _sliver1 of 100 000 m² follows Reference Schiefer, Menounos and WheateSchiefer and others (2008) and is higher than corresponding values used in other studies. For example, Reference Bolch, Menounos and WheateBolch and others (2010) and Reference Frey, Paul and StrozziFrey and others (2012) used thresholds of 50 000 m² and 20 000 m², respectively.

The above parameters were kept fixed during the next sequence of steps initiated by defining the radius R in terms of the flow accumulation, F. A linear relationship between R and F proved unfeasible due to the large spread of the flow accumulation values. In the case of the Juneau Icefield area, the flow accumulation of the obtained pour points ranged between <100 and >428 000 cells, which suggested an exponential or logarithmic relationship between R and F. We ultimately used the exponential Eqn (1) with parameters a and b. Threshold c was introduced concurrently to constrain the maximum value of R. We set the c threshold to 3500 m and proceeded with determining the parameters a and b. By visual inspection, we derived two preliminary a and b values. In practice, we ran the glacier separation algorithm until step 3, ‘flowshed calculation’ (Fig. 2). Next we measured minimal radii required to merge flowsheds that belong together, and maximal radii allowed to keep separate flowsheds apart. This resulted in pour points that had flow accumulation, radius and category (‘maximum’, ‘minimum’ ) values allocated. Visually fitting a curve to the resulting point cloud led to the preliminary values for a and b. For the actual calibration of a and b, we used a workflow similar to the error assessment. We first created a set of reference glacier outlines by running the glacier separation algorithm using the preliminary a and b values. We then visually checked and manually adjusted these outlines. Next, we ran the code multiple times, varying the parameters a and b. Each time, we determined the number of errors using Eqn (2). Figure 6 shows the resulting error surfaces with the a and b parameters on the x- and y-axes. The panels distinguish the error categories. A smaller P _glac, i.e. smaller a and b values, leads to a decreasing number of errors in the category ‘merged incorrectly’ (Fig. 6a). In fact, as P _glac decreases, the number of mergers decreases overall, including wrong mergers. Concurrently, the number of errors in the category ‘split incorrectly’ increases (Fig. 6b), because a smaller P _glac leaves flowsheds split that should be merged. Figure 6c (total misclassifications) illustrates this compromise involved in choosing optimal a and b pairs. We selected 14.3 m for parameter a, and 0.5 for b, which corresponds to a local minimum on the error surface in Figure 6c. Using the entire set of calibrated parameters (Table 2), 2.0% of the 1283 glaciers were misclassified in the Juneau Icefield area (Table 3).

Fig. 6. The number of misclassifications obtained in the Juneau Icefield test area as a function of the applied a and b values: (a) category ‘merged incorrectly’, (b) ‘split incorrectly’ and (c) ‘total misclassified’ . Crosses indicate individual test runs. Error surfaces are interpolated from the number of errors at each cross. Contours are added in (c) as an additional reference. The black circles show the combination of parameters a and b (14.3 m and 0.5) that yielded a minimum of 26 errors.

3.2. Validation

The calibrated parameters (Table 2) were ultimately applied to other regions of Alaska (Western Chugach Mountains, Western Alaska Range, Central Alaska Range, Eastern Alaska Range, Stikine Icefield) (Fig. 5a) and southern Arctic Canada (Fig. 5b) to validate the algorithm. For the error assessment of each area, we created a set of visually checked reference outlines. The automatically derived outlines were then compared to the reference outlines, and errors were determined according to Eqn (2).

As an example, Figure 7 shows the individual glaciers for a subregion of the Western Alaska Range, indicating that glacier basins are properly separated. Overall, 1.9% out of 8121 glaciers were misclassified in Alaska’s validation regions (Table 3). For southern Arctic Canada, 2.9% of the 7537 glaciers were misclassified.

Fig.7. Color-coded glacier outlines automatically derived for a subregion of the Western Alaska Range. The white dots are the pour points. Cross-hatched polygons indicate the P _glac used to merge separate flowsheds. The enlarged insets (b) and (c) have white 50 m contours added that can be checked against the glacier outlines for validation. See Table 1 for sources of the original glacier complex outlines.

There are five cases where the algorithm fails and leads to misclassified glaciers no matter the quality of DEM or glacier complex outlines. All these errors are considered in the error assessment (Table 3), and typical examples are illustrated in Figure 8. In Figure 8a, the algorithm incorrectly splits the glacier complex into two glaciers because the radii R are too small and do not cover both flowsheds. The opposite is illustrated in Figure 8b. Here the algorithm incorrectly merges two glaciers because R is too large and merges two flowsheds that should remain separate. Figure 8c illustrates that glaciers ending within nunataks are not separated at all because the algorithm does not identify pour points along the borders of nunataks. The corresponding misclassifications contribute errors to the category ‘merged incorrectly’ in Table 3. In Figure 8d, the algorithm fails to split the glacier complex into two glaciers because both glaciers have a shared accumulation area while, at the same time, the northern glacier reaches the southern glacier within the distance of buffer 1. In this case, no flowsheds are calculated for the northern glacier individually as its flow is captured by the pour points of the southern glacier. These misclassifications also contribute errors to the category ‘merged incorrectly’ . Figure 8e illustrates another reason why the algorithm may incorrectly split a glacier complex into too many glaciers. Columbia Glacier, shown in yellow, drains into multiple watersheds. If this occurs in the ablation area, such ice masses are generally not considered to constitute separate glaciers. The algorithm is not able to capture these cases because it treats all topographical divides equally, whether or not they are located in the glacier ablation areas.

Fig. 8. Glacier basins (color-coded) for selected subregions of the Eastern Alaska Range (a, b), the Juneau Icefield (c) and the Western Chugach Mountains (d, e). The arrows indicate cases of misclassification. White dots indicate pour points, and the cross-hatched polygons are the P _glac used to merge the flowsheds of individual glaciers. The black polygons are manually adjusted glacier polygons from the reference glacier outlines. The white 50 m contours in (e) illustrate the topographical divide within the ablation area of Columbia Glacier. Table 1 states the sources of the original glacier complex outlines. All the maps have the same scale.

Glacier separation also fails if small portions of glacier complexes reach over mountain ridges into different watersheds. In this case the algorithm correctly divides the complex into separate bodies, and sliver polygons are created. Although the problems with slivers are not technically an algorithm failure, we consider them in the error assessment. According to our validation, undetected slivers are the main contributor to the category ‘split incorrectly’ and thus an important contributor to the total number of errors (Table 3).

Manual intervention is finally required if the used DEM has poor quality. In this case, the derived topographic divides do not coincide with the true topographic divides. Consequences are most noticeable in flat accumulation areas, where small DEM errors have large impacts on the position of the ice divides (Fig. 9). Our error assessments and thus Table 3 exclude this error source because the errors are not related to a failure of the algorithm but rather to quality issues with the input data. Our assessments also exclude errors that occur due to flow divides that are not identical to the true topographic divides. As with errors due to inaccurate DEMs, flat accumulation areas are most susceptible to this last source of error.

Fig. 9. Failed automatic glacier separations in the Juneau Icefield area. Colors indicate individual glaciers separated by our algorithm using the ASTER GDEM2 while black lines define the glaciers using the same algorithm with a more reliable DEM (SRTM). Large discrepancies occur in the eastern part of the domain (annotated by arrows) due to the poor quality of the ASTER GDEM2 in this part of the accumulation area.

4.1. Algorithm performance

With success rates of ~98% (Alaska test areas) and ~97% (Arctic Canada test area), the algorithm shows a good overall performance (Table 3). The percentage of misclassified glaciers ranges between 1.2% (Western Alaska Range) and 2.9% (southern Arctic Canada). It is remarkable that two of the validation areas (Eastern and Western Alaska Range) have higher success rates than the calibration area, the Juneau Icefield. This is mainly because the Juneau Icefield contains a high number of complex geometries (such as glaciers branching in the ablation area) which lead to more misclassifications. Also, the DEM used in the Juneau Icefield area (ASTER GDEM2) has a relatively low quality. The lowest success rate, in southern Arctic Canada, is most likely due to the prevalent complex glacier geometries, in conjunction with the lack of pronounced topographic relief. Moreover, DEM and glacier complex outlines for this region have the lowest quality of all the DEMs and outlines used in this study.

For most regions, the number of incorrectly split glaciers is disproportionately higher than the number of incorrectly merged glaciers (Table 3). At the same time, the median area of the incorrectly split glaciers is much lower than the median area of the incorrectly merged glaciers. This is consistent with our finding that undetected sliver polygons are one of the main contributors to the total error number.

We found that the algorithm fails to split glacier complexes into separate glaciers if pour points are not identified properly (failure of step 2, Fig. 2). This can occur if one glacier’s tongue reaches another’s within the distance of buffer 1, or if glaciers are located within nunataks along which the algorithm does not identify pour points. Although both cases are small contributors to the total error number, their elimination would further improve the performance of the algorithm.

A glacier complex is typically split into too many glaciers if the calculated flowsheds do not comply with the typical definition of a glacier (failure of steps 3 and 5). This occurs, for example, because small portions of glaciers often reach into neighboring watersheds, resulting in sliver polygons. The fact that slivers are the main contributor to the total number of errors indicates that our approach of using thresholds A _sliver1, A _sliver2 and C _sliver is only partially successful in detecting sliver polygons. However, using A _sliver2 in conjunction with C _sliver, we were able to reduce the number of slivers compared to the number generated by using A _sliver1 only. We found that many slivers occur due to artifacts in DEM or glacier complex outlines. In the future, as DEM and glacier complex outlines continue to improve, we anticipate a reduction in the number of slivers. A glacier complex is split into too many or too few glaciers if the flowsheds are merged incorrectly (failure of step 4, Fig. 2). This may occur if the derived pour points have atypical flow accumulation values. Cases of atypically low flow accumulation values and thus small R values are found for glaciers that have a high number of pour points. These cases are rare in our test areas, in part because the gutter significantly reduces the number of pour points. Nevertheless, the occasional failure of step 4 shows that our algorithm cannot deal correctly with the complex nature of all possible glacier shapes and topographies. At this stage, the algorithm needs a small amount of manual correction to obtain optimal results.

Significant manual intervention may be needed if the DEM used is of low quality. Consequences are most pronounced in flat accumulation areas, where small DEM errors have large impacts on the position of ice divides. While certain DEM products have sufficient quality to obtain reasonably accurate flowsheds, other DEM products should be used with care as input for our tool. The Shuttle Radar Topography Mission (SRTM) DEM, for example, produced reasonably accurate divides in our test areas. According to Reference Frey, Paul and StrozziFrey and others (2012), however, the use of the same DEM product is more problematic in other glacierized regions of the world. As the future will bring higher-quality DEMs, the influence of this error source will likely decrease. It is particularly promising that techniques such as interferometric synthetic aperture radar (InSAR) yield good DEMs even in low-contrast glacier accumulation areas (Reference Frey and PaulFrey and Paul, 2012).

Our error assessment does not determine the algorithm’s performance for different glacier types. However, it is plausible, and confirmed by visual inspections, that valley glaciers are most easily identifiable for our algorithm while ice caps pose the most challenges. Figure 10 shows a mixed ice-cap/valley-glacier complex in southern Arctic Canada as separated by our algorithm. Reference Svoboda and PaulSvoboda and Paul (2009, fig. 7b) provide a manually derived version of the same area. We compare the two solutions, although the DEMs and glacier complex outlines used in the two studies are not identical. While their manual approach derives four distinct glaciers overall, our approach derives seven glaciers. Reference Svoboda and PaulSvoboda and Paul (2009) combine glaciers 1, 2 and 3, which illustrates the fact that our algorithm divides ice caps into more sections than most manual approaches. The lobes of the ice cap are sufficiently large for our algorithm to detect separate glaciers. Reference Svoboda and PaulSvoboda and Paul (2009) also combine our glaciers 5 and 6 (Fig. 10). Neither solution is ‘wrong’ ; however, they illustrate the subjectivity inherent in glacier separation and the difficulties of quantitatively assessing the performance of an automatic algorithm.

Fig. 10. Glacier separation of a mixed ice-cap/valley-glacier complex in southern Arctic Canada. This figure corresponds to figure 7b in Reference Svoboda and PaulSvoboda and Paul (2009). In their manual solution, the glaciers annotated with 1, 2 and 3 as well as 5 and 6 are merged.

4.2. Transferability of parameters

Our method of applying one parameter set (Table 2) to entire regions proved to be robust for our test areas in Alaska and Arctic Canada. Although our test areas comprise a wide range of glacier geometries, more testing is required to determine transferability of parameters.

If optimization of parameters should become necessary for a new application area, six of the ten parameters should be considered for adaptation: a, b, c, A _sliver1, A _sliver2 and C _sliver. However, algorithm success is not equally sensitive to all six parameters. Figure 6, for example, illustrates that a small perturbation of b considerably changes the number of derived glaciers in our calibration area, while the same perturbation of a has a smaller influence. Parameter c has a limited influence if varied within a range (i.e. tens of meters) around the value specified in Table 2, because the R values of only a few large flowsheds are actually affected. Finally, it is plausible that the variation of the A _sliver and C _sliver parameters can significantly change the number of derived glaciers. Notably, C _sliver is sensitive to the shape of glacier complex outlines. Jagged outlines derived directly from raster data have lower compactness values compared to similar outlines with a smoothing filter applied. If outlines are derived entirely by hand, the level of digitized detail can also influence the compactness values. As a consequence, the same value for C _sliver can identify different numbers of sliver polygons only because different techniques have been used to derive the glacier complex outlines.

We recommend keeping the parameters DEM resolution, buffer 1, buffer 2 and L _gutter fixed, because they are strongly interrelated. L _gutter, for example, is optimized for the specified buffer widths. A change of the DEM resolution would also implicitly require a recalibration of a and b, because these parameters are optimized for the spatial resolution of 40 m. Varying the spatial resolution changes the flow accumulation, F, at the pour point and thus affects the R values computed in Eqn (1). We used a DEM resolution of 40 m because this cell spacing retains large-scale relief features and allows for deriving reasonably accurate flowsheds. Even if the used DEM had excellent quality, a higher spatial resolution would not significantly improve the quality of the glacier basins while the processing time increased considerably. As most large-scale DEMs used as input for this tool have a native resolution on the order of 40 m, any resampling to a higher spatial resolution would be impractical, leading to oversampling.

4.3. Error assessment

Our error assessment determined the number of errors semiautomatically. The results are reproducible provided reference outlines are available. Clearly, the availability of accurate reference outlines is the main limitation on our approach. In the present study, the reference outlines were obtained by visually scrutinizing outlines from the same algorithm that was eventually assessed using the checked reference outlines. As a consequence, the used reference outlines are not fully independent. The conducted visual checks are also partly subjective, as dividing glaciers can be subjective, despite the explicit definition of ‘a glacier’ in the literature (e.g. Reference Racoviteanu, Paul, Raup, Khalsa and ArmstrongRacoviteanu and others, 2009). Given the large size of the test areas, errors may also be missed during visual checks. We aimed at reducing subjectivity by checking the reference outlines repeatedly and carefully. In the present study, we did not assess errors associated with technician interpretation. Ideally, however, future work should incorporate results from multiple interpreters.

We did not use available independent outlines because of the different techniques and standards used during derivation of these reference outlines. For example, Reference Le Bris, Paul, Frey and BolchLe Bris and others (2011) published semi-automatically derived and manually checked outlines for the Western Chugach Mountains containing a high number of very small glaciers. These small glaciers occur because this study uses a different approach to address sliver polygons. Using their outlines for reference yields 350 misclassified glaciers (19.6% of the total number of glaciers), although larger glaciers are separated nearly identically. Clearly, the high number of errors does not reflect the general agreement between the two datasets.

Using the error area (i.e. the summed area of the misclassified glaciers) instead of the error number would be another way of quantifying the algorithm’s performance. However, we consider the error number more useful because it better reflects the amount of work required to adjust the dataset. The error area may be useful in the above case of the Western Chugach Mountains to show that the two datasets generally have a good agreement. In this particular case, the small overall area of the misclassified glaciers would reflect the good agreement. In general, however, the error area is not considered suitable to quantify the algorithm performance, mainly because it is very sensitive to outliers. One large incorrectly split glacier significantly increases the error area, not reflecting the amount of work required to adjust this error. In fact, misclassifications of large glaciers are particularly easy to identify through visual inspection and thus straightforward to repair.

Our error assessment excludes misclassifications due to DEM inaccuracies that cause derived topographic divides to be different from the true topographic divides. Although adjustments of these errors can be very time-consuming, we exclude these errors because they are not directly related to a failure of the algorithm. Moreover, a reproducible quantification of these errors would be difficult as the required adjustments involve a shift of glacier boundaries only, as opposed to the merging and splitting of entire glaciers required to adjust the errors included in our error assessment. Also, to obtain suitable reference divides, either a high-quality DEM or velocity fields would be required. To date, none of these data are available for our test areas.

4.4. Previous algorithms

Quantitative comparisons of our algorithm to previous algorithms are complicated by the fact that these algorithms are only briefly described as part of broader inventory studies. Moreover, none of the publications contain extended error assessments or examples of the raw tool output that could be compared to the output of our algorithm.

Reference Manley, Williams and FerrignoManley (2008) published the first semi-automatic separation algorithm and applied it to the glacier complexes of the Eastern Alaska Range. His work, also published on the GLIMS website (http://glims.colorado.edu/tools/icedivide_algorithms), established ideas that have been used in more recent publications, including the algorithm presented here. To summarize, Reference Manley, Williams and FerrignoManley (2008) uses a DEM to calculate the median elevation of each individual glacier complex. Every gridcell below the median glacier elevation is considered a pour point. The Reference Manley, Williams and FerrignoManley (2008) algorithm uses a now superseded scripting language, and we therefore were unable to test it against our approach. However, based on our understanding of the method, we speculate that misclassifications would occur because pour points located above the median glacier elevation would be missed. Moreover, as the median glacier elevations of glacier complexes vary, the glaciers of neighboring complexes would be treated differently. For example, a tongue would be identified as such if it is part of a glacier complex that has a high median elevation. A similarly shaped tongue would be missed if it is part of another complex that has a lower median elevation. These potential limitations suggest that other ways to determine pour points (e.g. by identifying local elevation minima (Reference Schiefer, Menounos and WheateSchiefer and others, 2008) or flow accumulation maxima (this study)) may be more promising.

Reference Schiefer, Menounos and WheateSchiefer and others (2008) identify pour points by searching for local elevation minima along the outlines of glacier complexes. If the relief between pour points exceeds a predefined elevation threshold, they are considered to belong to separate glaciers. Reference Schiefer, Menounos and WheateSchiefer and others (2008) optimize the elevation threshold for their study area in British Columbia, western Canada, which is characterized by pronounced topography. Their 250 m threshold represents a compromise between a lower threshold that identifies multiple termini along undulating glacier tongues and a higher threshold that fails to detect smaller glaciers. Manual checks carried out within the present study suggest that the resulting glacier outlines are very sensitive to the choice of the elevation threshold. In addition, the DEM quality is important, as one erroneously high cell is sufficient to raise the relief above the threshold. We speculate that the performance of our approach may be less susceptible to local topography and DEM quality than the approach of Reference Schiefer, Menounos and WheateSchiefer and others (2008).

A third approach, by Reference Bolch, Menounos and WheateBolch and others (2010), is built around a fully automated watershed function from the ArcGIS software package. This function determines pour points automatically and outputs the watersheds. In their study area in western Canada, Reference Bolch, Menounos and WheateBolch and others (2010) obtain best results by running this function not on the outlines of the glacier complexes directly, but on a buffer 1000–1500 m outside the outlines. This has the advantage that the derived basins are not linked to the glacier morphology for a given period of time, so they can be used for different time periods and glacier extents. However, the large buffer distance also implies that the derived basins can contain more than one glacier. Two studies that applied this algorithm in western Alaska (Reference Le Bris, Paul, Frey and BolchLe Bris and others, 2011) and in the western Himalaya (Reference Frey, Paul and StrozziFrey and others, 2012) found that the algorithm can generate numerous artificial polygons that have to be merged manually. Our algorithm addresses these problems using L _gutter and P _glac, which should reduce the amount of manual intervention overall.