Mapping changes in the glaciers of the eastern Tienshan Mountains during 1977–2013 using multitemporal remote sensing

Abstract Because it has similar spectral characteristics as glaciers, snow around glaciers is a complicating factor in glacier inventory from satellite images. It is rare to acquire snowless and cloud-free satellite images in perennial snow mountain regions, and, moreover, glaciers are also usually shaded by mountain shadows. Therefore, a monotemporal satellite image can hardly map an intact glacier boundary. An object-oriented image segmentation method is proposed to map glaciers and their changes with multitemporal Landsat imagery. First, a “global-local” iteration segmentation method was used to delineate the snow and ice boundaries with each single phase image. Second, the mountain shadows of multitemporal Landsat images with different sun angles were mapped by methods of terrain analyses, and some parts of the glaciers inside these shadow regions can be identified. Finally, the complete glacier boundaries were restored with multitemporal images, and then these boundaries were intersected to get the minimum glacier extents. The method was tested with multitemporal Landsat images during 1977–2013 in the eastern Tienshan Mountains of Central Asia, and the spatial patterns and temporal process of these glacier changes in the last years were also analyzed. The results showed that the proposed method performs well in mapping glaciers in rugged alpine regions. Combining multiple images in different seasons provides a more effective way of removing disturbing factors during the process of glacier extraction. The total glacier area in the eastern Tienshan Mountains has decreased to 106.83     km 2 in 2013, with a loss of 21.5% since 1977, and the rate of glacier shrinkage accelerates in recent decades. The relationship between glacier recession and river runoffs was also analyzed, and the results showed that small glaciers whose areas are about 1.1     km 2 provide less runoff than those larger glaciers with about 2.5     km 2 .


Introduction
Continental glaciers in arid regions of Central Asia play an important role in water cycle and water resource management. 1 These glaciers provide important and valuable water resources for their downstream areas during drought years. Recent studies showed that glacier areas in the Tienshan Mountains have rapidly decreased in the last decade, leading to the decline of river discharge. 2 In addition, glaciers in the Tienshan Mountains have experienced dramatic changes and remain sensitive to temperature and precipitation variations. 3,4 Therefore, mapping glaciers and their changes will help us to understand the regional climate changes and water cycles in arid regions. Remote sensing provides us with more satellite images with various spatial and temporal resolutions, making glacier investigations more feasible and cost effective. 5 Many academic and research institutions have launched plans to map and monitor the glacier changes with remote sensing technologies. For example, a global glacier inventory was produced by global land ice measurements from space (GLIMS) in 2002. 6 Although numerous studies have focused on continental glacier changes in the western China and Central Asia, 7 few remote sensing images with good quality are available for glacier mapping due to cloud cover, snow cover, and mountain shadows. 8 Most of the continental glaciers in Central Asia are mainly distributed throughout the tops of rugged mountains, where clouds and snow persist throughout the year. Hence, the real boundaries of these glaciers are usually covered by snow or hidden by clouds, and they are hardly and precisely delineated with only one phase satellite image. Especially, for applications in glacier change investigations, snow and clouds will become the constraint factors to improve the mapping accuracies. In addition, continental glaciers are usually located at the top of mountains with rough terrains, and some parts of the glaciers might be inside the shadow of the mountains due to terrain shading, making them invisible on satellite images. At present, numerous remotely sensed glacier mapping methods are based on monotemporal satellite images, and classification or threshold segmentation-based methods 9 are the most common ways of glacier delineation. [10][11][12] Multitemporal and multispectral optical remote sensing images, such as Landsat TM, ASTER, provide an ideal data source for detecting glacier changes at appropriate spatial scales. [13][14][15] Normal difference snow index (NDSI) is a widely used indicator to enhance snow and ice while depressing other ground features, providing a simple but effective way of developing automated segmentation methods with multispectral images. 16 However, snow cannot differentiate from glaciers due to their similar spectral characteristics, and the thresholds for glacier segmentation can hardly get a stable value for different types of glaciers in various glacial environments. In order to map the glacier centerlines, slope information derived from digital elevation models (DEMs) is also used to track maximum local slope directions. 17 In the rugged mountain areas, glaciers are usually located in the shaded regions, and some parts of the glaciers are hidden by mountain shadows. The extent and spatial distribution of shadows vary with seasons due to the changing angles of the Sun. Therefore, the effects of mountain shadows to glacier mapping should be evaluated by terrain analysis with DEM data.
The objective of this paper is to develop a glacier mapping method based on multitemporal Landsat images. It combines Landsat images with different seasons to remove the effects of snow cover and cloud cover, to get the minimum glacier extents. Meanwhile, mountain shadows derived from multitemporal Landsat images with different sun angles are intersected to identify some parts of the glaciers inside the shadow regions. The method was also tested with multitemporal Landsat images during 1977-2013 to map the glacier changes in the eastern Tienshan Mountains, and the spatial patterns and temporal process of these glacier changes in the last decade were also analyzed.

Study Area
The study area of this paper is located at the Mt. Karlik, Hami District, Tienshan Mountains ( Fig. 1), extending from 42.95°N to 43.20°N and 94.15°E to 94.65°E, covering most of the glaciers in this area. According to the glacier inventory by GLIMS, there are a total of 156 glaciers covering 136.8 km 2 . The mountain runs from northwest to southeast, and Tomort Peak is its highest mountain peak with an altitude of 4886 m. Most of these glaciers are valley glaciers or hanging glaciers whose surrounding terrains are rugged and steep. 18 Some parts of the glaciers are easily shaded by great shadows of high mountains.

Data
Multiple Landsat remote sensing images are used to map the glacier changes, including Landsat 2 multispectral scanner (MSS), Landsat 5 thematic mapper (TM), Landsat 7 enhanced thematic mapper plus (ETM+), and Landsat 8 operational land imager (OLI) ( Table 1). Most of these images are snow-free or cloud-free images. These Landsat products are released by the U.S. Geological Survey (USGS) and can be downloaded from the USGS Global Visualization Viewer in Ref. 19. They are orthorectified Landsat products, and the quantity and quality of these valuable resources are sufficient for monitoring environmental change.
There are three to four scene images of not more than 2 years for each glacier mapping. Using the Landsat 5 TM image of August 24, 2001, as a base, all other images are co-registered and projected to the WGS-84 ellipsoid and the Albers Equal Area Conic projection system. The mountain shadows are obtained from a DEM product with a resolution of 90 m of the shuttle radar topography mission (SRTM) that is most likely sufficient for glacier mapping elsewhere in the world.

Methodology
The glaciers in the study area are classified as continental glaciers. These glaciers are developed in a continental climate with little precipitation and weak glacial movement. The fastest retreat rate in the eastern Tienshan Mountains occurred in the NO. 4 glacier of the Sigong River on the Bogda Mountain, and the retreat rate of the glacier length was 13.3 m∕annual during 2006-2009. 20 According to this shrinking rate, the glacier can be considered to be invariant over 2.26 years when it mapped with a 30 m resolution Landsat image. As a result, if a glacier is mapped by the Landsat images within 2 years, its boundary can be taken as stable without any changes.
The schematic flow of the multiangle method is shown in Fig. 2. In this paper, the glacier extraction method is based on the following two assumptions: 1. The part of the glacier boundaries hindered by mountain shadows can be identified and mapped by multitemporal Landsat images with multiple sun angles (Fig. 2, Step 2). 2. The glacier boundary during a short period (2 years) is invariant, mapping glaciers by combining the multiple temporal Landsat images get the minimum extents of snow and ice, hence the glacier boundary covered by snow can be identified (Fig. 2, Step 3).

Mountain shadows
The rugged terrain of Mt. Karlik has a strong shadow effect on the glaciers on the Landsat images, and some parts of the glaciers are invisible due to the shields against the sunlight [ Fig. 3(a)]. Different sun angles have different shadow distributions on the image. The regions of mountain shadows can be calculated by topographic models with SRTM DEM; therefore, the influence of mountain shadows can be eliminated by viewing the snow/ice in images with different sun angles. The glaciers on the Tienshan Mountains are primarily located in the highest relief area at the altitude of more than 4000 m, and most of the glaciers are in the northern parts of the mountain, as less sun or heat irradiates directly to avoid rapid glacier melting. The distribution of mountain shadows is determined by the incident angle of the Lambert cosine law, which is determined by calculating the normal lines of the digitized slope surface elements. In the terrain modeling tools of ENVI (The Environment for Visualizing Images), the SRTM DEM data and the sun angle are imported to generate the corresponding shaded relief image with a range of R ∈ ½0;1. The mountain shadows [Figs. 3(b) and 3(c)] are extracted from the shaded relief image while R ¼ 0. The shaded relief image uses a shadow principle to describe the terrain. When R ¼ 0, it means the area is completely sheltered by the mountain. 21 The Landsat image is overlaid with red vectors of mountain shadows [ Fig. 3(d)] to show that mountain shadows extracted from DEM are accurate.

