The length of the glaciers in the world a straightforward method for the automated calculation of glacier center lines

Glacier length is an important measure of glacier geometry but global inventories are mostly lacking length data. Only recently semi-automated approaches to measure glacier length have been developed and applied regionally. Here we present a ﬁrst global assessment of glacier length using a fully automated method based on 5 glacier surface slope, distance to the glacier margins and a set of trade-o ﬀ functions. The method is developed for East Greenland, evaluated for the same area as well as for Alaska, and eventually applied to all ∼ 200000 glaciers around the globe. The evaluation highlights accurately calculated glacier length where DEM quality is good (East Greenland) and limited precision on low quality DEMs (parts of Alaska). Measured 10 length of very small glaciers is subject to a certain level of ambiguity. The global calculation shows that only about 1.5 % of all glaciers are longer than 10 km with Bering Glacier (Alaska/Canada) being the longest glacier in the world at a length of 196 km. Based on model output we derive global and regional area-length scaling laws. Dif-ferences among regional scaling parameters appear to be related to characteristics of 15 topography and glacier mass balance. The present study adds glacier length as a central parameter to global glacier inventories. Global and regional scaling laws might proof beneﬁcial in conceptual glacier models.


Introduction
Glacier length is one of the central measures representing the geometry of glaciers. 20 Changes in climate have a delayed but very clear impact on glacier length and the advance or retreat of glaciers is frequently used to communicate observed changes to a broader public. In a scientific context, glacier length change records are interpreted with respect to variations in climate (e.g., Hoelzle et al., 2003;Oerlemans, 2005). But despite of being of scientific relevance and easy to communicate, glacier length is 25 2492 basic processes controlling the geometry of glaciers. But defining the longest flow line on a glacier is a non-trivial task because ice forming in the upper accumulation area travels close to the glacier bed towards the tongue. Thus the longest flow line is located somewhere close to the bottom of the glacier while flow trajectories of surface particles do never extend over the full length of a glacier. Then again glacier length generally 10 refers to the glacier surface represented on a map, a satellite image or in reality. Hence, glacier length as a surface measure can only have an indirect representation in the three-dimensional process of glacier flow.
In the past, glacier length was determined manually in a laborious way. Automated computation of glacier length has gained new relevance with the advent of the Ran- 15 dolph Glacier inventory (RGI), a worldwide data set of glacier polygons (Pfeffer et al., 2014). While other geometric parameters such as area, elevation and slope can be automatically derived from glacier polygons and Digital Elevation Models (DEMs), until recently no automated approach existed to measure glacier length. The criteria for such an approach should be based on the definition of glacier length given above but 20 also need to address practical issues: the method should (i) mimic glacier flow, (ii) be computationally efficient, (iii) be fully automated and (iv) needs to be able to deal with inaccurate DEM data. Thereby the requirements (ii) to (vi) result from application to large-scale glacier inventories and limitations in quality of the input data.
Two recent studies by Le Bris and Paul (2013) and Kienholz et al. (2014) presented 25 semi-automatic approaches to derive glacier length and demonstrated the methods in local or regional applications. Le Bris and Paul (2013) suggest determining the highest and lowest point of a glacier. A line connecting the two points is drawn in a way that distance to the glacier margins is maximized and downhill flow is respected. Kienholz  al. (2014) introduce the term "center line" and base their approach on the same principle of maximizing the distance to the glacier margin. Elevation is considered as a second criterion and both conditions are combined by minimizing the costs on a cost grids.
The approach by Le Bris and Paul (2013) has the advantage that limitations in DEM 5 quality have little influence on the center lines. Disadvantages are the restriction to only one center line per glacier that does not necessarily correspond to the longest one. Finally, the method does only work well on certain glacier types. The approach by Kienholz et al. (2014) performs well on most glacier types and each branch of a glacier is represented by its own center line. Calculating several center lines per glacier in- 10 creases the likelihood that the actually longest center line is chosen to represent total glacier length.
Here we present a third approach for the calculation of center lines. In contrast to the aforementioned studies, we design a fully automated approach and apply the method globally. Thereby we aim at closing an important gap in glacier inventories by calcu- 15 lating a length attribute for all glaciers in the world. Our concept firstly relies on hydrological flow which, in the past, was considered to be of limited value to calculate glacier length (Schiefer et al., 2008;Paul et al., 2009;Kienholz et al., 2014). In fact, hydrological flow is a good predictor for glacier length when combined with centrality as a second condition. In this paper, the two conditions of maximizing surface slope 20 angle and centrality are combined in a non-hierarchical way and their weights are flexibly controlled by trade-off functions. The methodology requires glacier polygons and a DEM for input. The output is a set of center lines covering every individual branch of a glacier.
Development and initial validation of the approach was performed for local glaciers 25 of central East Greenland, followed by calculating glacier length for all Alaskan glaciers and comparing our results to Kienholz et al. (2014). Eventually the method is applied to all glaciers of the globe and length characteristics are analyzed on a regional and global scale.
2494 2 Area of application and input data

