Comparison of multiple glacier inventories with a new inventory derived from high-resolution ALOS imagery in the Bhutan Himalaya

Digital glacier inventories are invaluable data sets for revealing the characteristics of glacier distribution and for upscaling measurements from selected locations to entire mountain ranges. Here, we present a new inventory of Advanced Land Observing Satellite (ALOS) imagery and compare it with existing inventories for the Bhutan Himalaya. The new inventory contains 1583 glaciers (1487±235 km), thereof 219 debris-covered glaciers (951± 193 km) and 1364 debris-free glaciers (536±42 km). Moreover, we propose an index for quantifying consistency between two glacier outlines. Comparison of the overlap ratio demonstrates that the ALOS-derived glacier inventory contains delineation uncertainties of 10–20 % which depend on glacier size, that the shapes and geographical locations of glacier outlines derived from the fourth version of the Randolph Glacier Inventory have been improved in the fifth version, and that the latter is consistent with other inventories. In terms of whole glacier distribution, each data set is dominated by glaciers of 1.0–5.0 km area (31–34 % of the total area), situated at approximately 5400 m elevation (nearly 10 % in 100 m bin) with either north or south aspects (22 and 15 %). However, individual glacier outlines and their area exhibit clear differences among inventories. Furthermore, consistent separation of glaciers with inconspicuous termini remains difficult, which, in some cases, results in different values for glacier number. High-resolution imagery from Google Earth can be used to improve the interpretation of glacier outlines, particularly for debris-covered areas and steep adjacent slopes.


Introduction
To assess the effects of ongoing climate change on glaciers, a comprehensive programme of glacier monitoring is underway throughout high-mountain Asia (e.g.Cogley, 2016).However, the inherent inaccessibility of Himalayan glaciers hampers an adequate number of in situ measurements, meaning that glacier response to climate change can only be inferred from a small number of accessible glaciers that tend to be small and/or located at lower elevations (Bolch et al., 2012;Cogley, 2012).Moreover, spatially heterogeneous shrinkage of Himalayan glaciers (i.e.rapid thinning in south-eastern/humid regions and slow thinning in northwestern/arid regions) has been identified via in situ measurements (see Bolch et al., 2012;Yao et al., 2012), remote sensing (Kääb et al., 2012;Gardelle et al., 2013), and numerical modelling based on climatic data sets (Fujita and Nuimura, 2011).More intensive long-term monitoring of Himalayan glaciers is therefore of great importance and will help reduce uncertainty in projections of global sea level rise (Radić and Hock, 2011), regional water availability (Kaser et al., 2010;Immerzeel et al., 2010), and risk assessment of glacier-related hazards (Richardson and Reynolds, 2000).All glacier-specific measurements of change (e.g.elevation, volume/mass, length, velocity) require a well-defined boundary of each glacier entity.However, assessment of glacier changes in high-mountain Asia is hindered by the timeconsuming process of delineating debris-covered glaciers manually.Using a data set created by others for change assessment is actually not recommended without adequate quality control.
A glacier inventory containing delineated outlines is a fundamental asset for assessing glacier change in these remote and extensive mountain regions.Existing glacier inventories for the Himalayas include those of the International Centre for Integrated Mountain Development (ICIMOD; Mool et al., 2001), the project of Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM; Nuimura et al., 2015), and the Chinese Glacier Inventory (Guo et al., 2015) of the Chinese Academy of Science.Subsequently, parts of each have been included in a global glacier inventory compiled by the Randolph Consortium, termed the Randolph Glacier Inventory (RGI) (Pfeffer et al., 2014).Several studies have utilised these data sets to clarify recent glacier change (e.g.Vaughan et al., 2013;Bajracharya et al., 2014a, b), glacier lake evolution (Gardelle et al., 2011), and relationships between precipitation and the distribution of glacier elevations (Sakai et al., 2015).The RGI was designed for and applied in several large-to global-scale assessments and modelling studies in support of IPCC AR5 (Vaughan et al., 2013), including past and future global sea level rise (Gardner et al., 2013;Radić et al., 2014;Marzeion et al., 2012) and ice-thickness distribution (Huss and Farinotti, 2012).
To date, the preferred method of generating multiple glacier outlines has involved a combination of automatic extraction and manual correction from mid-resolution multispectral satellite imagery, such as Landsat-7 and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery.Glacier ice is extracted using a band ratio classification of visible and near-infrared bands (e.g.Racoviteanu et al., 2008;Paul and Andreassen, 2009;Bajracharya et al., 2014a, b).Watershed analysis, using flow directions derived from digital elevation models (DEMs), can be employed for drainage division (i.e.dividing an ice surface into multiple glaciers along the ridge line) (e.g.Bolch et al., 2010;Kienholz et al., 2015) and for consideration of debris-covered areas (e.g.Paul et al., 2004;Bajracharya et al., 2014a).Following these automatic processes, visual quality checking and manual correction are usually required; complete visual interpretation and manual delineation using Landsat imagery have only been performed for the GAM-DAM inventory (Nuimura et al., 2015).
In contrast, high-resolution satellite imagery such as Quickbird, Ikonos (Paul et al., 2013), Panchromatic Remotesensing Instrument for Stereo Mapping (PRISM) onboard Advanced Land Observing Satellite (ALOS; Narama et al., 2010;Nagai et al., 2013), and Corona (Bolch et al., 2008;Narama et al., 2010) enable more detailed scrutiny of terrestrial features with spatial resolutions of a half to a few metres.For accurate delineation, visual perception by an investigator is more important than physical values obtained from multispectral satellite observation.Although manual delineation does not always result in precise and repro-ducible glacier outlines, it is superior to automatic classification when many debris-covered and shaded-ice areas require interpretation (Paul et al., 2013).For example, Nagai et al. (2013) used manual delineation and ALOS PRISM imagery to distinguish between debris-covered glaciers and adjacent hill slopes in the Bhutan Himalaya, where automatic classification is complicated by complex mountain topography.
Glacier inventories are a basis with enormous potential for use at a variety of scales and accuracies, such as glacierdistribution assessments over entire mountain ranges, calculations of climate-driven runoff change in specific river basins, and decadal reporting of terminus retreat.However, it is not yet clear whether existing inventories yield the same or different results for these various applications.The aim of this study, therefore, is to identify consistencies and inconsistencies in glaciological variables (e.g.area, number, elevations, slope gradient, and aspect) among the different glacier inventories and to discuss the potential causes of any disparities, thereby helping to identify the optimum inventory for each respective purpose and helping to improve future work.

