3.1. Criteria to identify surge of CYG

There are several distinctive surface features for identifying a surge-type glacier. These include repeated loops and folds in medial moraines, chaotically distributed and rapidly developing surface crevasses, potholes or surface lakes, and rapid advances of glacier snouts (Reference Meier and PostMeier and Post, 1969; Reference SturmSturm, 1987; Reference Barrand and MurrayBarrand and Murray, 2006; Reference CoplandCopland and others, 2011). These surface features make it easy to distinguish a surge-type glacier from a non-surge-type glacier by remote sensing. Many studies have used remote-sensing techniques to identify surge-type glaciers, such as with early use of aerial photographs (Reference PostPost, 1969) and recent studies based on satellite images (e.g. Reference Rolstad, Amlien, Hagen and LundenRolstad and others, 1997; Reference Murray, Strozzi, Luckman, Pritchard and JiskootMurray and others, 2002; Reference Jiskoot, Murray and LuckmanJiskoot and others, 2003; Reference Barrand and MurrayBarrand and Murray, 2006; Reference CoplandCopland and others, 2011). CYG was confirmed as a surge-type glacier based on some specialized surface features, i.e. contorted longitudinal foliations, rapid development of surface crevasses, and advances of glacier termini (see Section 4).

3.2. Improved cross-correlation feature-tracking (CCFT) method

The CCFT method, which is widely applied in remote-sensing-based glacier studies (e.g. Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992; Reference Rolstad, Amlien, Hagen and LundenRolstad and others, 1997; Reference Strozzi, Luckman, Murray, Wegmuller and WernerStrozzi and others, 2002; Reference Dowdeswell and BenhamDowdeswell and Benham, 2003; Reference Luckman, Murray, Jiskoot, Pritchard and StrozziLuckman and others, 2003; Reference Scherler, Leprince and StreckerScherler and others, 2008; Reference Mayer, Fowler, Lambrecht and ScharrerMayer and others, 2011; Reference Quincey, Braun, Glasser, Bishop, Hewitt and LuckmanQuincey and others, 2011), was used to collect image control points and extract glacier surface velocities. Most CCFT methods perform transformations on a satellite image by Fourier, wavelet or other methods prior to image matching, and obtain sub-pixel-scale accuracy by fitting the correlation coefficient surface around the matching pixel in the master image with two-dimensional (2-D) functions (e.g. Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992; Reference Strozzi, Luckman, Murray, Wegmuller and WernerStrozzi and others, 2002; Reference Leprince, Barbot, Ayoub and AvouacLeprince and others, 2007). Although this type of image-matching method can provide 1/10 to 1/20 pixel-matching accuracy, their algorithms are relatively complicated and have lesser efficiency.

We developed a more straightforward and simplified image-matching method in this study, using Interactive Data Language (IDL) programming tools, to extract matching points from a large number of Landsat image pairs more efficiently. It can also provide sub-pixel-scale image-matching accuracy. The maximum correlated matching window was first searched in one image band, without performing any transformation. We used a sub-pixel interpolation method rather than 2-D function fits to improve the accuracy of matching-point coordinates. The best matching-window result in the first step was divided into 1 0 x 1 0 sub-pixel positions, and the pixel value at each position was calculated by a bilinear interpolation method. The sub-pixel position having the maximum correlation coefficient within the matching window in the slave image was taken as the final matched position.

Although this method only finds the maximally matched position at 1/10 pixel scale, it is simpler and more efficient to attain a large number of image matches. The nominal 1/10 pixel-scale positioning accuracy is reasonable, so the method is capable of capturing small movements greater than 1/10 pixel size.

3.3. Landsat images and their orthorectification accuracy validation

Images acquired by TM and ETM+ on board Landsat-5 and Landsat-7 satellites, respectively, were used to study the CYG surge. All images were obtained from the US Geological Survey (USGS), and were orthorectified based on the Shuttle Radar Topography Mission digital elevation model (SRTM DEM) and ground control points from the Global Land Survey 2005 (GLS2005) dataset (USGS, 2011). In total, 19 Landsat images were used (Table 1), 12 of which are ETM+ images partially influenced by missing data strips due to failure of the satellite scan-line corrector.

