A monitoring system for mountain glaciers and ice caps using 30 meter resolution satellite data

We developed a monitoring system for deriving outlines of mountain glaciers and ice caps (MG&IC) at a 30 m horizontal resolution from Landsat Thematic Mapper (TM) and Landsat Enhanced Thematic Mapper plus (ETM+). Location and area information at 30 m resolution was obtained using a band ratio (TM4/TM5) and a threshold value of TM3 with a 9 by 9 pixel average filter. The total area and number of MG&IC were 449482 km and 414258, respectively. The glacier outlines were similar to previous satellite-derived products for different regions. Although the derived glacier area was similar to previous estimates at regional scales, it was overestimated in some parts of Scandinavia where available satellite images are limited and only snowy season images can be used, and was underestimated in the western Himalayas and Caucasus where the glacier outlines are derived with difficulty from satellite images because of the effect of debris cover. Our system to monitor MG&IC has potential application in global hydrological and land-surface models and estimates of global sea-level rise.


INTRODUCTION
Loss of volume in mountain glaciers and ice caps (MG&IC) has been reported around the world in response to atmospheric warming (IPCC, 2007), significantly affecting terrestrial water resources (Immerzeel et al., 2010;Kaser et al., 2010) and contributing to global sea-level rise (Rignot et al., 2003;Radić and Hock, 2010;Church et al., 2011). To calculate current and future changes in glacier volumes, one of the most important data sets is a complete global glacier inventory with topographic attributes (e.g., minimum and maximum elevation, slope and aspect) for each glacier. However, recent calculations of sea-level rise due to MG&IC melting relied on estimates of the total volume of land ice scaled up from incomplete glacier inventories (e.g., Radić and Hock, 2010;Mernild et al., 2011). Hence, a globally complete and detailed glaciological data set is urgently required.
The World Glacier Inventory (WGI) compiled informa-tion on MG&IC from aerial photography, map and satellite imagery acquired during the 1960s and 1970s (WGMS, 1989). However, the WGI is not complete with respect to detailed glaciological information; only about 44% of MG&IC and 23% of their area are inventoried (Dyurgerov and Meier, 2005). Furthermore, the WGI makes change detection for individual glaciers nearly impossible because glaciers in the WGI are represented by point data rather than glacier outlines (WGMS, 1989). For this reason, the Global Land Ice Measurements from Space (GLIMS) project was designed to generate a global survey of digital glacier outlines from satellite data (Raup et al., 2007;Armstrong et al., 2012). To date the GLIMS Glacier Database, hosted at the National Snow and Ice Data Center (NSIDC) in Boulder, CO, USA, contains records for over 96792 glaciers with a total area of 325968 km 2 (Armstrong et al., 2012). Although glacier inventories for British Columbia, Caucasus, Heard Island, Iceland, Irian Jaya, and Switzerland are complete in the GLIMS Glacier Database, the inventories for other regions are not yet completed. Further, some inventoried glaciers in the GLIMS Glacier Database are lacking topographic parameters and some data are based on somewhat old maps (Armstrong et al., 2012). In addition, Cogley (2010) compiled an extended version of the WGI, called WGI-XF, through compilation of existing inventories including several older regional inventories that have been documented in WGMS (1989) but not in the WGI, and new inventories in Canada and the Sub-Antarctic. The WGI-XF contains records for over 131000 glaciers, covering approximately half of the MG&IC area (Cogley, 2010). Although glacier outlines can be derived from multispectral satellite data (e.g., Bayr et al., 1994;Sidjak and Wheate, 1999;Paul and Andreassen, 2009), these studies mainly focused on regional scales. Such studies allow estimation of the fluctuations in glacier volumes in these regions, but are insufficient for estimating global changes. A globally complete inventory of glacier outlines, the Randolph Glacier Inventory (RGI) has been published recently (Arendt et al., 2012). RGI primarily contains glacier outlines with little additional metadata for each record, and does not give the exact date of the determination of all glacier outlines. There is no question that global glacier outlines in different time series rather than snapshot and with date information are valuable.
Here we present a developed semi-automatic system to derive outlines of MG&IC (excluding glaciers in Greenland Correspondence to: Yong Zhang, Institute of Engineering Innovation, School of Engineering, the University of Tokyo, 2-11-16 Yayai, Bnkyo-ku, Tokyo 113-8656, Japan. E-mail: zhyong@sogo.t.u-tokyo.ac.jp ©2013, Japan Society fo Hydrology and Water Resources. and Antarctica) from global satellite data with 30 m horizontal resolution. The system consists of semi-automatic data archiving and an automated mapping system running in a Linux/UNIX system. With the automated mapping system, our data set can be updated as new satellite data become available and will form a time series of glacier outlines, rather than a snapshot of location and extent. Furthermore, apart from enhancing glacier inventories for regions where few direct measurements exist, the 30 m global database of MG&IC, could be integrated into global hydrological and land-surface models.