Study site
Our study focuses on glaciers located in the Bhutan Himalaya, an eastern part of the Himalayan mountain range.The study site comprises a rectangular area (89.10-91.90• E, 27.50-28.70• N; Fig. 1) in which the main east-west Himalayan range forms the international boundary between Bhutan and China.The highest peak in our study site reaches 7500 m above sea level (a.s.l.) and forms the headwaters of the Brahmaputra River.The climate is dominated by the Indian monsoon, resulting in greater precipitation on southern side (> 2000 mm yr −1 ) than on northern side slopes (< 1000 mm yr −1 ; e.g.Bhatt and Nakamura, 2005;Bookhagen et al., 2006;Houze et al., 2007).
Glaciological characteristics also exhibit marked disparities between southerly and northerly aspects, such as the rapid retreat of debris-free southern glaciers relative to northern glaciers (Karma et al., 2003).Indeed, the latest ICI-MOD inventory suggests that glaciers on southern flanks have lost 23.3 ± 0.9 % in total area over the last 3 decades (Bajracharya et al., 2014a).Because southern glaciers possess more extensive debris cover than their northern counterparts, owing to the greater severity of diurnal freeze-thaw cycles on south-facing slopes (Nagai et al., 2013), southfacing glaciers are characterised by steep headwalls, lowangle debris-covered areas, and slower flow rates.In contrast, north-flowing glaciers typically exhibit gentle slopes, fewer debris-covered areas, and higher flow rates (Kääb, 2005;Scherler et al., 2011a).

Data sets
We utilised ALOS satellite imagery.First, we used PRISMderived stereo-pair panchromatic imagery with a spatial resolution of 2.5 m, acquired between 2007 and 2011 (Table S1 in the Supplement), to generate DEMs.The DEMs were then used in conjunction with DSM (digital surface model) and Ortho-image Generation Software for ALOS PRISM (DOGS-AP; Takaku and Tadono, 2009; Fig. 2) to ortho-rectify the original images.This software was developed specifically for generating DEMs and ortho-rectified imagery from PRISM Level 1B1 standard imagery, which supports the rigorous sensor models of PRISM, ALOS orbit data, and ALOS attitude data (Takaku and Tadono, 2009).In this study site, we reported an accuracy of 0.37 m rootmean-square error for the horizontal component in positioning (Ukita et al., 2011).Because all PRISM images included cloud and seasonal snow cover, we selected the best multitemporal PRISM images available for each glacier, typically those collected between October and February.We selected 58 PRISM images from the 1364 available (Table S1).Multi-spectral visible and near-infrared imagery, with a spatial resolution of 10 m, was acquired via an additional ALOS sensor, the Advanced Visible and Near Infrared Radiometer Type 2 (AVNIR-2; Table S2).The AVNIR-2 images acquired between 2007 and 2011 were ortho-rectified using the 90 m Shuttle Radar Topography Mission (SRTM3 v3) DEM in DOGS-AP.Composite true-colour images were generated from bands 3, 2, and 1, corresponding to red, green, and blue, respectively.
We used the second version of the Global Digital Elevation Model by ASTER (ASTER GDEM2) to perform glacier delineations and spatial analysis.This model was provided by the Earth Remote Sensing Data Analysis Center, Japan Space Systems, with a cell size of 1 arcsec (∼ 30 m) and in geographic coordinates covering the global surface (Tachikawa et al., 2011).The absolute length of 1 arcsec varies depending on the latitude, therefore GDEM2 pixels were resampled to 30 m cell size on the Universal Transverse Mercator (UTM) coordinate system zone 46.We chose the ASTER GDEM2 over SRTM3 because the latter exhibited less consistency with Ice, Cloud, and land Elevation Satellite (ICEsat)-derived data in the high mountains of Asia, including a negative bias of −99 m and a large analytical uncertainty.A lower bias of +40 m was shown by ASTER GDEM2 (Nuimura et al., 2015).Precipitation data from the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) were used to explore the relationship between glacier distribution and climate potential for glacier accumulation.These data were generated using a calibrationbased sequential scheme that combines precipitation estimates from multiple satellites and rain gauge data at a pixel spacing of 0.25 • (∼ 25 km; Huffman et al., 2007).Yamamoto et al. (2011) used in situ measurements in the Nepal Himalaya to demonstrate that the TMPA data have superior consistency to other precipitation products.Nonetheless, the 0.25 • spacing means that both sides of a high mountain range are included within a single pixel and, therefore, that stark differences in precipitation (i.e.greater rain/snowfall on the windward, upstream side and less precipitation on the downstream side) are poorly resolved.For this reason, we avoided spatially detailed assessments in favour of providing an overview of the distribution of annual precipitation amounts.A monthly product, TRMM 3B43, was obtained for the period 1998-2014 and averaged to provide mean annual precipitation data.
In addition, we calculated the annual flux of solar radiation for each glacier.Intensity of received solar radiation depends primarily on slope inclination and aspect, and these data are provided via spatial analysis of the ASTER GDEM2 for each date and time.Using a basic tool in ArcGIS, solar radiation was estimated from the DEM (30 m cell size) every 0.5 h for each month, and those data were summed to provide an annual flux in GJ m −2 .This process considers the latitude of the study site homogeneously set at 28 • N, as well as the effects of shading by adjacent mountains, but does not consider actual cloud cover and other meteorological factors.Calculated values are averaged within individual ALOS-derived glacier outlines to produce a single value for each glacier.