Table 1. Landsat images used and their geometric accuracies. A l l accuracies were calculated by improved CCFT method, using Landsat TM image from 6 August 2010 as common reference image

Before retrieval of glacier surface speed, orthorectifica-tion accuracies of all images must be validated. The Landsat TM image from 6 August 2010, which has no bad strips and minimal snow and cloud cover among the images, was used as a common reference. All Landsat images around Yulinchuan glacier were made into subsets (spatial extent shown by the window in dark purple in Fig. 1), and at least 25 evenly distributed control points on them were collected using the CCFT method and minimum correlation coefficient 0.8. The first band of all Landsat TM and ETM+ images (0.45-0.52 |im; 30 m spatial resolution), which has greater contrast over normal land surfaces, was selected for the ground control point search. The search window in this process has a size of 51 x 51 pixels, whereas the matching window size varied from 17 to 35 pixels. Maximally correlated sub-pixel locations of all matching window sizes were preserved, and optimally matched locations were manually selected from the results. Results located on the glacier surface were also manually deleted. For ETM+ images, matched windows with bad strip area larger than 1/3 were skipped to minimize their influence. Residual error in the x and y directions and root-mean-square error (RMSE) of different images were calculated based on coordinates of the 25 matched control points.

Validation results (Table 1) illustrate that all Landsat images used have higher and more stable orthorectification accuracies, relative to the results of widely used manual orthorectification processes (half to one pixel size was generally achievable). Most RMSEs are less than half-pixel size, and some even reach 1/6 to 1/10 pixel. The larger RMSEs of two ETM+ images from 11 October 2008 and 4 March 2009 are due to thick snow cover, which may have caused failure in some of the matching points searched. Therefore, this study applies these USGS orthorectified Landsat images directly to the study of CYG’s surge, without additional corrections. However, we include the validated orthorectification errors in error assessment.

3.4. Detection of glacier area and length changes and error assessment

Changes in glacier area were based on glacier outlines. Outermost glacier outlines were automatically extracted by the widely used band ratio method (band 3/band 5; Reference Paul , Kaab, Maisch, Kellenberger and HaeberliPaul and others, 2002), adding some manual improvements for imagery with extensive snow and/or cloud cover. All these outlines were then divided by topographical ridgelines of Muztag mountain. The ridgelines were generated by procedures that automatically extract boundaries of small drainage basins using ArcInfo hydrologic analysis tools, and delete boundary sections in flat topography with reference to the aspect difference of pixels along both sides of the boundary, using a set of IDL procedures (Reference Guo, Liu, Yu and XuGuo and others, 2011). Glacier lengths at different dates were calculated by intersecting corresponding glacier outlines and two glacier flowlines, which were manually digitized by reference to the Landsat TM image of 6 August 2010.

Errors in glacier area and length and their changes were assessed by incorporating Landsat image resolution (30 m), relative orthorectification error (RMSE and x,y residuals; Table 1), and accuracy of the improved CCFT image-matching method (nominally ± 3 m). Errors in glacier length and its changes are calculated by

(1)

where EL is error in glacier length induced by image resolution (counted as half the pixel size, A) and orthorectification errors relative to the Landsat TM image of 8 August 2010 (er, corresponding to RMSE in Table 1), and ELC are errors in glacier length changes that propagate from errors in glacier length in the pre-change (subscript 1) and post-change (subscript 2) images. The largest errors in glacier length and its change via Eqn (1) are ±28.6 m and ±35.5 m, respectively, which are related to the Landsat image of 11 October 2008 with the maximum orthorectification error.

Errors in glacier area and its change are calculated following principles similar to the error estimation method of Reference Bolch, Menounos and WheateBolch and others (2010):

(2)

(3)