Satellite data
Landsat TM and Landsat ETM+ images were provided in an orthorectified version (UTM projection, WGS84 datum) by the Global Land Cover Facility (GLCF). The images were selected with as low snow and cloud cover as possible to be suitable for glacier mapping. Each image covers 185 km × 185 km and contains reflection intensity information for each band. Satellite images including the locations of MG&IG obtained from Hirabayashi et al. (2010) or the GLIMS Glacier Data were analyzed. Overall, we selected 334 Landsat TM images acquired during 1985-1997 and 415 ETM+ images acquired during 1999-2003.

ASTER GDEM
Specific topographic inventory parameters such as minimum, maximum, mean, and median elevation, mean slope, and mean aspect derived from a digital elevation model (DEM) were used for analysis of the obtained glacier outlines. For this study the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) with a horizontal resolution of 1 arc-second (~30 m) was used. It was released in July 2009 by the Ministry of Economy, Trade and Industry (METI), Japan, and NASA, covering the earth's surface between 83°N and 83°S.
As recent studies focusing on comparing a pre-release version of the ASTER GDEM and the DEM from the Shuttle Radar Topography Mission (SRTM) have shown (Hayakawa et al., 2008;Frey and Paul 2012), the ASTER GDEM has higher resolution, fewer missing data, and better topographic representation than the SRTM DEM, particularly in high mountainous areas with steep slopes. A report about typical artifacts and errors in the ASTER GDEM is given in METI/ NASA/USGS (2009), which indicates that the overall accuracy of the global ASTER GDEM can be taken to be approximately 20 m at 95% confidence interval, which is adequate for deriving glacier details for our global inventory.

Glacier inventory data
Glacier area on a continental scale from Radić and Hock (2010), the WGI-XF and the RGI were used to validate our estimate in each region (Table I). Radić and Hock (2010) calculated total glacierized area per region with complete glacier inventory from the WGI-XF data set and with incomplete glacier inventory from a global grid cell of glacierized area data set. The WGI-XF is an extended version of the WGI, in which a set of explicit specifications has been used to enhance the usefulness of the WGI data by eliminating low-level inconsistencies (Cogley, 2010). The RGI contains 170000 glacier outlines, which is a combination of both new and existing published glacier outlines (Arendt et al., 2012). In addition, glacier inventories for British Columbia, Caucasus, Switzerland, and the eastern and western sides of the Himalayas (excluding India and Pakistan) are complete in the GLIMS Glacier Database, and were used to further evaluate our estimates.

Mapping of MG&IC
Among methods for automated glacier mapping from multispectral satellite data, the band ratio method is considered the most efficient and convenient method (Paul, 2002;Raup et al., 2007;Paul and Andreassen, 2009). It involves creating a thresholded ratio image from the raw digital numbers of bands TM3 (0.63-0.69 μm) or TM4 (0.76-0.90 μm) and TM5 (1.55-1.75 μm), relying on the spectral contrast between ice and snow and other surfaces. As expected from the spectral properties, glacier ice has a moderate reflectivity in the TM4 band (fine-grained snow has high reflectance) and very low reflectance in the TM5 band. The division of TM4 by TM5 yields high and consistent values over ice and snow, and very low values for the surrounding terrain. The glaciers can be therefore distinguished from other surface types where the TM4/TM5 ratio is larger than 3.6 (Supplement Figure S1). Because of Table I. Total area of MG&IC in each region from the study of Radić and Hock (2010) (RH10), the WGI-XF (Cogley, 2010), the RGI (Arendt et al., 2012), and estimates from this study. spectral similarity with the surrounding terrain, some nonglacier areas at the edge of the image were identified as glaciers after applying the threshold TM4/TM5 > 3.6 (Supplement Figure S2). Also, some forested areas were identified as glaciers (Supplement Figure S2b). Hence, an additional threshold of TM3 (DN > 50) was applied to improve glacier mapping, exploiting the very high reflectivity of ice and snow in the TM3 band. The errors at the edge of the image were removed after applying the additional threshold in TM3 (Supplement Figure S2c). As a last step of the main processing, a 9 × 9 pixel average filter was applied to remove misclassification of pixels and noise, i.e. to eliminate isolated pixels and fill small gaps in glacier-covered regions (Supplement Figure S3). Finally, outlines of glaciers were obtained as an individual mass of contiguous glacier pixels. See Supplement Information S1 for the details.

Topographic information
Satellite-derived glacier outlines combined with the ASTER GDEM were used to derive topographic parameters for analyzing the glaciers, which include hypsometry, maximum, minimum and mean elevations. Note that the ASTER GDEM is the same spatial resolution as Landsat data (30 m), and also uses the same UTM coordinate system. However there are still some points that do not match between the two datasets. Although it is necessary to develop a mosaic from different scenes from different years, our study aimed to extract glacier outlines globally, and does not focus on the point or small catchment scale. This makes matching these points more difficult globally; hence we did not match them in the edges of imageries. As previous studies suggested (e.g., Frey and Paul, 2012;Hebeler and Purves, 2009), although some grid points that may have an influence on the hypsography of small glaciers are not equal to each other, local differences in these data have only a small influence on the hypsography of large glaciers.

RESULTS
In total, 414258 mapped glaciers with a total area of 449482 km 2 were identified in the 30 m global glacier database (Supplement Figure S4). Of all the regions defined (Table I), High Mountain Asia contains the largest ice mass outside of Antarctica and Greenland, representing ~42.6% of the total area of MG&IC, while the area in the Caucasus, Central Europe and New Zealand is relatively small, representing less than 1% of the total area.
The distribution of glacier area with elevation in Central Europe, the Caucasus, the Himalayas, North America and South America is depicted in Figure 1. Most of the glacier area occurs between 4500 m and 6000 m a.s.l. in the Himalayas, and glaciers in Central Europe and the Caucasus are distributed mainly in the altitudinal range of 2000 to 4000 m a.s.l. (Figure 1). In North America and South America, the elevation of most glaciers is relatively low, about 1500 m a.s.l., and some ice is also found between 3000 m and 5000 m a.s.l. (Figure 1).
The glaciers of each region are classified into size bins, and the total area per size bin is determined. Similar to previous studies (e.g., Raper and Braithwaite, 2005;Radić and Hock, 2010) we assign the upper boundaries for each area size bin to be 2 n km 2 with n ranging from −3 to 15. This means that the smallest size bin contains glaciers of less than 0.125 km 2 while the largest size bin contains glaciers between 16384 km 2 and 32768 km 2 . The cumulative number distribution of MG&IC for each size bin in Central Europe, the Caucasus, High Mountain Asia, North America and South America is shown in Figure 2. Our estimates were larger than those in the WGI-XF (Figure 2). Note that the WGI-XF does not cover all glaciers in these regions, and its average map year is 1964 with standard deviation of eleven years and a time range from 1901 to 1993 (Cogley, 2010;Radić and Hock, 2010), whereas our glacier database is derived from Landsat TM or ETM+ images acquired during recent decades. Furthermore, our results for small glaciers are more complete than other glacier inventories (e.g., WGMS, 1989;Cogley, 2010). In accordance with our estimation, small glaciers with area ≤ 1 km 2 in the Central Europe, Caucasus, High Mountain Asia, North America and South America account for 69% of the total number of MG&IC ( Figure 2). Especially in High Mountain Asia small glaciers account for ~44% of the total number (Figure 2e).