Manual digitisation of glacier outlines
We designed a practical procedure for delineating glaciers and debris-covered areas in PRISM and other satellite imagery (Fig. 2).Because automatic classification of glacier surfaces from a multispectral image cannot yet be applied to such high-resolution panchromatic imagery (i.e.single-band image), time-consuming manual delineation is currently the only way to take advantage of the high-resolution images needed for identifying complex features on glacier surfaces.Additionally, whereas Nagai et al. (2013) delineated outlines for only debris-covered glaciers, this study makes delineations for all glaciers, including debris-free glaciers, and provides attribute information of glaciological parameters for each.Our Bhutanese study follows the definition of a glacier proposed by Global Land Ice Measurements from Space (GLIMS) (Raup and Khalsa, 2007), with the sole exception that we excluded adjacent snow-covered slopes on which glacier flow is not expected (e.g.Figs. 3, 4).
First, glacier outlines were generated as vector-format polygons traced onto the PRISM images, over which ASTER GDEM2-derived contour lines (20 m interval) were laid using ArcGIS.Poor visibility owing to shadows, which are particularly pronounced on the accumulation zones of northflowing glaciers in winter, was mitigated by changing brightness and contrast.Empirically, it is important to begin glacier delineation from the terminus in order to identify the total number of glaciers.Thus, in this study, we began glacier polygons at the terminus, formed a boundary between the glacier (including debris-covered areas) and surrounding icefree terrain, and terminated again at the terminus (Figs. 3,  4).Having delineated the outer line, internal bedrock and nunataks were removed from the polygon domain.Points forming glacier polygons were plotted automatically with 10 m intervals to form natural, realistic shapes.Additionally, all PRISM-derived outlines were checked through interpretation of composite visible colour in 28 AVNIR-2 scenes (Table S2) to find misidentification caused by PRISM's  3).(T) Heavily snow-covered surface, on which ice flow is expected, is identified as part of the glacier, whereas the regions (U) and (V) are snow-covered areas with poor snow depth for accumulation that are excluded from the glacier because surface features suggest no ice accumulation on these surfaces.The isolated ice mass (W) is included as part of the glacier because its mass probably provides ice to the glacier.panchromatic imagery.PRISM data were unavailable for 32 glaciers located in the northernmost part of the study site, for which AVNIR-2 imagery alone was used.Delineated outlines smaller than 0.01 km 2 in area (n = 20, with a total area of 0.15 km 2 ) were removed to minimise the risk of including snow patches.
Manual delineation of debris-covered areas requires careful observation and interpretation of surficial features in high-resolution images.Easily identifiable thermokarst-like features, such as small ponds, ice cliffs, and surface relief with numerous small rises and inclines, were distinguished in PRISM images from surrounding bedrock and moraine (Fig. 4).Additionally, in numerous instances, stream channels emanating from the margins of debris-covered termini revealed key differences in internal structure between supraglacier debris sediment and bedrock terrain, thereby aidwww.the-cryosphere.net/10/65/2016/The Cryosphere, 10, 65-85, 2016 ing the delineation of glacier outlines (e.g.Q in Fig. 4).Boundaries between debris-covered and debris-free ice surfaces (i.e. the onset of debris cover) can vary among PRISM scenes of different dates owing to snowfall.In such cases, the most expansive debris-covered area was adopted and delineated.Several rock glaciers, flowing masses of boulders with less ice content, are identifiable in the study site by their steep front and furrowing and creeping forms in high-resolution satellite imagery.In most cases, they do not have the thermokarst-like features found on a typical debriscovered glacier, though gradual transitions between debriscovered glaciers and rock glaciers can exist.This study does not delineate rock glaciers because they are not defined as a type of glacier (Raup and Khalsa, 2007).Snow-covered rock walls adjacent to the ice surface were also excluded from our delineations of glacier area (Fig. 4d).First, we focused on the difference in the concentration of contour lines between glacier surfaces and rock walls, the former being relatively gentle compared with the latter.This difference results in relatively sparsely distributed contour lines on the ice surface and more concentrated contour lines on the walls.The contact between the two represents the edge of the glacier, even if it is obscured by snow in the PRISM imagery.Partly snow-covered walls, with exposed and identifiable bedrock, were excluded.Careful division was required for those glaciers that are connected to deeply snow-covered rock walls lacking bedrock exposures.Slopes exhibiting unbroken snow cover and that appeared to be thick and stable ice masses, and which suggested regular flow rather than avalanche feeding, were included in glacier outlines (e.g.T in Fig. 4).However, the specific threshold value of a DEMderived slope gradient was not defined for this delineation since hanging glaciers can have steeper gradients than the surrounding bedrock and must be included.
Ice surfaces connected to multiple glacier termini were separated along a hydrological boundary, such as a watershed or ridge line, by delineating orthogonal lines against GDEM2-derived contour lines (Fig. 5).Unlike the second Chinese Glacier Inventory (CGI2) (Guo et al., 2015), we did not perform automatic watershed separation because we intended to compare the results of automatic separation against visual interpretation in complex terrain.Typically, the generation of DEMs for snow-covered surfaces is problematic, as control points can be obscured by snow.However, we did not observe any unnatural contour lines, which hamper manual delineation, at this site.Careful interpretation is required to distinguish connecting glacier termini (Fig. 5).In the case that separation was difficult, they were integrated into a single glacier.Similarly, any detached ice masses located just above a glacier were included in the delineation if ice redistribution was expected.Shapes for these features are recorded and registered in the inventory as a single row in the attribute table, with the same ID, area, and other variables.
Upon delineation of glacier outlines and debris-covered areas, outlines were corrected by review of exported glacier polygons into Google Earth, which includes high-resolution Quickbird images superimposed upon the SRTM3 topography (Fig. S1 in the Supplement).Detailed surface features and 3-D visualisation are invaluable for identifying delineation errors (e.g.Raup et al., 2014).Whenever a significant error was identified in cloud-free, high-quality Google Earth imagery, the original outline was corrected manually in ArcGIS (in another screen).Especially outlines of north-facing accumulation areas neighbouring steep mountain ridges are carefully checked because of mountain shadows in the PRISM and AVNIR-2 imagery acquired in the winter season.To validate the final corrected version, this manual delineation process was repeated four times.Moreover, in order to validate manual delineations based on ALOS imagery, each author independently delineated 10 debrisfree glaciers and six debris-covered glaciers, as well as their debris-covered areas.Four outlines were provided for each glacier entity (as well as one debris-covered area) and consistencies among them were evaluated using the overlap ratio.

Glacier attributes
Statistical information for delineated glaciers was calculated by ArcGIS tools and recorded in an attribute table, provided here as supplementary material.A local identification (ID) number was assigned for each glacier polygon.Spatially disconnected polygons, regarded as one glacier entity according to GLIMS guidelines (e.g.W in Fig. 4), were merged and assigned a single ID.Scene ID and acquisition date for the primary PRISM image and the most appropriate AVNIR-2 image were also noted in the attribute table.
We calculated the areas of debris-free and debris-covered glacier surface from their horizontally projected polygons.
Elevations of maximum, minimum, mean, and median elevation were obtained for individual glaciers using the ASTER GDEM2.Median elevation represents the altitude at which a glacier polygon is separated into upper and lower halves of equal area.This point constitutes an index for denoting glacier distribution relative to a balanced-budget equilibrium line altitude (e.g.Braithwaite and Raper, 2009;Sakai et al., 2015).
We divided the study site into 11 river basins (Fig. 1) that were generated automatically using the ASTER GDEM2 and ArcGIS watershed analysis tools.On the southern flanks, basin names follow those reported by Mool et al. (2001).On the north side of the range, rivers flowing to the Tibetan Plateau were named Nianchu, Langkamu, and Subansiri in reference to the attribute table of the ICIMOD inventory (Bajracharya and Shrestha, 2011).
Slope gradients were obtained for each glacier by averaging the surface gradients of the ASTER GDEM2 for each glacier.Values for upper and lower halves were calculated separately and divided by the median elevation, and are assumed to represent slope gradients for the accumulation and ablation zones, respectively.The glacier aspect was calculated as the mean value of surface aspect on a glacier surface.The aspect value of a GDEM2 pixel (0-360 • ) is converted to a unit vector (x axis for longitudinal range and y axis for latitudinal range).Finally, the averaged unit vector is converted back to an aspect value (0-360 • ).