where E _A is error in glacier area from pixel resolution (E _a&ɛ , counted as half the pixel area) and the relative mis-orthorectification (E _a&ɛ , calculated summing x, y residual errors (&ɛ_rx , &ɛ_ry ) of each image in Table 1, and multiplying by pixel size), and includes area errors of all pixels along the glacier outline (with n the total, excluding pixels along ice divides). EAC is error in glacier area change, which is also calculated by area error in the pre-change (subscript 1) and post-change (subscript 2) images. The maximum area and area change errors from the above formulas are ±0.11 km2 and ±0.16 km2, respectively, which also correspond to the Landsat image of 11 October 2008.

3.5. Glacier surface speed derivation and error assessment

Glacier surface velocities were also derived by the improved CCFT method. Horizontal displacements of regular grid-points over the glacier surface with 5 x 5 pixel steps were extracted, based on Landsat images acquired on two successive dates. The near-infrared band (band 4; 0.76-0.90 &µ.m) of TM and ETM+, which has lower reflectance for glacier surfaces, was selected as the main data source in feature tracking. We calculated cross-correlation coefficients of a series of moving windows, with size increasing from 11 to 27 pixels in two-pixel increments. The maximum correlation coefficient was searched in a 45 x45 pixel window for every moving window size, and the best-matching window for each moving window size (with highest correlation coefficients) recorded. Selection of the searching and matching window sizes was dependent on manual measurement of maximum horizontal displacements of the glacier surface between 14 December 2008 and 16 February 2009 (with about 300 m or 10-pixel size). The larger matching window is more flexible for capturing surface features, but gives a smaller coefficient because of large glacier surface change during surge. In contrast, the smaller matching window can give a higher coefficient, but is less reliable. Thus, using multiple sizes of matching window is more appropriate for attaining the most accurate matching results. The sub-pixel match was then performed on the best-matched pixel-scale window of each size.

The maximally correlated sub-pixel position does not always correspond to the best point in representing the correct direction of glacier movements, because of the larger variance of the glacier surface during surge (sometimes it is even located on the upstream side of the glacier). Therefore, results of the CCFT method with different matching window sizes require further filtering. We used the glacier surface aspects calculated from SRTM to filter resulting points at each location, according to the fact that most movements of mountain glaciers are downslope. An index for each matching result was calculated to solve this problem:

(4)

where I is the index, R is the correlation coefficient of the matching window of different size, Rmax and Rmin are maximum and minimum coefficients among all matching window sizes, D is the absolute value of the difference between the direction of motion and aspect of the pre-move pixel, which was constrained within the range 0–&π, D _max and D _min are maximum and minimum D for all matching window sizes, and W_R and W_A are weights of the correlation coefficient and normalized angular difference, which were set to 0.6 and 0.4, respectively. The image-matching result with the largest I value was selected as the final matched location.

We also did sequential tracking of surface features, which are mainly composed of surface crevasses with certain spatial patterns, to obtain velocity and overall motion distance on different parts of the glacier surface. In total, 15 surface features were tracked, most of them formed by surface crevasses developed during the surge. Each was tracked by the same CCFT method and criteria as above, taking previously tracked positions on the post-move image as inputs to the next tracking process. The CCFT method is more sensitive to the structure of the matching window. Therefore, rather than absolute positions of surface features, inputs for the first track are roughly defined by visually determined centers of distinctive surface features formed by crevasses.

Accuracies of horizontal displacement and corresponding velocities were calculated as follows:

(5)

(6)

(7)

where E _d and E_D are displacement errors between two tracking dates and all tracked dates, respectively, E_v is accuracy of motion velocity, E _r1 and E _r2 are relative orthorectification accuracies of pre-move and post-move Landsat images (Table 1), u is uncertainty in the feature-tracking method (±3 m) and d is the number of days between two image acquisition dates.

4.1. Identification of CYG surge

There is one pre-surge surface feature that can be used to identify CYG as a surge-type glacier. Some contorted longitudinal foliations on the lower glacier snout (near the diverging point of CYG) can be distinguished from a fine-quality Landsat ETM+ image from 14 September 2004 (Fig. 2b), implying that at least one of its branches has surged in past decades. Furthermore, two types of evidence during the surge directly prove that CYG is a surge-type glacier. The first is that its termini advanced during 2008-10. As mentioned in Section 2, the united Yulinchuan glacier disintegrated between 1972 and 2004 (Fig. 2a and b). However, Landsat images acquired after 2007 (Fig. 2d-f) show that the termini of CYG have clearly advanced and reconverged into NYG and WYG, indicating that CYG has surged.