Error assessment
The accuracy of the mapped glacier outlines is an important but difficult topic in glacier mapping from satellite data (Raup et al., 2007;Racoviteanu et al., 2010), especially in regions with shadow, clouds, seasonal snow, proglacial lakes and debris cover. A comparison of glacier area from different databases indicated that the total area of MG&IC is fairly similar to previous estimates with the exception of Scandinavia and North and East Asia (Table I). Although images were selected with acquisition dates prior to the onset of heavy snowfall (approximately October 1 in the Northern Hemisphere and April 1 in the Southern Hemisphere) as late in the melt season as possible, only a few satellite images are acquired in May or June in Figure 1. Glacier area-elevation distribution in the Caucasus, Himalayas, Central Europe, North America and South America in 100 m elevation bins Scandinavia and North and East Asia. Consequently, seasonal snow may be identified as glacier pixels, resulting in overestimation in these regions. The relatively simple assumption of the onset of snowfall may lead to the overestimation of glacier area where snow season is out of the assumed periods, such as in the low latitude regions of the Southern Hemisphere where dry season is from May to August (Sicart et al., 2011). In addition, the area of debriscovered glaciers was underestimated in the western Himalaya and Caucasus where debris-covered areas are large, and are excluded from glacier area by this method. Although several methods for debris-cover mapping have been applied to different debris-covered glaciers (e.g., Bishop et al., 2000;Paul et al., 2004a;Racoviteanu et al., 2010), there is no single best method for debris-cover mapping that can be applied to large regions without some manual corrections of the resulting outlines. Combination of our method and state-of-the-art techniques to extract spatial distribution of debris using satellite information (e.g., Zhang et al., 2011) may improve our inventory over debriscovered glaciers.
The method used in this study is efficient in delineating clean ice in a timely manner, but might not work properly for ice in shadow. Ice in shadow needs to be delineated manually or using algorithms customized and tested for such cases (Racoviteanu et al., 2010). In addition, we applied a 9 × 9 pixel average filter to remove misclassification of pixels and noise, i.e. to eliminate isolated pixels and fill small gaps in glacier-covered regions, with the consequence that this process maybe lead to an underestimate of the area of small glaciers.

Area-altitude distribution for MG&IC
The altitude distribution as a function of area for glaciers Figure 2. Cumulative number of glaciers for each size bin in (a) North America, (b) South America, (c) Central Europe, (d) the Caucasus, and (e) High Mountain Asia. Blue and red colours denote our glacier database and the WGI-XF inventory has been modelled in previous studies to estimate glacier mass balance according to altitude (e.g., Raper et al., 2000;Paul et al., 2004b;Raper and Braithwaite, 2006;Hirabayashi et al., 2010). Because of limited hypsometric data for glaciers, these studies assumed relatively simple shapes to estimate the area-altitude distribution of glaciers. Raper et al. (2000) and Raper and Braithwaite (2006) assumed a linearly increasing function from the terminus to mean altitude and a linearly decreasing function above that altitude to approximate the area-altitude distribution (hereafter called the Triangle distribution). Hirabayashi et al. (2010) assumed that the altitude distribution of a glacier has a normal distribution curve that is shifted by one third of the range between the maximum and minimum altitudes. We evaluated the existing area-altitude models using our glacier database by calculating correlations for our estimate and the two areaaltitude assumptions (the Triangle distribution and 1/3 shifted normal distribution) for each 50 m altitude interval (Supplement Figure S5). The altitude distribution was obtained for the entirety of each glacier. The average correlations between the altitude distribution of glaciers and each assumed shape distribution in Central Europe, the Caucasus, the Himalayas, North America and South America shows that correlations are low in North America and South America (Table II). Overall, there is a better match with the Triangle distribution than with the 1/3 shifted normal distribution in these regions (Table II).

CONCLUSIONS
A 30 m global outline of MG&IC was derived from Landsat TM and ETM+ images, which are available for immediate use. Our database contains records for over 414258 glaciers with a total area of 449482 km 2 . Although there are errors due to limitations in satellite images and applied methods in some regions such as Scandinavia, North and East Asia, and the western Himalayas and Caucasus, our database matches well with previous satellite-derived estimates at regional scales, which is the appropriate scale needed for monitoring future fluctuations in MG&IC. Evaluation of published area-altitude models showed that the Triangle distribution resembles those obtained from our database in the Alps, the Caucasus, and the Himalayas. In particular, the developed 30 m glacier outlines can easily connect with other gridded data such as temperature or precipitation data to aid in the development of land surface modelling. We expect that our glacier outlines, which can be updated from recent satellite record in a semi-automatic way, can be used for a wide range of scientific purposes.