Skip to main content Accessibility help


  • Access
  • Cited by 16
  • Cited by
    This article has been cited by the following publications. This list is generated based on data provided by CrossRef.

    Moon, T. Sutherland, D. A. Carroll, D. Felikson, D. Kehrl, L. and Straneo, F. 2018. Subsurface iceberg melt key to Greenland fjord freshwater budget. Nature Geoscience, Vol. 11, Issue. 1, p. 49.

    Mölg, Nico Bolch, Tobias Rastner, Philipp Strozzi, Tazio and Paul, Frank 2018. A consistent glacier inventory for Karakoram and Pamir derived from Landsat data: distribution of debris cover and mapping challenges. Earth System Science Data, Vol. 10, Issue. 4, p. 1807.

    Idalino, Filipe Daros Rosa, Kátia Kellem da Acuña, Francisco Ferrando Kozhikkodan Veettil, Bijeesh Simões, Jefferson Cardia and Souza, Enoil 2018. Recent glacier variations on Mount Melimoyu (44°50'S-72°51'W), Chilean Patagonia, using Sentinel-2 data. Geocarto International, p. 1.

    Bolch, Tobias and Loibl, David 2018. Comprehensive Geographic Information Systems. p. 112.

    Mir, Riyaz Ahmad 2018. Recent changes of two parts of Kolahoi Glacier and its controlling factors in Kashmir basin, western Himalaya. Remote Sensing Applications: Society and Environment, Vol. 11, Issue. , p. 265.

    Rastner, Philipp Strozzi, Tazio and Paul, Frank 2017. Fusion of Multi-Source Satellite Data and DEMs to Create a New Glacier Inventory for Novaya Zemlya. Remote Sensing, Vol. 9, Issue. 11, p. 1122.

    Chudley, Thomas R. Miles, Evan S. and Willis, Ian C. 2017. Glacier characteristics and retreat between 1991 and 2014 in the Ladakh Range, Jammu and Kashmir. Remote Sensing Letters, Vol. 8, Issue. 6, p. 518.

    Huber, Jacqueline Cook, Alison J. Paul, Frank and Zemp, Michael 2017. A complete glacier inventory of the Antarctic Peninsula based on Landsat 7 images from 2000 to 2002 and other preexisting data sets. Earth System Science Data, Vol. 9, Issue. 1, p. 115.

    Paul, Frank Bolch, Tobias Briggs, Kate Kääb, Andreas McMillan, Malcolm McNabb, Robert Nagler, Thomas Nuth, Christopher Rastner, Philipp Strozzi, Tazio and Wuite, Jan 2017. Error sources and guidelines for quality assessment of glacier area, elevation change, and velocity products derived from satellite data in the Glaciers_cci project. Remote Sensing of Environment, Vol. 203, Issue. , p. 256.

    Earl, Lucas and Gardner, Alex 2016. A satellite-derived glacier inventory for North Asia. Annals of Glaciology, Vol. 57, Issue. 71, p. 50.

    Wu, Kun-peng Liu, Shi-yin Guo, Wan-qin Wei, Jun-feng Xu, Jun-li Bao, Wei-jia and Yao, Xiao-jun 2016. Glacier change in the western Nyainqentanglha Range, Tibetan Plateau using historical maps and Landsat imagery: 1970-2014. Journal of Mountain Science, Vol. 13, Issue. 8, p. 1358.

    Zhao, Liyun Ding, Ran and Moore, John C. 2016. The High Mountain Asia glacier contribution to sea-level rise from 2000 to 2050. Annals of Glaciology, Vol. 57, Issue. 71, p. 223.

    Pitcher, Lincoln H Smith, Laurence C. Gleason, Colin J. and Yang, Kang 2016. CryoSheds: a GIS modeling framework for delineating land-ice watersheds for the Greenland Ice Sheet. GIScience & Remote Sensing, Vol. 53, Issue. 6, p. 707.

    Shahgedanova, M. Nosenko, G. Kutuzov, S. Rototaeva, O. and Khromova, T. 2014. Deglaciation of the Caucasus Mountains, Russia/Georgia, in the 21st century observed with ASTER satellite imagery and aerial photography. The Cryosphere, Vol. 8, Issue. 6, p. 2367.

    Pfeffer, W. Tad Arendt, Anthony A. Bliss, Andrew Bolch, Tobias Cogley, J. Graham Gardner, Alex S. Hagen, Jon-Ove Hock, Regine Kaser, Georg Kienholz, Christian Miles, Evan S. Moholdt, Geir Mölg, Nico Paul, Frank Radić, Valentina Rastner, Philipp Raup, Bruce H. Rich, Justin and Sharp, Martin J. 2014. The Randolph Glacier Inventory: a globally complete inventory of glaciers. Journal of Glaciology, Vol. 60, Issue. 221, p. 537.

    Paul, Frank and Mölg, Nico 2014. Hasty retreat of glaciers in northern Patagonia from 1985 to 2011. Journal of Glaciology, Vol. 60, Issue. 224, p. 1033.




      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        A new semi-automatic approach for dividing glacier complexes into individual glaciers
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        A new semi-automatic approach for dividing glacier complexes into individual glaciers
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        A new semi-automatic approach for dividing glacier complexes into individual glaciers
        Available formats