The second piece of evidence from the surge stage is from surface crevasse development. Similar to the CYG termini, the glacier surface has experienced dramatic changes, as clearly demonstrated by the Landsat image sequence (Fig. 2d-h). The surface was dominated by continuous and smooth ice before 2008, with some longitudinal features (see Fig. 2b and c). However, some transverse crevasses appeared just below the CYG equilibrium line on 11 October 2008, similar to those found on some Svalbard surging glaciers (Reference Sund, Eiken, Hagen and KaabSund and others, 2009). The surface crevasses then entered a stage of rapid development, and their area expanded quickly between 11 October and 14 December 2008. After that, crevasses developed more slowly, especially after 4 March 2009. The surface became stable after 7 September 2010 (Fig. 2f); the maximum range of surface crevasses was >6km, and the area reached 8 km² (Fig. 3a).

Fig. 3. Surface crevasse development sequence of CYG in different periods (a), and their distributions along (b) east and (c) west flowlines. Crevasse regions were manually digitized, referenced with corresponding Landsat images. Dates in legend are year/month/day.

From the planar, longitudinal and altitudinal distribution of crevasse developments (Fig. 3), it was discerned that the surface crevasses first developed on the middle glacier (∼5610 ma.s.l. and 5.9 km from the top), and then rapidly expanded to the upper and lower glaciers, especially along the west flowline development profile (Fig. 3c). The profile along the east flowline shows some difference (Fig. 3b). The regions around and above the ELA in this figure (just below the convergence point of the small northeast branch) were not seriously crevassed until 14 December 2008, although the first crevasse appeared near this region. Another difference is with the uppermost elevation of the crevasse distribution. The highest elevation of crevasses on the larger east branch of CYG was close to 5850 m a.s.l. (Fig. 3b), while the corresponding elevation of the west branch is ∼5750 ma.s.l. (Fig. 3c).

4.2. Changes in glacier area and length

As mentioned above, advance of a glacier terminus during surge is a primary indication of most surging glaciers. The length changes of CYG provide detailed information on its terminus advances (Fig. 4a). Glacier lengths at the two termini gradually increased between July 2007 and October 2008 (with mean rate of increase 0.07 ± 0 . 0 9 and 0.01 ± 0.09 m d–1 at termini A and B, respectively). Then there was a period of rapid increase until March 2009, in which termini A and B advanced about 472 ± 39 m (3.28 ± 0.3 md–1) and 1 0 9 ± 3 9m (0.76 ± 0.3 md–1) , respectively. Total changes of CYG length at the two termini reached 590 ± 26 m (5.1 ± 0 . 2% of the 2004 length) and 1 8 2 ± 2 6m (1.8 ±0.3%), respectively.

Fig. 4. Changes of (a) length and (b) area extracted by glacier outlines and two main flowlines of CYG during 2007–11. Dates are year/ month/day. Error bars (short black bars) were calculated by Eqns (1–3).

Changes in glacier area had temporal patterns similar to glacier lengths. In total, glacier area increased ∼1.41± 0.11 km2 (4.6 ± 0 . 4% of 2004 glacier area). The largest increase was between October 2008 and March 2009, ∼ 0 . 8 8 ± 0 . 1 8 km 2 (Fig. 4b). The area changes show some differences with lengths, mainly within stages before October 2008 and after March 2009. The glacier area showed a larger increasing trend, whereas glacier termini had relatively small advances. Comparing Landsat imagery on different dates, we found that these phenomena were induced by glacier expansions near the divergence point before October 2008, and widening of the glacier snout near terminus A after March 2009.

4.3. Changes in glacier surface velocities