Other glacier inventories
In the present study site, existing glacier outline data sets are available from GAMDAM, ICIMOD, and RGI, and were referenced onto the Universal Transverse Mercator (UTM) coordinate system (zone 46).Contained outlines smaller than 0.01 km 2 in area were not analysed in this study to minimise the risk of including snow patches.We did not include the GLIMS database because the latest version, released on 2 December 2014, includes considerable duplication on the northern side of the range and because no glaciers were delineated on the southern side.
Comparison between the GAMDAM and ALOS-derived inventories highlights the differences between the types of satellite imagery used, as similar manual delineations were performed by similar authors.The GAMDAM inventory contains more than 83 000 glaciers distributed throughout the high mountains of Asia (Nuimura et al., 2015), of which we examined those lying within our study site.These glaciers were delineated manually using two Landsat-7 scenes acquired on 28 December 2000 and 20 November 2001, along with SRTM-derived contour lines (20 m interval) for glacier separation and differentiation of adjacent headwalls.Delineations were then verified via high-resolution Google Earth imagery.While the methods and criteria of the GAMDAM inventory are similar to our own, the ALOS-derived outlines used in this study were delineated solely by the first author, who did not participate in the Bhutan domain of the GAM-DAM inventory.Therefore, we can evaluate consistency between the ALOS and GAMDAM outlines objectively.
The first ICIMOD glacier inventory, published in 2001, utilised Landsat-5, Indian Remote Sensing (IRS)-1D, and Satellite Pour l'Observation de la Terre (SPOT) imagery acquired in the 1990s, in addition to a 1 : 50 000 scale topographic map depicting glacier distribution in the 1950s (Mool et al., 2001).The latest ICIMOD glacier inventory, published in 2014, includes glacier outlines for the Bhutan Himalaya that were generated using Landsat-5/-7 imagery for the periods 1977-1978, 1990, 2000, and 2010 with semi-automatic classification (Bajracharya et al., 2014a).We utilised the 2010 outlines, in which glacier ice was classified by a threshold of normalised difference snow index (NDSI).Erroneously incorporated shadows, water bodies, vegetation, bare rock, and debris outside of glacier surfaces were removed by multiple filters including the normalised difference vegetation index, a land and water mask, and mean hue with respective threshold values (Bajracharya et al., 2014a).Slopes > 60 • and surfaces < 4600 m a.s.l. were removed using the SRTM3 DEM.Debris-covered areas in the elevation range 3000-6000 m a.s.l. were extracted from the remaining terrain via the SRTM3-derived slope gradient (< 25 • ), normalised difference vegetation index, NDSI, and land/water mask.Finally, outlines of debris-covered and debris-free glaciers were superimposed onto high-resolution satellite images in Google Earth and corrected manually (Bajracharya et al., 2014a).
The RGI has a global coverage, for which glacier outlines were imported and merged from various sources (Pfeffer et al., 2014).In producing the RGI, the original intention was to compile the best possible data available at that time in order to permit global-and regional-scale analyses for the IPCC report.For our study site, glacier outlines were extracted from part of Region 13 Central Asia and Region 15 South Asia East of version 4.0 (RGI4.0)and version 5.0 (RGI5.0).The RGI4.0 includes glacier outlines from the first Chinese Glacier Inventory, which was based on aerial photographs and topographic maps from 1980 (Shi et al., 2009), and from a previous version of the ICIMOD inventory (Mool et al., 2001).The RGI4.0 is now outdated, but was frequently used for earlier comparisons.To improve quality, the RGI5.0 contains glacier outlines from the CGI2 for the northern flanks and from the GAMDAM inventory for the southern side.The CGI2 includes glacier outlines derived from topographic maps published between 1966 and 1971 and in 1980, as well as from Landsat TM (Thematic Mapper) images acquired in 2007, 2009, and 2010(Guo et al., 2015)).Six glacier outlines were based on topographic maps of unknown age.Automatic delineation of band ratio segmentation and ice-divide extraction by SRTM3 v4 were performed, followed by manual editing.Those RGI inventories are thus combined data sets, however delineation origins are not noted on individual glacier outlines.

Overlap ratio
To evaluate the geometric consistency among the multiple inventories, we propose an index of overlapping glacier outline between two inventories (r ov ), which is defined as where C is the overlapping area, and S A and S B are the total areas of a glacier polygon in two inventories of A and B, respectively.A higher r ov value means a larger proportion of the overlapping area is shared by two inventories, and an r ov of zero signifies no overlap.For example, outlines from ALOS (1.24 km 2 ) and RGI5.0 (0.94 km 2 ) produce an overlapped area (0.85 km 2 ), whereas an outline without overlapping area causes zero of r ov (Fig. 6).If one polygon in the A inventory is overlapped by several polygons in the B inventory, the largest r ov is adopted as the ratio of the targeted glacier in the A inventory.
Between two polygons from different inventories for one glacier, relative omission and commission parts are incorporated into the value of overlap ratio.If one polygon (S A ) is smaller than another (S B ), the shared part (C) is also smaller than S B , which causes a small value for r ov even if C equals S A .This means that a large area of omission causes small r ov values.On the other hand, because this is a relative comparison of two delineated polygons from different sources, absolute omission/commission errors based on actual glacier entities cannot be described.Multiple combinations should be discussed to evaluate consistencies of each delineation.Thus, the overlap ratio enables assessment of the consistency among glacier outlines, including location shifting in various sizes, which would be difficult to perform when using with absolute value of delineated areas.
To assess the consistency between two inventories, the mean value of the overlapping ratio (R B A ), in which the A and B inventories are defined as base and target inventories, respectively, is calculated as where N A is the number of glaciers in the A inventory.This value can vary according to the number of non-overlapping polygons (r ov = 0) in the two respective inventories, for example.

Validation of the ALOS glacier outline
Ten debris-free glaciers, six debris-covered glaciers, and their debris-covered areas were delineated from the ALOS PRISM/AVNIR-2 images by four authors operating independently.Overlap ratios were calculated for all combinations of the four authors (six ways).Mean values and standard deviations of the ratios calculated for the 16 glaciers are summarised against their ALOS-derived sizes (Fig. 7a).All debris-free glaciers give mean values for r ov greater than 0.7, with six exceeding 0.9.Debris-covered glaciers give values greater than 0.5, with larger glaciers giving higher (∼ 0.8) values.Similarly, debris-covered areas give higher values with increasing area.
To assess delineation uncertainties quantitatively, we compared and analysed the glacier areas of the four interpreters in closer detail.Mean values and standard deviations of the delineated area were calculated for each glacier to obtain their respective uncertainty ratios (i.e. standard deviation divided by mean area as a percentage) (Table 1).The uncertainty ratio varies from 1.97 to 48.62 %, with larger values implying greater uncertainty.Plotting the ratios against glacier size, uncertainties tend to be smaller for larger glaciers regardless of debris cover (Fig. 7b).The highest delineation uncertainty (as much as 50 %) is expected for the smallest (0.01-0.1 km 2 area) debris-covered glaciers and the lowest (< 10 %) for larger (1.0-10 km 2 ), debris-free glaciers.In a case of relatively larger debris-covered glacier (Fig. 7c), uncertainty ratio is indeed smaller (< 20 %), however absolute area of delineated glacier outlines greatly ranges from 42.51 to 63.78 km 2 (Table 1, Fig. 7c).Complicated topography with shadows and ice distribution at the surrounding steep slopes has the potential to cause large area differences with inconsistent interpretations.