Export citation


Many glaciological and hydrological studies require outlines of individual glaciers rather than total ice cover. Here we develop a new semiautomatic algorithm that uses a digital elevation model (DEM) and outlines of glacier complexes to calculate the extents of individual glaciers. The algorithm first applies hydrological modeling tools to a modified DEM to calculate flowsheds. It then merges flowsheds that belong to individual glaciers using a distance-based approach, whose required empirical parameters are derived from the Juneau Icefield area in Alaska. In this region, 2% of ∼1300 glaciers were misclassified. The algorithm was validated on >25 000 km2 of ice in other regions in Alaska and on >40 000 km2 of ice in Arctic Canada, resulting in ∼2% and ∼3% misclassified glaciers, respectively. Results indicate that the algorithm is robust provided the DEM and the outlines are of good quality.

1. Introduction

Considerable progress has been made towards delineating the extent of the Earth’s glaciers, mainly within the scope of the Global Land Ice Measurements from Space (GLIMS) program (Kargel and others, 2005; Raup and others, 2007). Glacier outlines are typically derived from optical space- or airborne data by (semi-)automated techniques or by manual digitization (Paul and others, 2002, 2004; Raup and others, 2007; Racoviteanu and others, 2009). The first step in assembling glacier outlines is to identify the presence or absence of glacier ice, which in turn allows the computation of total ice cover. The resulting outlines often include glacier complexes, i.e. collections of glaciers that meet at ice divides.

Many studies require not only knowledge of total regional ice extent, but also the outlines and areas of individual glaciers. For example, individual glacier outlines are required for accurate extrapolation of point geodetic or in situ data to arrive at glacier-wide mass balances. Volume–area scaling techniques (Bahr and others, 1997) used to assess regional or global ice volumes (Radić and Hock, 2010) or to project future glacier evolutions (Radić and Hock, 2011) also require accurate glacier areas. Individual glacier outlines are also necessary for physically based approaches used to calculate regional (Linsbauer and others, 2012) or global ice volumes (Huss and Farinotti, 2012). Beedle and others (2008) quantify the effects on glaciological applications that arise from using different glacier outlines.

Manual glacier separation by visual inspection of ice flow patterns is time-consuming and difficult, especially in terrain with low slope. To date, three semi-automatic algorithms have been developed to separate glacier complexes, all of which rely on a digital elevation model (DEM) and hydrological modeling tools (Manley, 2008; Schiefer and others, 2008; Bolch and others, 2010; see Racoviteanu and others, 2009, for a review). Because these algorithms are part of broader inventory studies, explanations of the algorithms and corresponding error assessments are brief. However, these approaches require manual intervention to some extent, in most cases to correct misclassified glaciers. Here we present a new algorithm that aims at minimizing the amount of manual intervention, with the ultimate goal of having a standardized, automated approach capable of coping with a wide range of glacier types. Similar to previous approaches, our algorithm requires a DEM and outlines of glacier complexes as input. Our method is novel in that it takes additional steps to automate the identification of individual glaciers. We test the method on glacier complexes in Alaska, USA, and southern Arctic Canada, and show that for these areas the algorithm performs well and requires little manual correction.

2. Methods

2.1. Overview