Test site East Greenland
The center line calculation is developed and tested on the example of a strongly glacierized area in East Greenland. The test site reaches from 68.0 to 72.5 • N and 21.5 to 32.5 • W and represents a transition zone between the Greenland ice sheet and local 5 glaciers. The area was chosen because among its 3950 individual ice bodies all possible morphometric types of glaciers are present: from small cirque glaciers to large valley glaciers and ice caps with marine terminating outlets. The total glacierized area is approximately 41 000 km 2 and comprises the Geikie plateau glaciation of roughly 27 000 km 2 where catchments of individual outlet glaciers reach up to 4200 km 2 in area.

10
North of the Geikie plateau smaller ice caps and valley glacier systems are dominant.
For the East Greenland test site we used the Greenland Mapping Project (GIMP) Digital Elevation Model (DEM) (Howat et al., 2014) at 90 m horizontal resolution. All sinks (i.e. grid cells or clusters of grid cells that are entirely surrounded by cells of higher elevations) were removed from the DEM using a sink-fill algorithm by Planchon 15 and Darboux (2001). The glacier polygons were obtained from Rastner et al. (2012).

Test site Alaska
For the purpose of model comparison we calculated center lines and glacier length for the same perimeter as Kienholz et al. (2014) and using identical DEM data and glacier outlines. Glacier outlines refer to the years 2000-2012 and the DEM with a resolution 20 of 30 m is a composite of the best regionally available data from the Shuttle Radar Topography Mission (SRTM), the Système Pour l'Observation de la Terre (SPOT), the GDEM v2 from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and Alaskan Interferometric synthetic aperture radar (IfSAR) (Kienholz et al., 2014)

Worldwide application
The Randolph Glacier Inventory provides digital outlines for all glaciers around the globe except the two ice sheets of Greenland and Antarctica (Pfeffer et al., 2014). In total, the RGI contains roughly 200 000 individual glaciers with a total area of 726 800 km 2 and is organised in 19 regions. By intersecting all glacier outlines with global DEMs a lo-5 cal terrain model and a glacier mask on a metric grid with a resolution of 25 to 200 m were derived for each individual glacier, see e.g., Huss and Farinotti (2012).
Here, we use the RGIv3.2 released in August 2013. Between 55 • S and 60 • N surface elevation is obtained from the SRTM DEM v4 (Jarvis et al., 2008) with a resolution of about 90 m. At latitudes above 60 • the ASTER GDEM2 v2 (Tachikawa et al., 2011) 10 was used. DEMs for glaciers and ice caps around Greenland are based on the GIMP DEM (Howat et al., 2014), and glaciers in the Antarctic Periphery are mostly covered by the Radarsat Antarctic Mapping Project (RAMP) DEM v2 (Liu et al., 2001) featuring a resolution of 200 m.
For some regions problems with the quality of the ASTER GDEM are documented, 15 in particular for low contrast accumulation areas of Arctic glaciers (Howat et al., 2014). In addition, the RGI still contains a certain number of polygons describing inaccurate or outdated outlines (Pfeffer et al., 2014). Although these limitations only apply to a relatively small number of glaciers (estimated as 1 % of the total) in confined regions they are expected to have a certain influence on calculated glacier length.

The center line model
In the following, the center line computation is explained in detail. The basic concept is schematized in Fig. 1 and visualized on the example of two small mountain glaciers of East Greenland in Figs. 2 and 3. The chosen settings for the model parameters are listed in Appendix Table A1.

Model input
The code performing the fully automated computation is written in IDL. The calculation of center lines is entirely grid-based and relies on two input grids, namely (i) a sinkfilled DEM, and (ii) a gridded mask of glacier polygons. The DEMs for the worldwide computation were smoothed to remove or suppress spurious small scale undulations 5 and subsequently sink filled. The glacier mask is of identical size and cell size (C) as the DEM and represents the glacier polygons in a rasterized form, i.e. each grid cell is assigned the ID of the overlaying glacier polygon. Grid cells outside the glacier polygons are given a no-data value.