Statistics of Bhutanese glaciers obtained from the ALOS glacier inventory
The ALOS-derived glacier inventory consists of 1583 glaciers, with a total area of 1487 km 2 and mean glacier area of 0.9 km 2 , and includes 219 debris-covered and 1364 debris- free glaciers (Table 2, Fig. 1).The area of the debris-free glaciers ranges from 0.01 to 5.9 km 2 (mean area 0.4 km 2 ), resulting in a total area of 536 km 2 , whereas that of debriscovered glaciers ranges from 0.06 to 79.6 km 2 (mean area 4.3 km 2 ), giving a total area of 951 km 2 (Fig. 8).The debriscovered area ranges from 0.01 to 21.2 km 2 , resulting in a mean area of 1.0 km 2 .The majority (n = 649) of ALOSderived glaciers belong to the 0.1-0.5 km 2 size category, whereas the greatest number of debris-covered glaciers (n = 102) belong to the 1-5 km 2 category.This category also contains the greatest glacier area (475 km 2 ), whereas the 10-50 km 2 size category contains the maximum area of debris- covered glacier (375 km 2 ).In larger size classes, debriscovered glaciers are both more numerous and more extensive than debris-free ice, representing 100 % of glaciers > 10 km 2 in area.
The anticipated uncertainty ratio for the 16 glaciers (Fig. 7b) was applied to the entire ALOS-derived glacier inventory.Uncertainties for the individual ALOS-derived glacier outlines are calculated by means of the trend lines shown in Fig. 7b.Summation of anticipated uncertainty reaches 235 km 2 (16 % of the total glacier area), where debris-covered and debris-free glaciers have anticipated uncertainties of 193 km 2 (20 %) and 42 km 2 (8 %), respectively (Table 2).In addition, they are classified into eight classes of size, where total area uncertainty is divided by total glacier area (Fig. 8c).Debris-covered glaciers in the smallest size class (0.05-0.1 km 2 ) exhibit the largest uncertainty ratio, at 52 %, whereas those in the largest size class (> 50 km 2 ) reveal the smallest uncertainty ratio, at 13 %.Similarly, debris-free glaciers exhibit the smaller uncertainty ratio for larger size classes (i.e.11 % for the 0.01-0.05km 2 class, 6 % for the 5-10 km 2 class).The area propor-tion of debris-covered glaciers versus debris-free glaciers is smaller/larger in smaller/larger size classes, with which influencing uncertainty ratio also shifts from ∼ 10 to ∼ 25 % for larger sizes.Integrating both types of glacier, therefore, the smallest uncertainty ratio is exhibited by the 0.1-0.5 km 2 class (10 %).The highest is exhibited by the 5-10 km 2 class (20 %) with which area proportion of debris-covered versus debris-free glacier is highest (221 to 12 km 2 ) (Fig. 8c).In the size classes larger than 10 km 2 , the uncertainty ratio is directly controlled by that of debris-covered glaciers without any debris-free glaciers.The maximum uncertain area is exhibited by the 1-5 km 2 size class (78 km 2 ) with the largest total area of debris-covered/debris-free glaciers.
Surface elevations of debris-covered and debris-free glaciers and of debris-covered areas were obtained from the ASTER GDEM2 and summarised as hypsographic curves with 100 m intervals (Fig. 9).Debris-covered glaciers range from 4000 to 7500 m a.s.l., with a maximum area of 77 km 2 at 5500 m a.s.l.Meanwhile, the debris-covered area ranges from 4000 to 6300 m a.s.l., exhibiting a maximum area of 23 km 2 at 4900 m a.s.l.Debris-free glaciers range from  4600 to 7000 m a.s.l., with a maximum area of 77 km 2 at 5400 m a.s.l.Although both debris-covered and debrisfree glaciers exhibit maximum areas at approximately the same elevation, the distribution range of the debris-covered glaciers is significantly wider (3500 m) with a larger total area (951 km 2 in Table 2) than that of debris-free glaciers (2400 m/536 km 2 ).
We compared the maximum, median, and minimum elevations of individual glaciers with their respective areas on a log scale (Fig. 10).The maximum and minimum elevations tend to be higher and lower, respectively, for larger glaciers, whereas median elevations show neither increasing nor decreasing trends.To assess elevation distribution in greater detail, we categorised glaciers by location.Northern basins include Nianchu, Langkamu, and Subansiri, while southern basins include Chamkhar, Mangde, Pho, Mo, Wang, and Amo.Traversal basins include Dangme and Kuri.An elevational comparison of glaciers of equal size revealed that the majority of glaciers located in southern basins occur at lower elevations than those in northern basins, while glacier elevation is widely distributed in the traversal basins (Fig. 11).
The spatial distribution of mean annual precipitation amount, obtained from TRMM 3B43 for the period 1998-2014, is shown in  the Tibetan Plateau.In contrast, southern regions experience in excess of 2000 mm precipitation per year, but glaciers do not exist there because the elevation is lower.The main Himalayan crest is characterised by west to east (W-E) oriented watersheds, occupied by numerous glaciers, and serves to hinder the northward passage of humid monsoon air masses; yet stark differences in precipitation are not apparent on either side of the range, as high peaks and valley bottoms are frequently included in the same pixel.Thus, even significant differences in local precipitation are not resolved fully by the 3B43 data.
The median elevations of individual glaciers were overlain onto a TRMM-derived spatial distribution of mean annual precipitation (Fig. 11), revealing a clear north-south contrast (i.e. higher values to the north and lower values to the south).Additionally, the spatial distribution of annual precipitation exhibits a north-south gradient, ranging from 462 to 1040 mm.Statistically, however, significant correlation was not observed for each glacier between annual precipitation and median elevation (r < 0.1), possibly because most parts of the significant north-south gradient of median elevation are distributed within a small area beyond mountain ranges, which is contained by one grid size of precipitation data (0.25 • ).
Mean slope gradients were calculated for both the upper and lower halves of individual ALOS-derived glacier surfaces, as well as for the entire surface (Fig. 12).The upper (lower) halves are assumed to comprise accumulation (ablation) zones.For the upper parts, the mean slope gradients reveal no correlations versus glaciers, whereas those for the lower parts denote negative correlations.This is related to the fact that a larger glacier, in many cases with a debris-covered ablation zone, tends to be developed on valley topography with a gentle slope (e.g.Scherler et al., 2011a).We observed no tendency in slope steepness between the differently categorised basins.
ALOS-derived glacier outlines were classified by aspect (Fig. 13).The number of debris-free glaciers on northfacing slopes is 3 times greater than on south-facing slopes, whereas debris-covered glaciers show little or no aspect dependency (Fig. 13a).Conversely, the area of the debriscovered glaciers exhibits considerable bias in both north and south basins (Fig. 13b).Mean values of median elevation for the eight aspects were calculated for each river basin (Fig. 14a).The highest ALOS-derived values occur in the Nianchu and Langkamu basins, followed by the Dangme and Kuri basins.Lower median elevations were observed in the southern basins.While south-facing glaciers have the highest median elevation in the Langkamu, Mo, and Wang river basins, in other regions we found no consistent inclination in the relationship between median elevation and surface aspect.
Total annual solar radiation was estimated using Ar-cGIS (Fig. 14b).We then calculated annual solar radiation (GJ m −2 ) at the ALOS-derived glacier surface using an Table 3. Surface area and overlap ratios for the whole glacier surface in the Bhutan Himalaya among the six inventories, in which one inventory represents one glacier surface.The numbers in italic font denote overlapped area (km 2 ).ASTER GDEM2 with 30 m cell size.Although the influence of surrounding topography was considered, the effects of cloud cover were impossible to determine.Results range from 4.0 to 9.6 GJ m −2 .The mean annual radiation flux for each glacier was then calculated and plotted against aspect.Our data show that north-facing glaciers received less radiation than south-facing glaciers.Additionally, the radiation flux over south-facing glaciers is apparently independent of slope gradient.

Consistency among glacier inventories
We first calculated overlap ratios for each of the five glacier inventories.The glacier surface was integrated such that one glacier area was determined for each inventory (Table 3).High ratios of > 0.85 were obtained for the ALOS, GAMDAM, ICIMOD, and RGI5.0 inventories, whereas the RGI4.0 inventory gave low ratios of < 0.75.Total overlap ratios for RGI4.0/5.0 suggest that while the RGI4.0 inventory contains larger areas that do not overlap with the other inventories, this situation has been improved in the RGI5.0 version.
Second, the number of glaciers is cumulated from those with a lower overlap ratio and cumulative curves for every combination of the inventories and is illustrated in Fig. 15.In comparison to the ALOS-derived inventory, the GAMDAM, ICIMOD, and RGI5.0 inventories reveal similar curves that increase towards larger r ov values (∼ 0.7), indicating that these data sets exhibit a high degree of well-matched outlines for ALOS-derived results with r ov > 0.7.In contrast, the RGI4.0 inventory exhibits constantly increasing curves, www.the-cryosphere.net/10/65/2016/The Cryosphere, 10, 65-85, 2016 on which the number of glaciers with low r ov (0.5) reaches 1000.Consequently, mean values for overlap ratios (R target base ) are highest for the ALOS-ICIMOD combination (0.63) and lowest for the ALOS-RGI4.0combination (0.35).Values for the ALOS-RGI combination show a marked improvement from RGI4.0 (0.35) to RGI5.0 (0.56).In analysing other combinations (Fig. 15b-e), the RGI4.0 inventory exhibits uniformly poor overlap with the other inventories.Figure 16a depicts area differences of the four inventories against our ALOS-derived inventory, in which the area difference among glacier outlines ranges from 0.1 to 10 times.The RGI4.0 in particular includes numerous glaciers for which the area exceeded 10 times that of the ALOS-derived inventory.Figure 16b shows variations in the standard deviation of area difference (ratios in Fig. 16a) calculated within six size classes.Smaller glaciers exhibit larger standard deviations, whereas larger glaciers show lower values.While each inventory shows this trend, the greatest standard deviations occur in the RGI4.0 inventory.

Glaciological variables among the inventories
Our inter-inventory comparison reveals that the RGI5.0 contains the lowest number (1266) of glacier outlines and that ICIMOD contains the highest number (1762).The minimum total area occurs in the GAMDAM inventory (1453 km 2 ), whereas the RGI4.0 exhibits the greatest, at 2261 km 2 (Table 2).We categorised these results into eight size classes (Fig. 8).The total glacier number and area differed widely among the inventories.However, in each, the majority of glaciers fall between 0.1 and 0.5 km 2 total glacier number (Fig. 8a).Similarly, the greatest total area lies in the 1-5 km 2 size class for each inventory, although the overall difference between RGI4.0 and GAMDAM reaches ∼ 300 km 3 in the size class (Fig. 8b).
We also compared the hypsometric curves of glacier outlines (Fig. 9).The maximum area of 187 km 2 occurs at 5600 ± 50 m elevation in the RGI4.0.The other inventories show a similar pattern.The ALOS, GAMDAM, ICIMOD, Glacier number increases from the lower r ov to the higher r ov , where a constant increase means a large number of low-r ov outlines, and a later high rate of increase means a large number of high-r ov outlines.Mean r ov values for each combination are denoted in the panels.and RGI5.0 outlines are located within the 4000-7500 m elevation range, whereas the lowermost outlines in the RGI4.0 occur at 2600 m a.s.l.
For each of the five inventories, we compared the maximum, median, and minimum elevations with their respective glacier areas (Fig. 10).Generally, a similar trend was observed: larger glaciers tend to have higher maximum and lower minimum elevations than smaller glaciers, although the median elevations uniformly cluster around 5500 m a.s.l.The slope gradients averaged over the upper half, entire surface, and lower half are depicted in on larger glaciers.Finally, we compared ALOS-derived total glacier numbers and areas in the eight aspect groups with the other inventories (Fig. 13a, b).Although we observed significant differences among the values, all five inventories exhibit similar trends.

Discussion
Our results suggest that although existing glacier inventories give different values for number, area, and topographic variables, they reveal similar statistical and topographical tendencies.The geometric differences in glacier outlines can be quantified using the overlap ratio, which enables relative evaluation among the inventories.In this section, we discuss glacier distributions in the Bhutan Himalaya obtained via the various inventories, examine the causes and effects of outline diversity using examples, and provide suggestions for making glacier delineations in areas of complex topography.

Glacier distribution in the Bhutan Himalaya
Each inventory indicates that glaciers with areas of 1-5 km 2 occupy the largest area, while those of 0.1-0.5 km 2 constitute the majority of glaciers by number (Fig. 8).Globally, the RGI5.0 suggests that this pattern is similar throughout the low-mid-latitudes, with the exception of the European Alps (i.e. a large number of small glaciers which do not affect areal proportion), whereas high-latitude regions are dominated by larger glaciers (Fig. S3).In the Bhutan Himalaya, glaciers exhibit a size distribution typical of low-latitude regions.
In this region, we observed that glaciers in southern basins have lower maximum, median, and minimum elevations than those in northern basins (Figs. 10,11).Also, we note that peak elevations are lower in southern than in northern basins (Fig. 9b in Nagai et al., 2013).Thus, we consider this topographic contrast to be a fundamental factor influencing the differences in maximum elevation between northern and southern basins (Fig. 10a).The north-south gradients in median elevation and precipitation in Bhutan show reasonable agreement, whereas north-south air temperature gradients are relatively minor due to the small latitudinal range (< 170 km) (Fig. 11).Therefore, we do not correlate the distribution of median elevations with glacier size, aspect, or debris cover, but to basin location, over which precipitation values differ significantly (Figs. 10b,11,14a) as found for the whole glaciers in high-mountain Asia (Sakai et al., 2015) and some glaciers chosen globally (Ohmura et al., 1992).Combinations of topography and precipitation may result in lower minimum elevations on southern glaciers, for which rapid retreat has been reported, and vice versa (Karma et al., 2003;Bajracharya et al., 2014a).
ArcGIS-based calculations of annual solar radiation give values of ∼ 9 GJ m −2 for south-facing glaciers regardless of surface gradient, whereas north-facing glaciers experience values of ∼ 8 GJ m −2 for gentle slopes and 5.5 GJ m −2 for steeper slopes (Fig. 14b).North-facing glaciers are more numerous than south-facing glaciers; this suggests a negative impact of solar radiation on small, debris-free glaciers on south-facing slopes.In contrast, large debris-covered glaciers occur on both southerly and northerly aspects (Figs. 1, 13b).The west to east orientation of the Himalayan range provides extensive north-and south-facing slopes (Fig. 1) favourable for the development of large valley glaciers.Potentially, debris cover is supplied to south-and north-flowing glaciers, although its intensity is greater on the south-facing slopes (Nagai et al., 2013).Total glacier area is dominated by debriscovered glaciers and is rather concentrated in the two aspects of the south and north, where the Himalayan main range forms larger valleys to develop large glaciers (Fig. 13).
In other regions, direction of the solar radiation and slope aspect of glacier distribution are not always correlated.In western China, similar areal concentration to the north is found for every subregion (Guo et al., 2015).In Tyrol, Austria, number and areal concentration are found in the north (Paul et al., 2011).However, in the Tian Shan Mountains, such concentration to the north is not found (Narama, 2010).In Alaska and north-west Canada, total glaciers denote homogeneous aspect distribution, whereas those in the northernmost part do areal concentration at the north (Kienholz et al., 2015).These variations suggest that glacier aspect distribution can be affected by local geography such as mountain range orientation, wind direction, precipitation homogeneity, etc. Therefore we cannot conclude that solar radiation from the south side is the only factor for the situation that the number of glaciers is smaller in the south-facing slopes, although it is a reasonable situation for the global trend that southfacing glaciers have higher elevation distribution than northfacing glaciers (Evans and Cox, 2005).