Tracking of distinctive surface features over most of the surge phase shows that horizontal displacements in the lower and middle parts of the glacier surface had a much larger motion distance than the primary glacier terminus (terminus A in Fig. 1). The maximum tracked surface displacement (feature point B in Fig. 5; 13 July 2007 to 6 June 2011) is 1657±297 m. Identifiable motion distances of most other tracked feature points (A, B, C, E, F, G and I in Fig. 5) on the lower and middle parts are also close to or greater than 1000 m.

Fig. 5. Surface horizontal displacements of CYG, derived by sequential tracking of surface features during surge. Black characters represent codes of surface features, and red numbers are their tracked motion distances. Dates in legend are year/month/day.

Velocities of each surface feature in different stages were calculated by dividing the horizontal displacements by time intervals. Results clearly show that the glacier surface experienced dramatic changes in velocity during its surge. The surging phase for all velocity profiles was between July 2008 and July 2009, with the maximum surge from October 2008 to March 2009 (Fig. 6). The maximum surface velocity among tracked features was at point A during 12-28 November 2008, reaching 13 ± 1.5 md- 1 . The maximum speed of advance of terminus A was 12.1 ± 1 . 8 m d – 1 , but this was before the date of maximum speed of other feature points (27 October to 12 November 2008) on the eastern half of CYG. This implies that some undetected movements propagated to terminus A before July 2008. On the other hand, it may also suggest that the surge was initiated near the terminus.

Fig. 6. Velocity profiles of tracked surface features in Figure 5 on (a) eastern half and (b) western half of CYG, calculated from motion distances in corresponding period. Dates are year/month/day. Horizontal lines indicate mean velocities of corresponding feature points.

The velocity profile of these mobile surface features does not represent velocity changes at a fixed position on the glacier. To evaluate changes of motion velocities on different parts of the glacier, we derived velocities at regular gridpoints on its surface (excluding upper accumulation areas of CYG unaffected by its surge) with 5 x 5 pixel steps, retrieving velocity fields of the surface during various surge phases. We performed visual inspection and remnant fault point deletion in the results. Then velocities along two main glacier flowlines and three transverse profiles on different parts of the glacier (Fig. 7b; we manually excluded the converging point of upper transverse line T-A at 1 km from the western edge, and the divergence point of lower transverse line T-C at 1.2 km from south) were interpolated at 200 m intervals using the inverse distance weighted (IDW) interpolation method. Nearly all velocity profiles are discontinuous, especially the profile before 11 October 2008. This was due to the lack of valid velocity points, largely induced by lower contrast of the glacier surface generated through snowfalls before the surge, significant changes of the glacier surface during surge, or by bad strips in ETM+ images. However, the velocity profiles still provide some information about the surge process of CYG (Figs 8 and 9).

Fig. 7. (a) Example of velocity fields tracked at fixed gridpoints on glacier surface with 5 x 5 pixel steps. Dates are year/month/day. (b) Points for retrieving surface velocity profiles along two flowlines (FL-E and FL-W) in Figure 8, and three traverse lines (T-A, T-B and T-C) in Figure 9. Number labels are distances from heads of glacier (black) and from its west side (red).

Fig. 8. Surface velocity profiles along west (FL-W) and east (FL-E) flowlines in different periods: (a) 14 October 2006 to 27 October 2008; (b) 27 October 2008 to 16 February 2009; (c) 16 February 2009 to 7 September 2010.

Fig. 9. Surface velocity profiles along upper (T–A), middle (T–B) and lower (T–C) transverse lines in Figure 7b for different periods: (a) 14 October 2006 to 27 October 2008; (b) 27 October 2008 to 16 February 2009; (c) 16 February 2009 to 7 September 2010.

In contrast to surface feature tracking, tracks of fixed gridpoints on the glacier surface yielded few but valuable valid velocity points on lower parts of the CYG snout before July 2008, where the surface had many surface features visible on Landsat images (Fig. 2b and c). Therefore, the longitudinal velocity profiles provide more detail about the early surge stage, demonstrating that the surge occurred between May 2008 and July 2009. The primary surge phase was between October 2008 and February 2009, when maximum profile velocities generally exceeded 5 m d – 1 (Fig. 8b). Initiation and termination phases of the CYG surge were May–October 2008 and March–July 2009, respectively, during which profile velocities were generally close to or less than 3 m d–1.