Computation principle
Center lines are computed from top to bottom. Starting points of center lines are selected automatically along the part of the glacier margin that is located in the accumulation zone, and on summits. 15 The first criterion requires knowledge of the Equilibrium Line Elevation (ELA) that separates accumulation and ablation zone. Because the ELA is not measured for the vast majority of the worldwide glaciers, it is approximated as the median elevation (z med ) which provides a reasonable representation of the ELA for glaciers where mass loss is mostly restricted to melt (cf. Braithwaite and Raper, 2009). On calving glaciers 20 z med lies above the actual ELA but this overestimation is unlikely to affect the automated measurement of glacier length. From all grid cells located at the glacier margin and above z med every n s th cell is picked as a starting points. The starting points are complemented by summit points, here defined as local topographic maxima within an arbitrarily chosen diameter of d s grid cells. Introduction Beginning from each starting point, a center line is calculated in a stepwise manner. In each calculation step the next center line point is chosen from the cells of a ringshaped and flexibly-sized search buffer (Fig. 2b). The width of the ring is always one grid cell and a circular shape is approximated as well as the currently applied radius allows. The radius of the search buffer depends on the distance to the glacier-margin 5 of the current center line point and is always chosen to be one grid cell smaller than the current distance to margin.
Center lines are continued until a point is reached where grid cells in all directions are either up-slope, at a zero angle or non-glacierized. Thus, center lines find their endpoint autonomously and do not progress to a predefined ending point as it is the case in the 10 approaches by Le Bris and Paul (2013) and Kienholz et al. (2014). An exception are tide-water glaciers where endpoints are automatically suggested, but not prescribed (see Sect. 3.4).