Snow/ice
Snow/ice has a low reflection in the mid-infrared band (wavelength of approximately 1.5 μm) and high reflection in the green band; therefore, the normalized difference snow index (NDSI) uses these two bands for ratio calculation. These NDSIs are known for their ability to enhance the snow/ice features and simultaneously depress other ground features. 16 The top-ofatmosphere (TOA) reflectance of bands 2 and 5 of the Landsat TM/ETM+/OLI images is used to calculate the NDSI as following equation: ρ B2 and ρ B5 are the TOA reflectance of bands 2 and 5, respectively, of the Landsat TM/ETM+/ OLI images. The NDSI algorithm and the thresholds have been verified to distinguish snow/ice  from clouds, vegetation, or water bodies, but a global threshold is required to produce different accuracies for diverse local areas. Moraine often covers the terminus of glacier, making the segmentation threshold unstable.
A "global-local" algorithm is applied to extract the snow/ice information from Landsat imagery. 22 To ensure that all the snow/ice pixels are delineated, an initial threshold (here the threshold T ¼ 0.4 is small enough to extract all glaciers) is applied to segment the NDSI image (Fig. 4), dividing the image into snow/ice and background pixels. The normalized difference water index (NDWI) threshold of −0.2 is set to remove the water misclassified into snow/ice. 16 ρ B2 and ρ B4 are the TOA reflectance of bands 2 and 4 of the Landsat TM/ETM+/OLI images, respectively.
In the previous section, all snow/ice boundaries are globally segmented with the same threshold, which leads to an inaccurate determination of the snow/ice boundaries due to the value variations of the different NDSI images. Figure 4 showed a flow chart of the "global-local" glacier segmentation method. Figure 4 is a glacier at the terminus of the ablation season in the image acquired on August 24, 2010 (the orbit number is 147/35, the center position is 43.4°N, 93.3°E); Fig. 4(c) is a snow-affected glacier in the image acquired on September 16, 2007 (the orbit number is 148/31, the center position is 41.7°N, 78.8°E). The red boundary is the snow/ice boundary obtained through global segmentation: the green boundary is the buffer boundary, the yellow boundary is the snow/ice boundary obtained through iterative local segmentation, and the pink boundary in Fig. 4(c) is the water boundary removed by NDWI.
A global threshold value (Global T ¼ 0.4) is applied to segment the glacier from other ground features, while some parts of the glacier pixels are missing and the glacier boundaries are not intact. When a buffer zone is built based on the initial segmentation results, each glacier is further segmented in its own buffer region. A histogram segmentation is applied, and each glacier has different local segmentation values. However, according to segmentation results, there are still some disturbing factors; for example, the snow/ice boundaries cannot be delineated when there are mountain shadows around its regions [ Fig. 4(c)].

Mountain Shadow Removal
Mountain shadows have similar spectral characteristics as snow/ice; therefore, it is difficult to differentiate them from the snow/ice in the spectral domain. 23 The glaciers in Mt. Karlik are mainly distributed in the rugged mountain regions where mountain shadows are densely distributed. Mountain shadows are intermixed with glaciers over the mountain, which leads to a large amount of manual editing work. 24,25 Therefore, the multiangle method to eliminate the interference of mountain shadows is required for large areas.  A GIS spatial overlay analysis is used to obtain the glacier in the total shadows area of the multiangle images: (1) The sum shadows of multitemporal images are calculated out to determine the total shadows area S max of remote sensing images with different sun angles SðiÞ; when there are N temporal images for excluding mountain shadows, SðiÞ represents the shadows of temporal i images, i ∈ ½1; N. (2) In order to extract the glaciers from the snow/ice in the total shadows area, the snow/ice G S ðiÞ of temporal i image in the total shadows area is calculated as follows: G t ðiÞ is the snow/ice extracted by the method of Part 4.1 of corresponding temporal i image. Then, the glaciers in the total shadows area G S can be obtained by (3) Finally, the mountain shadows of each image are eliminated. The snow/ice GðiÞ excluding the shadows of temporal i image is as follows: SðiÞ ¼ ϕ, it means the multiangle images can completely eliminate the mountain shadows. If not, mark the remaining shadows not eliminated by the multiangle images.