Inconsistency in glacier identification
Identification of glacier termini affects the consistency among glacier inventories.The largest discrepancy in the number of non-overlapping glaciers (n = 471) occurs between the RGI4.0 and GAMDAM inventory, whereas the smallest (n = 34) is found between the RGI5.0 and GAM-DAM inventory (Fig. 15).In some cases, small glacierlike polygons larger than 0.01 km 2 were found where no glaciers exist (Fig. 17a).Since glaciers are unlikely to disappear from one inventory to the next, we suggest that temporary snow cover is the main reason for this variability and, thus, that multiple images acquired on different dates should be checked.Additionally, snow-covered bedrock surfaces can present a similar appearance to glacier ice (e.g.Fig. 17b).However, while details of surface roughness in high-resolution Google Earth imagery suggest that these are not ice surfaces, such identification is not possible using the coarse resolution of Landsat-5/-7 imagery (30 and 15 m, respectively).

Inconsistency in glacier distribution
The entire RGI4.0 glacier surface overlaps poorly with those of other inventories (Table 3).Specifically, while the elevation of maximum glacier area lies around 5400 m a.s.l., similar to other inventories (Fig. 9), minimum glacier elevation in the RGI4.0 is considerably lower at 3000 m a.s.l.Individual glacier outlines in the RGI4.0 also overlap poorly with those in the other target inventories, for which mean values are no greater than 0.35 (Fig. 15d).We found significant distortion in RGI4.0 outlines when compared with other data sets (Fig. 18), which implies inaccurate cartography and low-quality georeferencing of the original topographic maps as well as a poor condition of the original imagery.This is one possible reason for the uniquely larger glacier area reaching 300 km 3 larger than that of the GAMDAM inventory (Fig. 8b).In contrast, RGI5.0 does not exhibit these inconsistencies, nor do the ICIMOD, GAMDAM, and ALOSderived inventories, all of which contain large numbers of well-overlapped glacier outlines with values > 0.5 (Fig. 15).The ICIMOD inventory was generated using semi-automatic processes, whereas GAMDAM and ALOS-derived outlines were generated via manual delineation.In all cases, the use of Google Earth imagery for manual correction is likely to result in outlines of similar quality.The RGI5.0 contains GAMDAM and CGI2 outlines for southern and northern sides, respectively.Figure 15e shows that the RGI5.0 comprises a large number of delineations that are highly consis- tent with GAMDAM data, of which 344 outlines of 1266 exhibit a perfect match.Importing GAMDAM outlines to the RGI5.0 inventory thus contributes to the consistency with ICIMOD and ALOS-derived inventories (Fig. 15d, e).
Excessive under-or overestimation of individual glacier areas produced results that range between factors > 0.1 and 10 of the ALOS-derived outlines, although the majority of results plot along the one-time line (Fig. 16a).This large difference might be due to the number of glaciers identified from one connecting glacier surface, rather than inaccurate delineation or poor georeferencing.While automatic GIS (geographic information system)-based separation of glaciers along ridge lines is necessary to resolve this issue, a remaining issue is indistinct glacier termini, which can reduce the consistency of glacier counting (Fig. 5).
High inter-inventory consistency at the minimum (i.e.terminus) elevation is requisite if glacier retreat is to be analysed using multitemporal data.In addition to the comparatively stable termini of debris-covered Himalayan glaciers (Scherler et al., 2011b), misidentification of the terminus can also affect interpretation of retreat rates for debris-covered glaciers (e.g.Fig. 4b).High-resolution Google Earth satellite imagery reveals numerous ice cliffs on a debris-covered area, suggesting a hidden ice body.Such details are difficult to identify in 30 m resolution Landsat imagery.Detailed visual inspection and manual correction with high-resolution satellite images is thus needed, particularly for smaller debriscovered glaciers, which, as described above, exhibit greater inconsistency in manual delineation than larger glaciers do (Fig. 7b).
Steep slopes surrounding the accumulation zone are excluded by ALOS-derived, GAMDAM, and ICIMOD inventories, but are included in the RGI4.0 and CGI2.Such headwalls constitute an additional source of glacier nourishment via avalanching (e.g.Benn and Lehmkuhl, 2000;Scherler et al., 2011a;Hewitt, 2011), which must be considered in studies of glacier mass balance and runoff in high relief regions.For example, glacier nourishment from surrounding slopes is especially relevant to the behaviour of Karakoram-type (i.e.Turkistan and Mustagh glaciers) (Hewitt, 2011).In contrast, headwall areas should be avoided when monitoring changes in glacier mass.Thus, the choice of inventory must be based on the purpose of study.
Finally, despite following the same definition of a glacier, differences in interpretation among investigators can result in considerable variability in glacier outlines where steep snowcovered slopes occur immediately adjacent to the glacier (Fig. 7c, e).This factor should be considered during the manual-editing process of automatic delineation by making a comparison with Google Earth imagery.Nonetheless, con-tinual updating of satellite imagery, temporally variable snow cover, and coarse resolution of images preclude straightforward glacier mapping.