Another noteworthy feature of all longitudinal velocity profiles is staged changes of velocities along two flowlines during the primary surge phase. The first stage was 11 October to 12 November 2008 (magenta triangle profile in Fig. 8a and dark yellow triangle profile in Fig. 8b). Surface velocities rapidly increased to the maximum of all profiles (11.6± 1.6 md–1, near 8.2 km of FL-E, dark yellow triangle profile). The second stage closely followed the first, through 4 March 2009 (theother profiles in Fig. 8b and gray profile in Fig. 8c). The surface velocity of the upper glacier (near 6.8 km of FL-E, dark blue hexagon profile in Fig. 8b) rapidly increased to the secondary maximum (10.6 ± 1.5 m d–1), and then the entire glacier entered a swiftly moving period. The average velocity of two profiles during 28 November to 14 December 2008 reached 4.8 ± 1.3 m d - 1 . After February 2009, surface velocities gradually decreased, returning to near stagnation in July 2009.

The third feature distinguishable from the longitudinal profiles is the surge wave propagation. The first wave appeared at the beginning of the initiating stage, when the maximum visible velocity occurred around 5.8 km on the east profile (dark blue triangle profile of FL-E in Fig. 8a). Then the wave propagated downward, with the maximum velocity arriving at 9.8 km on the east profile during 11-27 October 2008 (magenta triangle profile of FL-E in Fig. 8a). This wave continued to develop between 27 October and 12 November 2008, leading to the first massive surge of the lower glacier (below 6.8 km of the east flowline, dark yellow triangle profile in Fig. 8b). The second wave was a rapidly developing one, and originated from the upper glacier (above 6.8 km of the east flowline, dark blue hexagon profile of FL-E in Fig. 8b) during 12-28 November 2008. It quickly developed into a larger, second massive surge involving nearly the entire glacier, from 28 November to 14 December (purple star profiles in Fig. 8b). After that, glacier movements gradually slowed, entering a new quiescent phase through July 2009, with average surface velocity of all profiles <0.2 ± 0.3 md–1.

The transverse velocity profiles furnish more detail about development of the CYG surge. The lower traverse (T-C) shows that the lower glacier snout was essentially stable before 11 October 2008. In contrast, the upper and middle traverses (T-A and T-B) had higher velocities during July through October 2008 (Fig. 9a). This pattern changed during October to early November 2008. Velocities at the middle and lower traverses exceeded those at the upper traverse (Fig. 9b). This may be explained by the surge wave propagation mentioned above. The second surge wave initiated from the eastern part of the middle glacier after early November 2008, and rapidly propagated toward terminus A. This was followed by a larger, massive movement of most parts of the glacier.

The transverse profiles illustrate a very distinctive feature of the CYG surge: the east and west diverging branches experienced different motion processes during the surge. The east branch was more active during the first stage, and the maximum velocity of all profiles appeared on its center line between 27 October and 12 November 2008. During the second stage, velocities of the two branches were more evenly distributed than in the first stage, although at the start of this stage (12–28 November 2008) the east branch remained more active than the west.

4.4. On the surge history of CYG

An important characteristic of a surging glacier is the repeatability of surge activities. Most surging glaciers will resurge every 101–102 years (Reference Meier and PostMeier and Post, 1969). Some surging glaciers have recurrence intervals of 500 years (Reference SolheimSolheim, 1991). We studied the surge history of CYG and tried to retrieve the recurrence period. However, there is no sign of another surge between 1972 and 2007. The earliest data collected, the MSS image and topographical maps of the study region, were all from winter 1972. As mentioned in Section 2, the three branches of Yunlinchuan glacier converged in 1972 (Fig. 2a). Subsequent Landsat TM/ETM+ images from 1972–2004 (not included in Table 1) show that all termini of the glacier retreated continuously. Therefore, we can conclude that the recurrence period of CYG surging is not less than 40 years.