Multitemporal Glacier Intersection
The distribution of snow varies with elevation, aspect, and other geography features; therefore, snow around the glaciers makes glacier mapping results uncertain. In the absence of meteorological data, remote sensing images cannot determine which locations are not disturbed by snow. 26 The retreat of the Tienshan glacier is an overall trend. The most dramatic state is deemed to be an ideal glacier boundary. 27 In this paper, the minimum boundaries G Multi of the multiangle snow/ice boundaries Gð1Þ; Gð2Þ; Gð3Þ; · · · ; GðNÞ are calculated as the ideal glacier boundary (Fig. 6). G Multi ¼ Gð1Þ ∩ Gð2Þ ∩ Gð3Þ ∩ · · · ∩ GðNÞ:

Result Assessments
Glaciers extracted from four scenes of Landsat images during 2006-2007 listed in Table 1 are used to verify the validity of this multiangle method. The SRTM DEM and the sun angles of the four images are used to extract the mountain shadows, and the shadow intersection of the four images is equal to a null set. The snow/ice boundaries of each scene are extracted by the "globallocal" threshold method of NDSI. After overlaying the multiangle snow/ice boundaries and the mountain shadows information, the snow/ice boundaries influenced by the mountain shadows are determined. Using the scene of September 7, 2006, as an example, the glaciers covered by the mountain shadow are detected by the information of the snow/ice and mountain shadows of three sun angles, and the four scenes of snow/ice boundaries are not affected by the mountain shadows. The snow/ice boundaries of the four scenes are intersected to obtain the minimum boundary as the glacier boundary G Multi . Glaciers disturbed by shadow and snow (Fig. 7) are used to assess the effect of the multiangle method on each temporal image (

Glacier Change Overall
The contiguous glaciers are divided into their drainage basins to obtain a glacier inventory. An automated approach is presented and derives basins from the SRTM DEM based on hydrological analysis around each glacier. 29 According to the results, the overall glacier area has shrunk from 136.84 km 2 in 1977 to 106.83 km 2 in 2013 (Fig. 8). Many of the observed changes (growing rock outcrops, tongue separations, formation of pro-glacial lakes, and albedo lowering) during 2007-2013 are related to positive feedbacks that accelerate further glacier disintegration. 30 It is found that glaciers of Mt. Karlik are retreating with varying rates (Table 3) in the different decades, with a faster glacier retreat occurring during 2007-2013. A total of 30.01 km 2 of glacier areas were lost during the past 36 years, with a 21.93% loss from 1977 to 2013. The overall retreat rate is estimated as 0.6% by year. Glacier samples are selected for the change analysis. The split glacier polygons in 1991 provide the basic input for the analysis with 80 glaciers as the glacier samples. Figure 9 shows the areas of glacier samples according to size classes. The areas, numbers and loss areas during 1977 and 2013 in Fig. 9 are classified in sizes, but loss areas are calculated based on individual glacier areas in 1977. The strong retreating is moving toward the class of 2 to 3 km 2 with an area loss that is 43.1% of total area loss from 1977 to 2013, and a class larger than 5 km 2 has a second area loss of 35.4%. This behavior implies that larger glaciers are still major contributors to the river runoff despite their number being limited (Fig. 9). The number of the class less than 0.5 km 2 Fig. 8 Glaciers of 1977 and 2013. The background is the DEM of this area. increases by being synchronously converted from the larger class, ranging from a minor advance to complete disappearance; furthermore, this class has a rapid response time to climate perturbations, and these systems exhibit instability with highly individual area changes compared to larger glacier systems, as their dynamical response time can be many decades. 31

Glacier Changes in the Sub-Basins
There are 12 glaciers with a mean area of 1.1 km 2 distributed in the 5Y822A Basin, and 9 glaciers with a mean area of 2.5 km 2 in the 5Y822B Basin in 1972. In this study, the sub-basins 5Y822A and 5Y822B reduced by 24.6% and 17.4%, respectively, from 1977 to 2013 (Fig. 10).
Although investigations of this region were conducted in Refs. 32 and 33 (Table 4), the trends of glacier shrinkage are obvious during the past three decades. It is unlikely that the recent trend of glacier wastage will stop in the near future.

Discussion
Climate-related factors on the regional scale mainly influence the glacier retreat, such as the temperature rise and the change of precipitation allocation between summer and winter.
Since the 1980s, the climate in northwest China has been changed to warmth and wetness. 34 For glacier-free basins, such as at the Toudaogou Hydrological Station, the annual variability of runoff is high. However, for glacier-covered basins, such as Baiji and Yushugou, although the retreat of the glaciers has caused variability of the runoff, the effect is relatively small (Fig. 11).
For the lower glacier-covered catchment, such as Baiji, the discharge is decreased because of the strong consumption of small glaciers during the past decades. For the catchment of Yushugou, the runoff tends to increase due to the larger glaciers.
The enhanced glacier melting produces increased river runoff; then the runoff starts to decrease after glaciers reach a critical size. Glacier cover is likely to control the interannual variability of runoff, and the variability is the lowest at a moderate percentage of glacier cover; the variability increases as the glacier cover both decreases and increases. 35 The main difference between glaciered and glacier-free catchments is that the runoff from a glacierfree basin is dominated by precipitation, whereas glaciered basins are temperature dominated. This behavior highlights the importance of glaciers because they produce the most water during hot and dry periods when the precipitation is lacking.
Generally, changes in summer temperatures appear to drive the increased rates of glacier loss, while the increase of the winter temperature and its contribution to glacier melting should be considered. The increase of the winter temperature prolongs the time that the glacial surface is close to 0°C. Because the snow layer temperatures are high, the same amount of melting glacier consumed a lower amount of energy, thereby increasing the sensitivity of the glacier to air temperature rise. 36 An increase in precipitation creates a favorable condition for the accumulation of glaciers. Because the main glacier ablation period is in summer, the increase of the summer temperature causes an increase in the liquid precipitation but not in the solid precipitation, thus causing the glacier mass income to decrease and the glacier retreat to excelerate. 37 Mountain glaciers can modulate the river runoff and maintain a balance of the hydrological cycle. Glaciers are directly able to delay the runoff by preventing precipitation from running off. Such storage of  precipitation occurs in summer with snow accumulating and melting on the glacier, and it moderates when glaciers encounter heavy rains. Thus, the water is stored as snow/ice in glaciers and is melted depending on climate. 38 Figure 12 shows an illustration that the temperature is rising at the studied region [ Fig. 12(a)] and the precipitations in Yushugou and Baiji nearly maintain their borders [Figs. 12(b) and 12(c)], while Toudaogou is slightly increased [ Fig. 12(d)]. This behavior may indicate that the runoff of Toudaogou expresses the increase with high variability.
Increasing temperature appears to account for recent glacier shrinkage. Thus, the strong ice mass loss leads to an increase of the interannual variability of river runoff. In the coming decades, a continuous increase in the glacial runoff can be reasonably expected in response to warming, especially in spring and early summer.

Conclusions
For the snow and mountain shadow covering, a multiangle method is applied to map the glaciers during the invariant period of glaciers in 1977, 1991-1992, 2001, 2006-2007, and 2013. Some conclusions of this study are as follows: 1. A "global-local" iteration extraction method of NDSI can increase the glacier boundary accuracy of the discrepant local area. The glacier and mountain shadow information obtained from multiple sun azimuths and sun elevations can determine the glacier boundaries affected by the mountain shadows. In addition, multitemporal glacier information can remove the snow-disturbed areas. In Table 2, during 2006 and 2007, the date disturbed by the minimum mountain shadow and the date of the minimum snow cover are different. The multiangle method can combine the two best observation conditions of remote sensing, thereby to extract the glacier boundaries that are closest to the real glacier boundaries. 2. The method is used to report the changes of the glacier boundaries on Mt. Karlik in eastern Tianshan from 1977 to 2013. The investigated glaciers witnessed climate warming during the past four decades and reduced their total area by 30.01 km 2 , or a relative change of −21.9%. A large quantity of the glaciers of 2 to 3 km 2 relatively lost more area. The larger glaciers mainly contributed to river runoff. 3. Glaciers in the upper stream of a river play an important role in collecting solid water in winter and releasing it in the form of water. Even a low fraction of glacier covered within a basin has a tremendous impact on the hydrology. The runoff is decreased when the glacier tends to become smaller or vanish. In regions where glaciers are primarily larger than 2 km 2 , the runoff may increase in the future with the present climate warming situation. 4. With different water supplies of upstream glaciers, a great disparity in the trends of river runoff appears. For glacier-free basins, such as the Toudaogou hydrological station, the annual variability of runoff is high. While for the glacier-covered basin, such as Baiji and Yushugou, the retreats of glaciers cause the variability of the runoff. For the lower glacier-covered catchment, such as Baiji, the discharge is decreased because of the strong consumption of small glaciers during the past decades. For the catchment of Yushugou, the runoff tends to increase due to the existence of larger glaciers.