Conclusion
We compared currently available glacier inventories for the Bhutan Himalaya with a newly created ALOS-derived glacier inventory.After manual delineation and quality control of glacier outlines, the ALOS-derived inventory contains 1583 glaciers with the total area of 1487 km 2 .Delineation uncertainty ratio is estimated as 10-20 % in area.We propose an index of overlap ratio to quantify the consistency among glacier inventories, which can be used to improve glacier delineation and compare glacier outlines from different sources.In terms of glacier outline, the GAMDAM and ICIMOD inventories were highly consistent with ALOSderived data, despite being generated from different sources (Landsat/ALOS), by different methods (automatic/manual), and by different agencies (ICIMOD/Nagoya University).Updating of the RGI4.0 has improved the consistency of the RGI5.0 with other inventories.Specifically, the historic data set from the first Chinese glacier inventory was replaced by newer data, in addition to incorporation of the GAM-DAM inventory.While specific values for glacier area, elevation, and slope gradient differ among the inventories, we observed broad similarities in common area classes, such as the largest number of glaciers (0.1-0.5 km 2 ) and area covered (1.0-5.0 km 2 ), elevation interval of the maximum area (∼ 5400 m a.s.l.), correspondence of hypsometry (Fig. 9) and surface gradient (Fig. 14) to their respective horizontal areas, and aspect distribution (Fig. 13).Although in situ measurements are difficult for the majority of glaciers in high mountain regions, digitised glacier outlines can be validated by comparing those from different sources.
The selection of a specific inventory depends on the purpose of an application.While any of them is appropriate when overview information only is required, the correct data set needs to be selected carefully if specific glaciological variables are of interest and individual glaciers are to be analysed.For example, glacier area can vary significantly owing to difficulties in the interpretation of debris-covered areas (Fig. 4) or the separation of individual glaciers (Fig. 5).
The Supplement related to this article is available online at doi:10.5194/tc-10-65-2016-supplement.

Figure 1 .
Figure 1.Glacier distribution of the ALOS-derived glacier inventory in the Bhutan Himalaya, in which glaciers are categorised into debriscovered and debris-free glaciers.Eleven river basins are categorised into (N) northern, (T) traversal, and (S) southern basins.The background is a hillshade image derived from the ASTER GDEM2.

Figure 2 .
Figure 2. Flow chart of data processing and manual delineation for the ALOS-derived glacier inventory.Control and tracing of delineation quality were conducted four times prior to finalising glacier outlines.

Figure 3 .
Figure 3. Examples of manually delineated outlines of glaciers and debris-covered areas.

Figure 4 .
Figure 4. Complex surface features on the termini of debris-covered glaciers in our study site, shown in (a) ALOS PRISM imagery, (b) Landsat-8 imagery, and (c) Google Earth (location shown in Fig. 3).(P) Small ponds, ice cliffs, and bumpy surface relief, as well as (Q) traces of water stream, are recognised in (a) and (b), whereas (R) bedrock and moraine surfaces are hard to identify in (b).(S) Debriscovered ice is partly exposed and identified.Delineation of a glacier with steep surrounding slopes located in our study site, shown in (d) ALOS PRISM, (e) Landsat 8 imagery, and (f) Google Earth (location shown in Fig.3).(T) Heavily snow-covered surface, on which ice flow is expected, is identified as part of the glacier, whereas the regions (U) and (V) are snow-covered areas with poor snow depth for accumulation that are excluded from the glacier because surface features suggest no ice accumulation on these surfaces.The isolated ice mass (W) is included as part of the glacier because its mass probably provides ice to the glacier.

Figure 5 .
Figure 5. Separation of connected glaciers in our study site.Basic separation is performed along ridge lines identified by contour lines.Glacier surfaces with vague termini were not separated.

Figure 6 .
Figure 6.Conceptual example of the overlap ratio (r ov ).The background is an ALOS AVNIR-2 image.

Figure 7 .
Figure 7. Validation of manual glacier delineation using ALOS imagery.(a) Mean values of r ov and their standard deviations, from six combinations among four outlines, are plotted against glacier area (or debris-covered area).(b) Deviations of glacier delineation in 16 sampled glaciers.(c) Debris-covered glaciers and debris-covered areas contain inconsistencies among four interpreters (coloured respectively), whereas (d) a simple shape of a debris-free glacier exhibits good consistency.(e) A debris-free glacier with steep slopes has inconsistencies on the higher parts of the outline (region in shadow).

Figure 8 .
Figure 8. Cumulative glacier (a) number and (b) area values for eight size classes.(c) Estimated area of delineation uncertainty ratio is shown on debris-covered, debris-free, and total glaciers from the ALOS-derived glacier inventory.

Figure 9 .
Figure 9. Hypsographic curves of glacier surfaces from the ALOS, GAMDAM, ICIMOD, RGI4.0, and RGI5.0 inventories.Also shown are the curves for the debris-covered area, and the debris-covered and debris-free glaciers, extracted from the ALOS inventory.
Fig. 11.Maximum precipitation (1136 mm) occurred at 27.625 • N, 91.375 • E and minimum precipitation (417 mm) at 28.375 • N, 89.625 • E, resulting in a mean value of 807±202 mm for our study site.Northern areas received lower amounts for each degree of longitude, and vice versa.The most arid basin is the Nianchu river basin on

Figure 10 .
Figure 10.Glacier (a) maximum, (b) median, and (c) minimum elevations in the ALOS inventory plotted against the log-scaled area.Glaciers are categorised as southern, northern, and traversal basins, as described in Fig. 1.Mean (lines) with standard deviation (vertical bars) of all inventories were calculated in four size classes.

Figure 11 .
Figure 11.Spatial distributions of median elevation for the ALOS-derived glacier outlines, overlain on the spatial distribution of annual precipitation from TRMM 3B43 data (1998-2014).

Figure 12 .
Figure 12.Slope gradients of the (a) upper half, (b) entire surface, and (c) lower half of ALOS-derived glacier outlines, categorised as northern/southern/traversal basins as shown in Fig. 1.Calculated mean values and standard deviations for four size classes are shown along with those of other inventories.

Figure 13 .
Figure 13.Cumulative glacier (a) number and (b) area of ALOSderived glacier outlines, classified by aspect.Also plotted are total glacier number and area data from other inventories.

Figure 14 .
Figure 14.(a) Median elevations in 9 of 11 river basins (Fig. 1) are plotted along eight aspects.(b) Estimated mean values of annual solar radiation for respective ALOS-derived glacier outlines, for which glaciers are grouped into three classes based on their mean slope gradients.

Figure 15 .
Figure 15.Cumulative glacier numbers against overlap ratios (r ov ) for the six glacier inventories covering the Bhutan Himalaya.Glacier number increases from the lower r ov to the higher r ov , where a constant increase means a large number of low-r ov outlines, and a later high rate of increase means a large number of high-r ov outlines.Mean r ov values for each combination are denoted in the panels.
Figure 16.(a) Area difference between the ALOS-derived glacier inventory and other inventories.(b) Standard deviations of area difference are summarised in seven area classes.

Figure 17 .
Figure 17.(a) Glacier outlines in part of our study site.Some outlines are delineated only in the (P) ICICMOD, (Q) RGI4.0, and (R) GAMDAM inventories.(b) A close-up image from Google Earth suggests that no glacier exists on the (S) snow-covered hill slope (location shown in box in a).(T) We found differences among the inventories in their separation of surrounding steep slopes.

Figure 18 .
Figure 18.Glacier outlines containing significant distortion in the RGI4.0.These have been improved for the RGI5.0 and are now consistent with other inventories.