Implementation of the basic conditions and trade-off functions
In each computation step the choice of the next center line point relies on two basic 15 principles: hydrological flow: maximize the slope angle from one center line point to the next.
-Distance to margin: maximize the distance to the glacier margins.
In most computation steps there is no grid cell in the search buffer where both centrality and downhill-slope are at their maximum. The importance of the two basic con-20 ditions also varies with glacier type and specific location on a glacier as explained in detail below. Consequently, the two basic conditions are flexibly weighted by applying a basic trade-off factor c 0 which is in every individual calculation step modified according to three trade-off functions T (c.f. T1 adds weight to slope in the accumulation area of a glacier and gives preference to centrality in the ablation area. The goal is to prioritize centrality on glacier tongues since 5 too much emphasis on slope would make the center lines drift towards the margins as discussed by Schiefer et al. (2008). In the accumulation area, slope receives more weight -otherwise center lines maximize centrality at the cost of running diagonally to elevation contours across the surface. The trade-off function expresses the elevation (z) of the current center line point through a dimensionless factor c 1 : if z is equal to 10 the maximum glacier elevation, then c 1 = f 1 , at median elevation c 1 = 0, and c 1 = −f 1 at minimum elevation. See the Appendix for a detailed description, a complete list of parameter values and explanations on the numerical example given in Fig. 2b.
T2 emphasizes centrality at locations where a glacier has a constant width and adds weight to slope where the glacier changes width. The function reduces the weight 15 of centrality where a glacier suddenly widens. This is, for instance, the case at the confluence of two glacier tongues. The function is also responsible for a more direct course towards the progressively more narrow glacier termini. A dimensionless factor c 2 is calculated from the ratio of the mean distance to margin of all up-hill grid cells (w u ) and the mean distance to margin of all down-hill grid cells w d (see the Appendix 20 for full details).
T3 emphasizes slope in sections where (i) the glacier steepens and (ii) the glacier surface is generally steep. The function achieves a more direct down-hill flow in steep glacier sections and a more direct and slope-controlled course where ice masses from an ice cap progress into outlet glacier tongues. Such locations are often associated 25 with a general steepening of a glacier. A dimensionless factor c 3 is calculated based on (i) the ratio of the mean slope α u and α d above and below, respectively, the current elevation of the flow path, and (ii) based on α d alone (see the Appendix for full details). The basic trade-off factor and the three functions are combined in a non-hierarchical way by c t = 3 i =0 c i . A minimum slope angle is then calculated according to α min = α d c t . Finally from all grid cells fulfilling the α min condition the one at maximum distance from the glacier margin is chosen as the next center line point. If two or more grid cells are at maximum distance-to-margin then the cell with the maximum slope angle is selected 5 as illustrated in Fig. 2b.

Suggesting glacier endpoints
The autonomous selection of endpoints (see Sect. 3.2) results in arbitrarily chosen endpoints on marine-or lake-terminating glaciers with a wide glacier front where even manual definitions of sole glacier endpoints are debatable. Figure 4 illustrates the is-10 sue on the example of two tide-water glaciers. Glacier length could be maximized by measuring at the margins but it appears more logical and consistent to end center lines in the middle of a calving front.
An automated approach is applied to approximate the middle of a calving front and suggest these points as endpoints. A glacier is assumed to be lake-or marine- 15 terminating whenever there is a certain number (n c ) of grid cells that are (i) located at the glacier margin and (ii) within a certain elevation threshold (z c ) of the lowest grid cell of the glacier (c.f. Fig. 4). If these conditions are fulfilled then potential glacier endpoints are determined by performing a neighborhood analysis where all grid cells fulfilling conditions (i) and (ii) are brought into groups of directly adjacent cells. For each group for 20 which the number of members n g exceeds n c the geometric center of the location of all group members is calculated. Finally, the group member closest in distance to the geometric center is chosen as a suggested endpoint. Since there can be several groups exceeding n c in members, a glacier can have more than one suggested endpoint.
Calculation of center lines for glaciers with suggested endpoints is identical to other 25 glaciers with the exception that as soon as a center line has moved to within a certain distance (D e ) of a suggested endpoint, the line is redirected to the endpoint and the center line is terminated (Fig. 1). Thereby D e is defined as the maximum distance-2500 Introduction to-margin value of all glacier cells located within a radius of C · n g from the potential endpoint.

Filtering
Finally two filters are applied to smooth the center lines: -F1: in four iteration steps points are removed that describe an angle of less than 5 θ 1 (θ 2 in the second to fourth iteration) with their two neighboring points.
-F2: a minimal spacing of D f meters between center line points is introduced by deleting points from sections of short spacing between points.
F1 and F2 are applied consecutively. F1 mainly smooths the center lines in sections where the minimal search radius of one grid cell in each direction is applied. Under 10 these circumstances a center line can find its next point only in eight directions. If the overall direction deviates from these eight angles, the center lines cannot progress in a straightforward way. Examples are illustrated in Fig. 2a close to the glacier terminus to the west where the unfiltered center line describes a zigzag pattern. F2 basically reduces the spatial resolution of the center lines. The filter is less important on good 15 quality DEMs as can be seen in Fig. 2a where only marginal changes results. However, the filter is useful in removing some of the irregularities in center lines calculated on low quality DEMs.

Calculating glacier length
For each glacier the same number of center lines is calculated as there are starting 20 points (Fig. 3). The length of each center line is calculated by summing up the distances between all individual points and the length of the longest line is chosen to represent glacier length.

Model calibration and evaluation for East Greenland and Alaska
The approach was calibrated on the example of East Greenland by varying the parameters of the trade-off functions and filters unless center lines were achieved that fulfill the following qualitative criteria: center lines should cross elevation contours perpendicularly,

5
flow strictly downhill, not cut corners, be in the center of the glacier below z med , and end at the lowest glacier point.
The last two criteria are relevant on typical valley-glacier tongues but can be mislead- 10 ing on certain glacier types such as ice caps without outlet glaciers, slope glaciations and cirque-type glaciers. Thus, they are only considered on glaciers where it is assumed that they are in agreement with the characteristics of the actual (imaginary) longest flow line. 15 The calibrated model needs to maintain flexibility to deal with inaccuracies in input data. Figure 2, for instance, exemplifies a very common problem: the western glacier tongue is shifted relative to the DEM. Thus, flow diagonal to contour lines needs to be tolerated although objecting the first quality criteria. The accumulation areas close to the Greenland ice sheet (Fig. 5) show spurious surface undulations. Under such 20 conditions the calculation of reasonable center lines requires a basic ability of "leapfrogging" across smaller undulations at the cost of violating strict downhill flow. The two examples illustrate the aforementioned need for flexibility, but the latter must also be limited to avoid erroneous results where input data is of good quality. 2502

Considering inaccuracies in input data
The two main compromises for flexibility are the following: (i) The search radius is always maximized (cf. Sect. 3.2) to reduce the influence of DEM irregularities. Side effects are an optimization of computation time but also coarse resolution of the center lines where DEM quality would permit a better resolution. (ii) The settings for the tradeoff between slope and centrality (Sect. 3.3) are applied in a way that centrality receives 5 a relatively high weight. For instance the basic parameter c 0 is set to 0.6 (Appendix Table A1). In case none of the trade-off functions takes effect, this means that the minimum required slope α min is only 60 % of the mean slope angle α d of all downhill cells (c.f. Fig. 2). Consequently, only few of the downhill-cells are excluded prior to selecting the cell with maximum distance-to-margin. Slope thus receives less influence 10 and center lines maintain a certain flexibility to move laterally.

Evaluation East Greenland
On average ∼ 22 center lines were calculated per glacier for the East Greenland site, resulting in a total of 88 000 center lines. The longest center line for each glacier of the Geikie plateau is shown in Fig. 5. Realistic center lines result even for glacier polygons 15 of highest complexity and the approach performs well on all types of glaciers including ice caps and marine terminating glaciers. Somewhat erratic center lines appear in the wide accumulation areas close to the ice sheet where DEM quality is comparably low (marked with 1 in Fig. 5). Certain center lines do not start at the apparently most distant point of a glacier (marked with 2) because the surface at that location drains 20 into an adjacent glacier. For the vast majority of glaciers, the center lines end where envisaged, but there are a few locations (one example is marked with 3) where the automatic approach suggests erroneous endpoints. The reason is the complete absence of topography at the glacier tongue which might be the result of different glacier terminus positions in the inventory and the DEM. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | randomly selected. The length of the 100 glaciers was then measured manually while automatic center lines were masked. The averages of automatic and manual glacier length are almost identical (Table 1) while the mean of all glacier specific length ratios R a/m = L a /L m is 1.02 and indicates a small positive bias. The linear regression of L a against L m yields a high correlation (Fig. 6a). Deviations from a perfect agreement 5 (R a/m = 1) are generally small as supported by a root-mean-square deviation (RMSD) of 0.1 (i.e. 10 %, Table 1). When divided into four glacier size classes largest scatter of R a/m is found for glaciers smaller than 0.5 km 2 . Deviations are small for larger glaciers and at minimum for glaciers > 10 km 2 (Fig. 6b). In total there are 14 glaciers (7 of them from the small-10 est size class) where |R a/m − 1| exceeds 0.1. Analyzing the reasons revealed that 6 cases can be attributed to erroneous automatic center lines, in one case a manual center line was deemed wrong upon reconsideration and for the remaining 7 glaciers both the automatic and the manual center lines appear to be equally valid solutions. From the 6 erroneous center lines, 4 are from the smallest size class and one results 15 from inconsistencies in the DEM.

Comparison to glacier length for Alaska
On the example of the Alaska glacier inventory automatically derived glacier length was compared to the semi-automatic measurements (L sa ) by Kienholz et al. (2014). The comparison was done for all 21 720 glaciers exceeding 0.1 km 2 in area and is 20 summarized in Table 2 and in Fig. 7. By average L a is a few percent smaller than L sa and variability of R a/sa = L a /L sa is larger than for East Greenland as indicated by a higher RMSD and more extreme minimum and maximum values of R a/sa . Figure 7a indicates that the largest scatter of R a/sa is found for glaciers smaller than 0.5 km 2 . Figure 7b visualizes  the smallest glacier size class which accounts for 55 % of all investigated glaciers. With increasing glacier size the distribution of R a/sa becomes increasingly centered around R a/sa = 1 (Fig. 7b). In the smallest glacier class 56 % of all R a/sa are within a range of 0.9 to 1.1, the same is the case for 68 % of the glaciers in the second size class, 82 % in the third size class and for 92 % of the glaciers > 10 km 2 . Values of R a/sa > 1.1 are 5 rare for larger glaciers, but there is a relevant fraction of glaciers with R a/sa < 0.9 in all size classes (e.g., 36 % of the smallest glaciers, 7 % of the glaciers > 10 km 2 ). The comparison thus highlights two features: (i) a considerable scatter of R a/sa for small glaciers and (ii) a more general tendency towards R a/sa < 1. Analyzing a number of randomly picked glaciers suggests the following reasons: the deviations on small 10 glaciers are often related to a general ambiguity in defining length. The agreement is highest for elongated features with their longer axis pointing downhill. Ambiguities are large for glacierets located in gently-sloping terrain or for polygons of irregular shape. Automatically measured length of small glaciers is somewhat shorter on average because the automatic approach often adheres more strictly to downhill flow whereas the 15 semi-automatic method has a tendency to cross small polygons diagonally.
These aforementioned issues are of limited relevance to larger glaciers where either L a agrees with L sa or underestimates actual glacier length due to DEM irregularities. In all cases where R a/sa 1, center lines get stopped half-way and L a is eventually measured from lines that do not represent the entire glacier perimeter. Besides these 20 underestimations there are a few cases where DEM irregularities do not stop a center line but force a detour resulting in an overestimation. Surface slope is an important variable in the automatic approach and the method thus suffers from the low DEM quality for certain areas of Alaska. around the globe. Computations are fully automated and no glacier-, glacier-type or region-specific adjustments were conducted. Based on the evaluation in East Greenland and Alaska we estimated typical uncertainties in calculated glacier length. The validation indicated that there is ambiguity in measuring the length of small glaciers and that our approach depends on DEM quality.

5
As a rule of thumb, uncertainty of glacier length of small glaciers, including any ambiguity, is approximately 20 %. On larger glaciers uncertainty in calculated glacier length depends mainly on DEM quality: uncertainty is estimated to be around 2-5 % where the elevation data are reliable and 5-15 % for regions of lower DEM quality. Furthermore, calculated glacier length can be meaningless where glacier polygons are erroneous as 10 it is for instance the case in some areas of Northern Asia (c.f. Pfeffer et al., 2014, for an in-depth discussion of limitations of the RGI).
At a worldwide scale, 3153 glaciers outside of the two ice sheets are longer than 10 km, and 223 are longer than 40 km (Fig. 8). The majority of long glaciers is located in the polar regions (Alaska, Arctic Canada, Greenland, Svalbard, Russian Arc-15 tic, Antarctic) However, there are also more than 500 glaciers with >10 km in length in High Mountain Asia -Fedchenko and Siachen Glacier are more than 70 km long. Bering Glacier, Alaska/Canada, is the longest glacier in the world (196 km). Glaciers in the periphery of Greenland and Antarctica also reach lengths of more than 100 km (Fig. 8). The maximum glacier length in regions dominated by smaller glaciers (Euro-20 pean Alps, Caucasus, New Zealand) is between 10 and 30 km.
Several studies have shown that there is a characteristic scaling between glacier area, volume and length (Bahr et al., 1997;Radic et al., 2008;Lüthi, 2009). In order to analyze the differences in glacier length between individual regions around the globe and to provide a simple mean for estimating glacier length from its area we derive 25 region-specific scaling relationships of the form where L (km) is glacier length along the centerline, A (km 2 ) is glacier area and k (km 1−2β ) and β (−) are parameters.
By least-square linear regression for all ∼ 200 000 area-length pairs k = 2.29 km 1−2β and β = 0.556 were determined. A typical glacier with an area of 1 km 2 thus has a length of 2.3 km and a 10 (100) km 2 glacier can be expected to be 8.2 (29.6) km 5 long. Scaling parameters were also evaluated for four large-scale regions integrating specific glacier properties: (i) Alaska with the largest valley glaciers, (ii) mid-latitude mountain glaciers, (iii) polar regions dominated by ice caps, and the South American Andes (Fig. 9). Correlation coefficients in the log-log space were between r 2 = 0.83 and 0.94. 10 As the differences in glacial morphology and climate are large among the regions the empirical scaling parameters between glacier area and length show some variability (k = [0.85, 3.40], β = [0.467, 0.606]). Expected length for a glacier area of 100 km 2 can thus vary between 18 and 47 km (Fig. 9) between the large-scale regions. The longest glaciers for a given area are found in Alaska. This might be explained by the 15 large elevation differences and the canalizing structure of mountain morphology. Interestingly, mid-latitude glaciers and polar ice caps almost show the same area-length scaling parameters although they strongly differ in shape (Fig. 9). Most likely, different effects of their morphology (average slope and width) and climatology (surface mass balance gradients) compensate for each in other in terms of the relation between 20 glacier area and length. Glaciers in the South American Andes are found to be shorter for a given area compared to the other regions. The Patagonian Andes are dominated by ice fields at comparably low elevations with relatively short outlet glaciers. At low latitudes very steep glaciers prevail that are rarely organized as distinct valley glaciers and both regions are subject to rather steep balance gradients (Warren and Sugden, 25 1993;Benn et al., 2005)

Computation principle
Model development was guided by the idea to mimic glacier flow with simple and computationally efficient algorithms using the full information available in a DEM. The first condition of maximizing downhill slope imitates gravitational pull while the condition of 5 centrality emulates the guiding effect of the surrounding ice masses. Calculated flow trajectories are not to be mistaken for actual flow lines because glacier flow is a threedimensional phenomenon and cannot be derived from surface information alone. Furthermore, glacier flow lines adjoin on glacier tongues in parallel flow while our trajectories unite in one central line. Then again our approach does not necessarily generate 10 center lines in a strict sense. Adherence to the glacier center can be flexibly varied across the glacier perimeter. With the here applied settings, centrality is less rigid than in previously published methods (Le Bris and Paul, 2013;Kienholz et al., 2014). We have presented one possibility of implementing the basic conditions and designing trade-off functions; alternatives to improve accuracy and efficiency certainly exist.
For instance, our implementation is based on slope angles between individual grid cells whereas an approach involving averaging over perimeters of a few grid cells might be less sensitive to small scale topographic features. One might also imagine involving additional conditions in future updates, such as the direct inclusion of observed glacier flow fields. 20 The approach is computationally efficient; 41 000 km 2 of glaciers at 90 m spatial resolution (East Greenland) are calculated within 15-20 min on an ordinary laptop computer. This efficiency is used to apply a "brute-force" method of start point sampling. The point density is simply set high enough so that every glacier branch receives several starting points. The probability that the most distant point is among the sampled 25 ones is high, independent on whether the point is located on a summit, a pass or elsewhere. However, establishing a geometric order of center lines as done by Kienholz et al. (2014) might be more challenging given the very large number of center lines. 2508 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Our approach is here applied in a fully automated manner, but is programmed to involve similar possibilities of manual intervention as illustrated by Kienholz et al. (2014). In view of the worldwide calculation none of these options were used however.

Model performance
Model evaluation indicates good performance for the East Greenland test site while the 5 comparison for Alaska shows stronger deviations between the methods.
The major reason for the success of the East Greenland calculation is seen in the quality of the GIMP-DEM which allows for accurately calculated center lines. Basic characteristics of topography are well represented in the DEM and the existing spurious surface undulations are small enough not to interfere with the calculations. A fur-10 ther reason for the good agreement is that manual drawing of center lines was strongly oriented on surface topography and only obviously erroneous topography features were ignored. The fact that the model was calibrated for the region is likely of limited relevance because the area comprises virtually any possible glacier type including very complex glacier shapes. 15 The laborious manual measurement allowed us to measure only 100 glaciers randomly picked from ten predefined size classes. Thus small glaciers are underrepresented compared to the full set of East Greenland glaciers. Furthermore, the average glacier area in East Greenland is larger by almost a factor three compared to Alaska (11.0 km 2 vs. 4.0 km 2 ). The influence of ambiguities and issues related to small glaciers 20 is thus underrepresented in the model evaluation for Greenland. Nevertheless, a visual review of most Greenland center lines (c.f. Fig. 5) confirms the good performance seen in the 100 glaciers sample.
In the Alaska comparison a good correlation is found where DEM quality is comparable to the East Greenland GIMP-DEM but strong deviations in measured glacier length 25 occur on low quality DEMs. A qualitative assessment with focus on the two DEMs most frequently used in the global calculation, suggests that performance is good on the SRTM DEM but worst on the ASTER GDEM v2. On the latter center lines are 2509 Introduction frequently terminated midway where DEM errors suggest that a glacier tongue flows uphill over a distance longer than the applied search radius. The approach by Kienholz et al. (2014) mostly maintains the ability to continue center lines because the method does not suppress uphill-flow for the major (lower) part of a glacier. Low-quality DEMs are also the main reason for different starting points chosen by the two approaches.

5
Basically our high sampling density would guarantee that nearly everywhere starting points are picked in close vicinity of the points chosen by Kienholz et al. (2014). However, DEM irregularities might block progress towards the glacier tongue and force center lines to end too early. Other lines starting lower down might become longer and be eventually chosen to represent total glacier length. 10 On the example of Alaska we compared our fully automated method to semiautomatic measurements that involved a total of ∼ 4300 manual interventions (Kienholz et al., 2014). Only a rough estimate can be provided how our method would perform when a similar number of manual corrections would be applied. We assume that DEM related issues could be alleviated only marginally, but manual definitions of glacier 15 endpoints would be of considerable benefit on piedmont-type glacier tongues where the automatic selection of endpoints often delivers unsatisfying results.
The difficulties with low-quality DEM data root in the strong computational involvement of surface slope. Relying on the latter, however, is an advantage on certain glacier types and more generally when calculating center lines on high-quality DEMs. Consid-20 ering surface slope is of particular importance on wide slope glaciations, on asymmetric glacier shapes, on ice caps, and where broad accumulation areas narrow into tongues of outlet glaciers. When relying strongly on centrality, center lines can run almost parallel to elevation contours. For such specific glacier types, and more generally wherever high-quality DEMs exist, our approach allows more strictly controlling unphysical lateral 25 flow.
Our approach has no strong dependency on catchment delineations due to the strong involvement of surface slope. If the trade-off between centrality and slope is recalibrated accordingly, approximate flow lines can be calculated for glacier complexes or ice caps without any catchment delineations. On the contrary, the importance of centrality in the approach by Kienholz et al. (2014) leads to a stronger dependency on accurate catchment delineations which are only possible when a certain level of DEM quality is given. Thus, both approaches depend on DEM quality, although to a varying degree and at different stages of the calculations. 5

Worldwide glacier length
Glacier length is an important, yet missing parameter in global glacier inventories (Paul et al., 2009). Here we provide a first globally complete assessment of the length of all ∼ 200 000 individual glaciers around the globe. Based on our data we investigate the relationship of glacier length and glacier area and calculate global and regional scaling 10 laws. Differences between the regions appear to be related to regional characteristics of topography and mass balance distribution. Due to large variability in glacier shapes, the scaling laws allow only a rough estimate of the length of individual glaciers. A particularly wide spread exists among small glaciers because of ambiguities in defining their length (c.f. Le Bris and Paul, 2013), but given their small size the application of 15 scaling laws involves only limited absolute errors. For reasons discussed above, our global data set of glacier length is subject to regionally varying quality. Low uncertainties can be expected where the SRTM or better quality DEM data was used, i.e. in-between 60 • N and 55 • S, as well as for the Greenland (GIMP DEM) and the Antarctic (RAMP DEM) periphery. Limited accuracy 20 is anticipated for most Arctic regions where the ASTER GDEM v2 had to be used, in particular for Arctic Canada and the Russian Arctic. Glacier length data from areas with RGI quality limitations should also be used with care.

Conclusions
We have presented a fully automated method to calculate glacier center lines based on surface slope, distance to the glacier margin and a set of trade-off functions. The approach was developed on the example of East Greenland, evaluated for the same area as well as for Alaska, and eventually applied to obtain the first global assessment 5 of glacier length. By calculating glacier length for each of the ∼ 200 000 glaciers in the RGI we add an important and previously unavailable parameter to global glacier inventories.
Our scaling laws and the differences in scaling factors among regions can be applied and investigated in the framework of conceptual glacier models. Global glacier 10 length data could potentially be used to assess changes in length for different regions or glacier types. The actual center lines might be beneficial to flow-line modeling approaches.
Our approach calculates center lines by mimicking glacier flow, is computationally efficient and fully automated. Thus three of the initially stated four conditions are met. 15 While accurate center lines result when using good quality DEM data, further research is needed to reduce the methods's sensitivity to DEM inaccuracies. With the upcoming TanDEM-X data in view (e.g., Martone et al., 2013), however, we believe that our basic concept is well suited for future use. Once these precise worldwide terrain data are available, we will aim at providing a high-quality data set of global glacier length.

Appendix A A1 Trade-off functions
In the following the three trade-off functions are explained in full detail. Figure 2b shows a numerical example of the use of the trade-off functions.

2512
Since t 2 is set to 0.35 (Table A1) the example in Fig. 2b does not fulfill the condition |w u /w d − 1| > t 2 and c 2 remains zero.
T3: The trade-off function firstly calculates the mean slope angle above (α u ) and below (α d ) the current elevation of the flow path by averaging slope to all uphill and downhill grid cells. If the ratio α u /α d exceeds a certain threshold t 3 and the mean slope 15 (α u + α d )/2 is within a range of 0.02 to 0.1 then a factor c 3 is calculated according to c 3 = f 3a (α u + α d )/2. In case that (α u + α d )/2 > 0.1 then c 3 = f 3a · 0.1 In addition α d is checked and if exceeding a threshold t 4 then c 3 = c 3 +α d ·f 3b . In the example of Fig. 2b all three conditions are met: α u /α d > t 3a , 0.02 ≤ (α u + α d )/2 ≤ 0.1 and α d > t 3b . Thus c 3 is calculated to be 0.67.