Like previous automated tools (Racoviteanu 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 (Manley, 2008; Schiefer and others, 2008; Bolch 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.


We begin with a DEM of a region encompassing all terrain surrounding the glacier complex of interest. As the native resolution and quality of provided DEMs varies, we resample the DEM to a fixed resolution and apply a smoothing filter. In order to constrain our subsequent watershed calculations to the glacierized terrain, it is necessary to extract only those DEM gridcells contained within the outline of the glacier complex. Before performing this extraction, we create a new outline to extend the glacier boundary several DEM gridcells beyond the outer edge of the glacier complex polygon, a procedure known as ‘buffering’ in GIS terminology (buffer 1, Fig. 3). The buffering of our glacier complex enables us to build an artificial ‘gutter’ of depth L gutter that extends beyond the outer edge of the glaciers, located at a fixed distance from the glacier complex (buffer 2), and that has a magnitude less than buffer 1. This process forces flow from upstream regions in the watershed to terminate in the gutter, which in turn reduces the number of pour points, and constrains the pour points to lower glacier elevations. We set the gutter some distance outside the actual glacier complex in order to make sure none of the gridcells of the gutter contain any ice. Thereby, we avoid changing the topography of the glacierized terrain.

Fig. 3. Example showing the DEM gridcells of a glacier terminus includingthe surrounding buffer region. Buffer 1 is used toextractthe DEM. Along buffer 2, the grid is lowered artificially by the distance L gutter to form a gutter. Colors indicate flow accumulation within the gutter. Circles mark pour points defined as gridcells with local flow accumulation maxima. A–B is a cross-profile illustrated by the inset diagram, showing the cross-sectional geometry of the gutter.

Identification of pour points

Next we use GIS tools to calculate the accumulated flow to each gridcell of the gutter (Fig. 3). The accumulated flow measures the number of upstream gridcells contributing flow to a certain cell, where high values indicate a large contributing area. The accumulated flow along the gutter is then used to locate pour points that are needed to identify individual glaciers. We exploit the fact that pour points correspond to local flow accumulation maxima. For each cell along the gutter we compare flow accumulation values to the neighboring flow accumulation values. If the local flow accumulation value is higher than the flow accumulation of its neighbors, the cell is considered a pour point.

Flowshed calculation

In the next step, contributing areas (watersheds) are calculated using standard GIS functions. These functions use calculations of flow direction, derived from the DEM, to identify as a single watershed all gridcells that contribute flow to a pour point. We iterate through the pour points in order of decreasing flow accumulation and calculate the watersheds for each pour point individually. During each iteration, we also check for potential pour points located within the calculated watershed. Such points are removed (not used to calculate watersheds) because they lie above the pour point that was used to calculate the watershed. This step prevents the final watersheds from being split unnecessarily. We then convert the final watersheds, represented as gridded maps, to vector polygons, and smooth the polygons to remove jagged edges that occur in the grid-to-vector conversion. As a final step, we overlay the watersheds which extend a short distance beyond the glacier edge with the original glacier complex polygon, yielding the glacierized portion of the watershed (i.e. the flowshed).

Allocation of flowsheds to individual glaciers

The fact that glaciers do not flow in the same manner as liquid water means that standard GIS watershed functions will often produce more than one flowshed per glacier (Fig. 1). Therefore, our next step is to develop a simple method for grouping flowsheds based on pour point topology. Experimenting with numerous distance-based approaches resulted in a workflow that begins with building circles around each pour point (Fig. 4). These circles are eventually used to decide whether adjacent flowsheds need to be merged. The radius of each circle, R (m), is defined as an increasing function of the flow accumulation, F, at that point,


This equation is based on the idea that we need larger radii to identify the flowsheds of larger glaciers. We exploit the fact that the flow accumulation increases with the size of the glacier. An exponential relationship with two scaling parameters a and b is used based on testing carried out during algorithm calibration. A threshold c constrains the radius R to an explicitly defined maximum size.

Next we determine the section of the circle that is covered by glacier ice, P glac (Fig. 4), and check whether one or more adjacent flowsheds are overlain by P glac. The overlain flowsheds are merged as we assume they belong to the same glacier. This procedure is carried out for every pour point separately by incrementally looping through the pour points. In some cases, P glac consists of multiple disconnected sections completely separated by unglacierized terrain (e.g. if R is sufficiently large to reach into disconnected flowsheds). In this case, we proceed only with the section of P glac whose perimeter is closest to the pour point and ignore any other disconnected sections of P glac.

Fig. 4. Location of pour points and corresponding flowsheds. Dashed circles indicate the radius R. The cross-hatched regions indicate the parts of the circles, P glac, that are covered by glacier ice. When two or more flowsheds are overlain by a particular P glac, those flowsheds are merged. For example, F1 and F2 are merged because they are overlain by P glac of P1. P glac of P3 overlaps F3 and F2 (F2 is already merged with F1 due to P1). P glac of P4 overlaps F4 and F5 (the same is the case for P5), so F4 and F5 are also merged. Summarized, F1, F2 and F3 as well as F4 and F5 are merged, which is the targeted result shown in Figure 1d.

Identification of sliver polygons

The previous step produces a dataset with flowsheds allocated to individual glaciers. This dataset typically contains ‘sliver polygons’, i.e. small, elongate features located mostly along mountain ridges. With their small area and characteristic shape, slivers are typically not considered distinct glaciers (e.g. Frey and others, 2012). In fact, sliver polygons are often artifacts due to DEM or outline inaccuracies, or due to relative shifts between the two datasets. To remove sliver polygons, we merge polygons identified as slivers with the neighboring polygon with the longest shared boundary. Previous approaches deal with sliver polygons by setting minimum area thresholds (e.g. Schiefer and others, 2008). We follow this approach and use a value of A sliver1 below which all polygons are identified as slivers. Our analysis of algorithm performance showed that not all slivers were identified this way. Therefore, we use an additional criterion, glacier compactness, which is defined as the perimeter of a circle with the area of the glacier divided by the actual measured glacier perimeter (Allen, 1998). In addition to the above criterion A sliver1, we identify polygons as slivers if they have an area less than the threshold A sliver2 and, concurrently, a compactness parameter less than the threshold C sliver. Using the thresholds A sliver2 and C sliver, we exploit the fact that sliver polygons typically have small areas along with low compactness values due to their elongate shape. Polygons lacking a shared boundary with a neighboring glacier or not fulfilling the above criteria are considered individual glaciers, which maintains the initial area of the glacier complex.

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.

Preparation of reference outlines

Quantifying errors in our separation algorithm requires a set of reference outlines against which we compare our output. While it would be best to use a fully independent dataset, none exist that would not introduce additional uncertainties due to, for example, differences in technician interpretation. Although they are not fully independent, we use reference data obtained by careful manual checking of our own ‘raw’ algorithm output, using contours, shaded relief DEMs and airborne/spaceborne imagery. We visually check and adjust outlines if a glacier tongue remains split, or if glaciers are merged erroneously. Sliver polygons, undetected by the algorithm, are also merged. However, we do not move glacier boundaries in the accumulation areas that are due to DEM errors rather than algorithm failure. The adjusted reference glacier outlines are used to derive a point dataset that contains one centroid point for every glacier, i.e. one point centrally located within the perimeter of each individual glacier. We then use these centroid points as a basis for counting numbers of errors in our dataset.

Error allocation

We combine the above point dataset with the raw glacier outlines and determine the number of points located within the perimeter of each glacier of the raw glacier outlines. This number of points allows us to derive the number of errors,


Where N error is the total number of errors, ni is the number of points within the ith glacier of the raw glacier outlines and z is the total number of glaciers in the dataset. A glacier that contains one point (n = 1) exists in both the raw and the reference glacier outlines, and no error is added to N error. A glacier that contains no points (n = 0) was merged with another glacier during the visual check and one error is added. A glacier that contains more than one point (n > 1) was split into multiple glaciers during the visual check and n − 1 errors are added to N error. For example, if the algorithm merged three glaciers that should be separate according to the reference outlines, three points (n = 3) are found within the perimeter of this specific glacier, and two errors are added to N error. All glacier outlines contributing errors are stored separately by assigning them to two datasets, ‘split incorrectly’ (if n = 0) or ‘merged incorrectly’ (if n > 1), which allows a statistical analysis of these outlines. To obtain relative errors, N error is divided by the total number of glaciers determined by the algorithm.

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. Application

The algorithm was applied to several glacier complexes in Alaska (∼30000 km2; Fig. 5a) and southern Arctic Canada (∼40000 km2; Fig. 5b). The regions include a variety of glacier types ranging from small mountain glaciers to large ice fields and ice caps. Details of the areas and datasets used are given in Table 1. Quality and native resolution of the used outlines and DEMs vary substantially among the regions. For example, Alaska’s National Elevation Dataset (NED) DEM has non-systematic horizontal and vertical offsets and is out of date due to substantial changes in glacier extent and topography that have occurred since its compilation in the 1950s.

Fig. 5. Location of glacier complexes to which the new algorithm was applied: (a) Alaska and (b) southern Arctic Canada. Both maps have the same scale.

Table 1. Test regions and their datasets. Parameters (Table 2) were derived for the region labeled ‘C’ (calibration), while the same parameters where then applied to the regions labeled ‘V’ (validation)

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 m2 follows Schiefer and others (2008) and is higher than corresponding values used in other studies. For example, Bolch and others (2010) and Frey and others (2012) used thresholds of 50 000 m2 and 20 000 m2, respectively.

Table 2. Calibrated parameters used in this study. Buffer 1 is used to clip the DEM, and buffer 2 locates the gutter. L gutter is the terrain lowering along buffer 2. a and b are the parameters used in Eqn (1) to calculate R. The parameter c constrains R to a maximum size. A sliver1 is the threshold area below which a polygon is merged with the glacier with the longest shared boundary no matter the compactness value.A sliver2 is the threshold area below which a polygon is merged with the glacier with the longest shared boundary if the compactness lies below the threshold C sliver

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.

Table 3. Total number of glaciers in the test regions and number/fraction of misclassified glaciers. The misclassified glaciers are attributed to the categories ‘split incorrectly’ and ‘merged incorrectly’ and the median areas are given for each category

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. Discussion

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 Frey 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 (Frey 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. Svoboda 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. Svoboda 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. Svoboda 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 Svoboda 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. Racoviteanu 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, Le 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.

Manley (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 (, established ideas that have been used in more recent publications, including the algorithm presented here. To summarize, Manley (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 Manley (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 (Schiefer and others, 2008) or flow accumulation maxima (this study)) may be more promising.

Schiefer 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. Schiefer 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 Schiefer and others (2008).

A third approach, by Bolch 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, Bolch 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 (Le Bris and others, 2011) and in the western Himalaya (Frey 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.

5. Conclusions

We have developed a new algorithm to separate glacier complexes into individual glaciers. The algorithm is based on hydrological modeling tools and identifies individual glaciers semi-automatically. Application of the algorithm to >65 000 km2 of ice in Alaska (∼98% success rate) and Arctic Canada (∼97% success rate) indicates that the method is robust, requiring only a small number of manual corrections. Most misclassifications are due to sliver polygons, which not only occur due to failure of the algorithm, but also due to inaccuracies of the glacier complex outlines or DEMs. Future refinements of the present algorithm, together with improved DEMs and outlines, are anticipated to further reduce the number of misclassifications. However, given the complicated nature of possible glacier geometries and inaccuracies of DEMs and glacier complex outlines, it will remain challenging to develop a fully automatic approach.

Sophisticated algorithms to split glacier complexes into single glaciers are a crucial link between the derivation of glacier complexes (e.g. from remote-sensing data) and the many applications that require individual glacier outlines as input (e.g. the compilation of glacier inventories). While there has been a wealth of research on both the automatization of glacier complex delineation and glaciological applications, research with regard to the actual glacier separation has been rare. Accordingly, glacier outlines have remained underived for many glacierized areas even though the corresponding glacier complex outlines are available. To help remedy this problem, our code is available for use as well as further development.


Thanks are due to J. Rich, S. Herreid and A. Gardner for providing outlines and DEMs. Constructive comments by H. Frey and two anonymous reviewers improved the manuscript. G. Vance and H. Jiskoot also provided helpful comments. Support for this work was provided by NASA under the Cryosphere Sciences Program (grants NNH10Z1A001N, NNX11AF41G, NNX11A023G) and the US National Science Foundation (grant EAR-0943742).


Allen, TR (1998) Topographic context of glaciers and perennial snowfields, Glacier National Park, Montana. Geomorphology, 21(3–4), 207216 (doi: 10.1016/S0169-555X(97)00059-7)
Bahr, DB, Meier, MF and Peckham, SD (1997) The physical basis of glacier volume–area scaling. J. Geophys. Res., 102(B9), 20 35520 362 (doi: 10.1029/97JB01696)
Beedle, MJ and 7 others (2008) Improving estimation of glacier volume change: a GLIMS case study of Bering Glacier System, Alaska. Cryosphere, 2(1), 3351 (doi: 10.5194/tc-2-33-2008)
Bolch, T, Menounos, B and Wheate, R (2010) Landsat-based inventory of glaciers in western Canada, 1985–2005. Remote Sens. Environ., 114(1), 127137 (doi: 10.1016/j.rse.2009.08.015)
Frey, H and Paul, F (2012) On the suitability of the SRTM DEM and ASTER GDEM for the compilation of topographic parameters in glacier inventories. Int. J. Appl. Earth Obs. Geoinform., 18, 480490 (doi: 10.1016/j.jag.2011.09.020)
Frey, H, Paul, F and Strozzi, T (2012) Compilation of a glacier inventory for the western Himalayas from satellite data: methods, challenges and results. Remote Sens. Environ., 124, 832843 (doi: 10.1016/j.rse.2012.06.020)
Gardner, AS, Moholdt, G, Arendt, A and Wouters, B (2012) Long-term contributions of Baffin and Bylot Island Glaciers to sea level rise: an integrated approach using airborne and satellite laser altimetry, stereoscopic imagery and satellite gravimetry. Cryos. Discuss., 6(2), 15631610 (doi: 10.5194/tcd-6-1563-2012)
Huss, M and Farinotti, D (2012) Distributed ice thickness and volume of all glaciers around the globe. J. Geophys. Res., 117(F4), F04010 (doi: 10.1029/2012JF002523)
Kargel, JS and 16 others (2005) Multispectral imaging contributions to global land ice measurements from space. Remote Sens. Environ., 99(1–2), 187219 (doi: 10.1016/j.rse.2005.07.004)
Le Bris, R, Paul, F, Frey, H and Bolch, T (2011) A new satellite-derived glacier inventory for western Alaska. Ann. Glaciol., 52(59), 135143 (doi: 10.3189/172756411799096303)
Linsbauer, A, Paul, F and Haeberli, W (2012) Modeling glacier thickness distribution and bed topography over entire mountain ranges with GlabTop: application of a fast and robust approach. J. Geophys. Res., 117(F3), F03007 (doi: 10.1029/2011JF002313)
Manley, WF (2008) Geospatial inventory and analysis of glaciers: a case study for the eastern Alaska Range. In Williams, RS Jr and Ferrigno, JG eds. Satellite image atlas of glaciers of the world. (USGS Professional Paper 1386-K) United States Geological Survey, Denver, CO, K424K439
Paul, F, Kääb, A, Maisch, M, Kellenberger, T and Haeberli, W (2002) The new remote-sensing-derived Swiss glacier inventory: I. Methods. Ann. Glaciol., 34, 355361 (doi: 10.3189/172756402781817941)
Paul, F, Huggel, C and Kääb, A (2004) Combining satellite multi-spectral image data and a digital elevation model for mapping debris-covered glaciers. Remote Sens. Environ., 89(4), 510518 (doi: 10.1016/j.rse.2003.11.007)
Racoviteanu, AE, Paul, F, Raup, B, Khalsa, SJS and Armstrong, R (2009) Challenges and recommendations in mapping of glacier parameters from space: results of the 2008 Global Land Ice Measurements from Space (GLIMS) workshop, Boulder, Colorado, USA. Ann. Glaciol., 50(53), 5369 (doi: 10.3189/172756410790595804)
Radić, V and Hock, R (2010) Regional and global volumes of glaciers derived from statistical upscaling of glacier inventory data. J. Geophys. Res., 115(F1), F01010 (doi: 10.1029/2009JF001373)
Radić, V and Hock, R (2011) Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise. Nature Geosci., 4(2), 9194 (doi: 10.1038/ngeo1052)
Raup, B and 11 others (2007) Remote sensing and GIS technology in the Global Land Ice Measurements from Space (GLIMS) Project. Comput. Geosci., 33(1), 104125 (doi: 10.1016/j.cageo.2006.05.015)
Schiefer, E, Menounos, B and Wheate, R (2008) An inventory of morphometric analysis of British Columbia glaciers, Canada. J. Glaciol., 54(186), 551560 (doi: 10.3189/002214308785836995)
Svoboda, F and Paul, F (2009) A new glacier inventory on southern Baffin Island, Canada, from ASTER data: I. Applied methods, challenges and solutions. Ann. Glaciol., 50(53), 1121 (doi: 10.3189/172756410790595912)