1.

Introduction

Continental glaciers in arid regions of Central Asia play an important role in water cycle and water resource management.¹ 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.² 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.⁵ 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.⁶

Although numerous studies have focused on continental glacier changes in the western China and Central Asia,⁷ few remote sensing images with good quality are available for glacier mapping due to cloud cover, snow cover, and mountain shadows.⁸ 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⁹ are the most common ways of glacier delineation.^10–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–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.¹⁶ 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.¹⁷ 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.

2.

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{\mathrm{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.¹⁸ Some parts of the glaciers are easily shaded by great shadows of high mountains.

Fig. 1

The location of the investigated area.

3.

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.

Table 1

Landsat data used in the study.

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.

4.

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\mathrm{m}/\text{annual}$ during 2006–2009.²⁰ 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:

Fig. 2

Schematic flow illustrating the multiangle glacier extraction method.

4.1.

Single-Scene Information Extraction

4.1.1.

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\in [\mathrm{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.²¹ The Landsat image is overlaid with red vectors of mountain shadows [Fig. 3(d)] to show that mountain shadows extracted from DEM are accurate.

Fig. 3

Automatic extraction of the mountain shadows. (a) DEM data. (b) Shaded relief image overlaid with red vectors of mountain shadows. (c) Landsat image combined with bands 7, 4 and 2 as RGB. (d) Landsat image overlaid with red vectors of mountain shadows.

4.1.2.

Snow/ice

Snow/ice has a low reflection in the mid-infrared band (wavelength of approximately $1.5\mu \mathrm{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.¹⁶ The top-of-atmosphere (TOA) reflectance of bands 2 and 5 of the Landsat TM/ETM+/OLI images is used to calculate the NDSI as following equation:

Eq. (1)

$$\mathrm{NDSI}=\frac{{\rho }_{\mathrm{B}2}-{\rho }_{\mathrm{B}5}}{{\rho }_{\mathrm{B}2}+{\rho }_{\mathrm{B}5}}.$$

${\rho }_{\mathrm{B}2}$ and ${\rho }_{\mathrm{B}5}$ 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.²² 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.¹⁶

Eq. (2)

$$\mathrm{NDWI}=\frac{{\rho }_{\mathrm{B}2}-{\rho }_{\mathrm{B}4}}{{\rho }_{\mathrm{B}2}+{\rho }_{\mathrm{B}4}}.$$

Fig. 4

Snow/ice extraction based on the “global-local” segmentation and their NDSI histograms. (a), (b) and (c) in the fist row are the Landsat images combined with band 7, 4 and 2 as RGB. The second row is NDSI images overlaid with green buffer vectors and red glacier vectors segmented by global thresholds. The third row is NDSI histograms and red vertical lines of local thresholds. The fourth row is Landsat images overlaid with yellow glacier vectors segmented by local thresholds.

${\rho }_{\mathrm{B}2}$ and ${\rho }_{\mathrm{B}4}$ 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(a) is a larger glacier in the image acquired on August 3, 2011 (the orbit number is 147/35, the center position is 35.7°N, 78.5° E); Fig. 4(b) 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 ($\text{Global}\_\mathrm{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)].

4.2.

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.²³ 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. Figures 5(a)–5(c) show the snow/ice and shadows of temporal 1, temporal 2, and temporal 3 images, respectively, with different sun elevation angles and sun azimuth angles: 28.3047, 160.1330; 45.0772, 149.8810; 51.3661, 138.7687 deg. Figure 5(d) shows red glacier vectors on September 25, 2001, removed shadow impact using images acquired on September 17, 2001 and August 24, 2001.

Fig. 5

Removal of mountain shadows in Landsat image acquired on (a) September 25, 2001, using images acquired on (b) September 17, 2001, and (c) August 24, 2001. (d) shows red glacier vectors on September 25, 2001, removed shadow impact overlaid on corresponding image.

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

  Eq. (3)

  $$S_{\max }=\bigcup _{i=1}^{N}S(i),$$
  when there are $N$ temporal images for excluding mountain shadows, $S(i)$ represents the shadows of temporal $i$ images, $i\in [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:

  Eq. (4)

  $$G_S(i)=G_t(i)\cap S_{\max }.$$

  $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

  Eq. (5)

  $$G_S=\bigcap _{i=1}^N[G_S(i)\cup S(i)].$$

• (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:

  Eq. (6)

  $$G(i)=G_t(i)\cup G_S.$$

If $\bigcap _{i=1}^{N}S(i)=\varphi$, it means the multiangle images can completely eliminate the mountain shadows. If not, mark the remaining shadows not eliminated by the multiangle images.

4.3.

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.²⁶ The retreat of the Tienshan glacier is an overall trend. The most dramatic state is deemed to be an ideal glacier boundary.²⁷ In this paper, the minimum boundaries $G_{\mathrm{Multi}}$ of the multiangle snow/ice boundaries $G(1),G(2),G(3),\cdots ,G(N)$ are calculated as the ideal glacier boundary (Fig. 6). Figure 6 shows glacier boundaries [Fig. 6(e)] calculated from the minimum snow/ice boundaries of four temporal images acquired on August 22, 2006 [Fig. 6(a)], September 7, 2006 [Fig. 6(b)], July 24, 2007 [Fig. 6(c)], and August 25, 2007 [Fig. 6(d)].

Eq. (7)

$$G_{\mathrm{Multi}}=G(1)\cap G(2)\cap G(3)\cap \cdots \cap G(N).$$

Fig. 6

Glaciers (e) calculated from snow/ice boundaries of four temporal images acquired on (a) August 22, 2006, (b) September 7, 2006, (c) July 24, 2007, and (d) August 25, 2007.

4.4.

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 “global-local” 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_{\mathrm{Multi}}$.

Glaciers disturbed by shadow and snow (Fig. 7) are used to assess the effect of the multiangle method on each temporal image (Table 2).²⁸

Fig. 7

The glacier of 2006–2007 extracted by the multiangle method.

Table 2

The multiangle glacier mapping effect.

Date	Glacier disturbed by shadow		Snow disturbed area
	Area (${\mathrm{km}}^2$)	Proportion (%)	Area (${\mathrm{km}}^2$)	Proportion (%)
August 22, 2006	13.63	13.1	45.48	12.9
September 7, 2006	31.83	30.6	41.89	30.3
July 24, 2007	0.104	0.1	39.81	28.8
August 25, 2007	17.16	16.5	63.86	46.2

5.

Glacier Changes

5.1.

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.²⁹ According to the results, the overall glacier area has shrunk from $136.84{\mathrm{km}}^2$ in 1977 to $106.83{\mathrm{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.³⁰ 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{\mathrm{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.

Fig. 8

Glaciers of 1977 and 2013. The background is the DEM of this area.

Table 3

Changes in glacier area of Mt. Karlik.

Year	Area (${\mathrm{km}}^2$)	Loss area (${\mathrm{km}}^2$)	Retreat rate ($\%/\text{annual}$)
1977	136.84	10.15	0.49
1992	126.70	5.85	0.51
2001	120.84	2.29	0.32
2007	118.55	11.72	1.67
2013	106.83

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{\mathrm{km}}^2$ with an area loss that is 43.1% of total area loss from 1977 to 2013, and a class larger than $5{\mathrm{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{\mathrm{km}}^2$ 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.³¹

Fig. 9

Area and count of glaciers classified by size (in ${\mathrm{km}}^2$: $<0.5$, 0.5–1, 1–2, 2–5, and $>5$).

5.2.

Glacier Changes in the Sub-Basins

There are 12 glaciers with a mean area of $1.1{\mathrm{km}}^2$ distributed in the 5Y822A Basin, and 9 glaciers with a mean area of $2.5{\mathrm{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).

Fig. 10

Glacier boundaries from 1977 to 2013 in (a) Basin 5Y822A, (b) Basin 5Y822B.

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.

Table 4

Glacier changes in different basins.

Sub-basins	Area (${\mathrm{km}}^2$)		Area change (%) (1977–2013)	Area change (%) (1972–2005)a	Area change (%) (1972–2005)b
	1977	2013
Basin 5Y822A (Baiji)	13.43	10.13	$-24.6$	$-14.1$	$-7.5$
Basin 5Y822B (Yushugou)	22.38	18.49	$-17.4$	$-11.7$	$-6.1$

^adata of the column in Ref. 32.

^bdata of the column in Ref. 33.

6.

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.³⁴

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

Fig. 11

Time series of the runoff deviation from the average during 1977 to 2010 at (a) Toudaogou, (b) Baiji hydrological stations, and during 1979 to 2010 at (c) Yushugou hydrological station.

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.³⁵ The main difference between glaciered and glacier-free catchments is that the runoff from a glacier-free 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.³⁶

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.³⁷ 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.³⁸

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.

Fig. 12

Time series of the temperature of Hami during the period of 1977 to 2013 and the precipitation of the Yushugou, Baiji, and Toudaogou hydrological stations. (a) Temperature of Hami, (b) precipitation of Yushugou, (c) precipitation of Baiji, (d) precipitation of Toudaogou.

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.

7.

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: