Four decades of glacier variations at Muztagh Ata (eastern Pamir): a multi-sensor study including Hexagon KH-9 and Pléiades data

Abstract. Previous in situ measurements have indicated a slight mass gain at Muztagh Ata in the eastern Pamir, contrary to the global trend. We extend these measurements both in space and time by using remote sensing data and present four decades of glacier variations in the entire mountain massif. Geodetic mass balances and area changes were determined at glacier scale from stereo satellite imagery and derived digital elevation models (DEMs). This includes Hexagon KH-9 (year 1973), ALOS-PRISM~(2009), Pleiades (2013) and Landsat 7 ETM+ data in conjunction with the SRTM-3 DEM (2000). In addition, surface velocities of Kekesayi Glacier, the largest glacier at Muztagh Ata, were derived from amplitude tracking of TerraSAR-X images (2011). Locally, we observed strong spatial and temporal glacier variations during the last four decades, which were, however, on average not significant for the entire massif. Some south-west-exposed glaciers fluctuated or advanced, while glaciers with other aspects rather experienced continuous shrinkage. Several glaciers such as Kekesayi indicate no measurable change at their frontal position, but clear down-wasting despite mostly thick debris coverage at low altitudes. The surface velocity of this debris-covered glacier reach up to 20 cm per day, but the lowest part of the tongue appears to be stagnant. The low velocity or even stagnancy at the tongue is likely one reason for the down-wasting. On average, the glaciers showed a small, insignificant shrinkage from 274.3 ± 10.6 km2 in 1973 to 272.7 ± 1.0 km2 in 2013 (−0.02 ± 0.1 % a−1). Average mass changes in the range of −0.03 ± 0.33 m w.e. a−1 (1973–2009) to −0.01 ± 0.30 m w.e. a−1 (1973–2013) reveal nearly balanced budgets for the last 40 years. Indications of slightly positive rates after 1999 (+0.04 ± 0.27 m w.e. a−1) are not significant, but confirmed by measurements in the field.


Introduction
The glaciers of high-mountain Asia revealed an average mass loss of −26 ± 12 Gt a −1 from 2003 to 2009, which, however, is affected by strong regional variations (Gardner et al., 2013).Heterogeneous glacier mass balances in Pamir and Karakoram are confirmed at least for the last decade (Bolch et al., 2012;Kääb et al., 2012;Yao et al., 2012;Gardelle et al., 2013).Glaciers in the Pamir continued to retreat and shrink on average, while several glacier surges were observed at the same time (Kotlyakov et al., 2008;Unger-Shayesteh et al., 2013).Since 1999, slight mass gain of +0.10 ± 0.16 m w.e. a −1 was measured in central Karakoram, and of +0.14 ± 0.13 m w.e. a −1 in west Pamir using SPOT and SRTM DEMs (Gardelle et al., 2013).However, Gardner et al. (2013) and Kääb et al. (2015) found presumably negative mass budgets in Pamir using ICESat laser altimetry data from 2003 to 2008/09.
In this study, we investigate four decades of glacier variations at Muztagh Ata (38 • 17 N, 75 • 07 E; 7546 m a.s.l., Fig. 1), situated in the most eastern Pamir, west of the Taklamakan Desert (Xinjiang province, China).This massif  Tian et al., 2006;Duan et al., 2007). of anomalously high topography reaches ∼ 1500 m higher than most neighbouring peaks in the eastern Pamir (Seong et al., 2009b;Yang et al., 2013).Its glaciers are seasonal to long-term water resources and play an important regulating role for downstream freshwater supply.Moreover, they act as valuable indicators of a changing climate (Seong et al., 2009b;Vaughan et al., 2013;Qiu, 2014).The cold valley glaciers are of continental type (Shih et al., 1980;Zhou et al., 2014).Analysis of gridded precipitation data and pre-cipitation measurements since July 2010 at elevations up to 4900 m a.s.l.indicate that the main accumulation seasons are spring and summer (unpublished data; Maussion et al. 2014).A roughly north-south-trending high ridge and watershed divides the massif into a windward area with small valley glaciers exposed towards the westerlies and an eastern leeward part with steeper slopes.Glacier meltwater drains southwards to the Taxkorgan River, a tributary of the Yarkand River, and northwards to the Gezhe River, being a tributary of the Kaxgar River.With a length of ∼ 18 km and an extent of 86.5 km 2 the debris-covered Kekesayi Glacier is by far the largest glacier of this massif (Shangguan et al., 2006;Seong et al., 2009b, a;Yang et al., 2013).Muztagh Ata is situated in one of the driest glacierized areas of China and one of the coldest environments in these low-and mid-latitude regions (Shangguan et al., 2006;Zhou et al., 2014).Its cold and semi-arid continental type climate is principally influenced by mid-latitude westerlies (Peel et al., 2007;Seong et al., 2009a;Yao et al., 2012).The Taxkorgan meteorological station (37 • 46 N, 75 • 14 E; 3091 m a.s.l.), situated ∼ 50 km south of Muztagh Ata, is the only station on the eastern Pamir Plateau above 3000 m a.s.l.(Shangguan et al., 2006).From 1957 to 2010 the mean annual temperature at this station was measured to be at +3.4 • C, the mean summer temperature at +15.1 • C (June-August), and the mean annual precipitation at ∼ 70 mm (Yan et al., 2013b;Yang et al., 2013).Summer precipitation at the glacier accumulation zone at Muztagh Ata is estimated to account for 30 % of the total annual amount, which was measured to be ∼ 300 mm at 5910 m a.s.l.(38 • 42 N, 75 • 01 E) (Seong et al., 2009b, a).
We used historical and state-of-the-art remote sensing data sets such as Hexagon KH-9 and Pléiades for a detailed reconstruction from 1973 to 2013.The main aim of this research is to fill a knowledge gap, since detailed glaciological studies in this region are scarce and show ambiguous results (Unger-Shayesteh et al., 2013).A further aim is to improve the knowledge of the reaction of debris-covered glaciers in this region by taking Kekesayi Glacier as an example.Therefore surface velocities were measured by TerraSAR-X amplitude tracking and compared to surface elevation changes.

Data
Imagery was acquired under the premise of having a minimum cloud and snow cover.Data are referenced to WGS-84 at UTM zone 43N and to the EGM-96 geoid.Stereo imagery employed for DEM extraction is summarized in Table 1 and described below.

Hexagon KH-9
Hexagon KH-9 was a photographic satellite surveillance system flown during 20 missions from June 1971 to April 1986 by the United States (mission numbers: 1201-1220).During 12 missions (1205-1216) approximately 29 000 photographs were acquired with its mapping camera and declassified in 2002 (Burnett, 2012).It is assumed that for the KH-9 mapping camera a similar design like for the NASA Large Format Camera (LFC) of 1984 was used (Mollberg, 1981;Surazakov and Aizen, 2010).This is a 23 × 46 cm frame format camera with 30.5 cm focal length.Photographs contain four fiducial marks with 1058 réseau crosses and provide ground cover-age of 250 × 125 km at 6-9 m resolution.Imagery was returned in single buckets of films from 171 km operational altitude.DEM extraction is made possible from triplet stereocoverage with 70 % overlap and a base-to-height (b/ h) ratio of 0.4 (Surazakov and Aizen, 2010;NRO, 2011;Burnett, 2012).
KH-9 photographs were scanned in two segments at 7 microns (3600 dpi) with about 1 cm of overlap and stored in an 8-bit TIFF file format.Four overlapping black and white scenes without any geo-corrections were purchased from the United States Geological Survey (USGS).Imagery with frame numbers 16-19 was recorded on 4 August 1973 during mission 1206-5.

Pléiades
Pléiades is a high-resolution satellite system developed by France.Pléiades 1A was launched in December 2011, followed by Pléiades 1B in December 2012.The spatial resolution of the panchromatic channel is resampled to 0.5 m with a pixel depth of 12 bits at acquisition.It provides an image swath of 20 km at nadir by flying at an operational altitude of 694 km.Pléiades offers in-track standard as well as tri-stereo capability with an additional quasi-vertical image.The location accuracy was measured to be 8.5 m for nadir-looking images of Pléiades 1A and 4.5 m of Pléiades 1B (Astrium, 2012;Gleyzes et al., 2012;Berthier et al., 2014).
We purchased a level-1A stereo bundle mosaic of three images that cover an area of 663 km 2 .The image file recorded on 20 June 2013 covers about two-thirds of the eastern part.The image of 19 June 2013 covers the western third, except for some small cloud-covered areas in the south-west that are covered by the image of 3 August 2013.Imagery was provided in GeoTIFF file format with b/ h ratios ranging from 0.2 to 0.3.

ALOS-PRISM
The Japanese satellite system ALOS (Advanced Land Observing Satellite) operated from January 2006 to April 2011 at an operational altitude of 692 km.Its PRISM (Panchromatic Remote-Sensing Instrument for Stereo Mapping) optical sensor consisted of three independent high-resolution panchromatic radiometers.These provided in-track tripletcoverage from backward-, nadir-and forward-looking directions.The b/ h-ratio is up to 1.0 when using forward-and backward-looking views with an inclination of ±23.8 • from nadir.It offered 2.5 m spatial and 8-bit radiometric resolution with a swath width of 35 km in triplet mode.The location accuracy amounts to be 8.1 m for nadir-looking images and 9.3 m for forward and backward views (Takaku et al., 2007;JAXA, 2008;Tadono, 2009).
We purchased a radiometrically calibrated triplet mode scene at level-1B1 in CEOS file format.

Landsat
Glacier delineation for the year 2000 and horizontal measurements of ground control points (GCPs) is based on imagery of the Landsat Enhanced Thematic Mapper 7 (ETM+).Orthorectification to 15 m pixel size for the panchromatic and 30 m for the multispectral band was conducted by use of the Global Land Surveys (GLS2000) data set.Terrain-corrected imagery at level-T1 was provided by USGS, acquired on 11 September 2000.

TerraSAR-X
TerraSAR-X is a German X-band radar satellite launched in June 2007.Data are available in Spotlight, Stripmap, and ScanSAR modes at all achievable polarization arrangements (HH, HV, VV, and VH).We obtained two Stripmap mode images acquired on 10 August and 1 September 2011 during the descending pass of the satellite.Image extent is 19.7 × 21.2 km on ground with a pixel spacing of 0.9 m in slant range (signal direction) and 3.0 m in azimuth (flight direction).The incidence angle at the scene centre is 44.3 • .The data were delivered by the German Aerospace Center (DLR) in single-look complex (SLC) format (Herrmann and Bottero, 2007;Eineder et al., 2008).
3 Data processing

KH-9 image preprocessing
We resampled the KH-9 photograph segments to 14 microns for simpler data handling, in consideration of the large file size.Prior to DEM extraction, it was necessary to reconstruct the original conditions of image geometry at the time of film exposure.This is indispensable to obtain accurate elevation information from KH-9 stereo photographs.Film distortions evolved over time due to duplication and storage during almost four decades.Such distortions were corrected by evaluating its réseau grid overlaid in the photograph which consists of 1058 crosses at 1 cm spacing.The original image geometry was reconstructed by a second-order bilinear interpolation.Based on a Python tool developed by Pieczonka et al. (2013), we automatically determined all réseau crosses in the imagery and resampled them back to their initial reference positions.Réseau crosses were expected to later confuse terrain extraction and were therefore eliminated using bicubic interpolation from surrounding pixels (Pieczonka et al., 2013).Prior to mosaicking, Wallis filtering with 51 × 51 pixels window size and histogram equalization was finally conducted for contrast enhancement (Surazakov and Aizen, 2010;Pieczonka et al., 2013).In most scanned photographs, there are unfortunately no fiducial marks.Hence, we assumed the image principle point to be identical with the central réseau grid coordinate of both corresponding mosaicked segments.This position was also considered as origin of initial reference for image geometry reconstruction and is therefore not affected by resampling.Besides the film distortion, also a rotation component appears around the principle point in the distortion vectors.This probably originates from an occasionally slightly rotated scan of a film segment (Fig. 2) (Holzer et al., 2012).

DEM extraction
Measuring ground control points (GCPs) proved to be challenging due to the remoteness of the region and the lack of accurate ground truth data.GCPs were situated at stable and plain terrain, ideally close to laser altimetry measurements of the Ice, Cloud, and Elevation Satellite (ICESat) which proved to be a valid elevation source (Nuth and Kääb, 2011).
The SRTM-3 DEM served as z-reference for one-third of the GCPs, since no ICESat information was available.x and y coordinates were measured from a pan-sharpened Landsat 7 ETM+ scene dating from 2000.All GCPs were finally crosschecked in Google Earth TM .Finding suitable GCPs was particularly difficult for Hexagon KH-9 due to its long temporal baseline when anthropogenic objects like road intersections and houses did not exist back in 1973.We measured 18 GCPs for KH-9 Hexagon, ALOS-PRISM is covered by six GCPs, and the Pléiades mosaic by 11 GCPs with at least four GCPs per scene.DEM extraction from Hexagon KH-9 photographs is based on a non-metric frame camera model using the Leica Photogrammetry Suite 9.2 (LPS).Inner orientation settings with 30.5 cm focal length were defined as fix for triangulation, but we used Brown's physical model to compensate for unknown lens and potentially remaining film distortions.The principle point offset was determined from the central réseau-cross coordinate to the midpoint of the image, which is defined by its extent.Due to the lack of ephemeral or analogue metadata information, the exterior orientation is solely based on 18 measured GCPs, by taking the earth curvature into account.The RMSE of bundle block adjustment proved to be 0.49 pixels.The DEM was extracted with adaptive au- DEM extraction from ALOS-PRISM and Pléiades stereo imagery was performed using the OrthoEngine of the PCI Geomatica 2013 software package.We used its Rational Functions model to derive DEMs with first-order RPC adjustment from attached ephemeral data and the measured GCPs.The GCP residuals of bundle block adjustment proved to be 0.36 m in x and 0.34 m in y direction for ALOS-PRISM as well as 0.18 and 0.12 m for Pléiades, respectively.Wallis filtering was applied to improve the image matching process.The DEM of ALOS-PRISM at 10 m resolution is based on epipolar pairs from the backward-and nadir-looking as well as from the nadir-and forward-looking views with their highest obtained score.For each of the three Pléiades scenes we derived a very high-resolution DEM at 1 m resolution from their stereo views (Fig. 4).The DEM of ALOS-PRISM does not cover the westernmost part of Muztagh Ata, and the DEM of Pléiades is affected by several gaps due to clouds in the south-west.

Clean-up of DEM areas with low quality
All DEMs including SRTM-3 were resampled to 30 m resolution by cubic convolution and to a common raster grid extent for exact cell alignment.PCI OrthoEngine provides an additional image which represents the stereo matching score for each extracted DEM pixel.We applied a threshold of 0.7 to exclude elevations of poor accuracy in the DEMs of ALOS-PRISM and Pléiades.The thematic point status image of LPS showed that most calculated DEM points from KH-9 Hexagon (76 %) were of fair quality with correlation scores ranging from 0.5 to 0.7.Of these points, 17 % proved to be of good accuracy and 7 % of excellent accuracy, with coefficients higher than 0.85.Large DEM parts, however, consisted of clearly erroneous elevations despite a quality indicated as fairly good.This was particularly observed in the case of poor contrast in the KH-9 imagery.We identified such elevations by both their hillshade and their difference image relative to SRTM-3, and set them to "no data".

DEM co-registration
Horizontal DEM co-registration to SRTM-3 was conducted analytically by minimizing the elevation error based on the relationship between elevation difference and aspect (Nuth and Kääb, 2011).The approach was based on a routine implemented in Python programming language by Pieczonka et al. (2013).For each DEM we calculated its difference image relative to SRTM-3 by excluding non-stable terrain such as (rock) glaciers, ice-cored moraines, and lakes.To allow for the slope dependency of the method, we excluded all terrain below a slope of 10 • .The initial spatial resolution of SRTM-3 (90 m) is coarser than that of the derived DEMs.This can lead to resolution-implicated biases at topographic extremes where curvature is strong (Paul, 2008;Gardelle et al., 2012a).To make considerations for outliers and such curvature effects, we first bounded valid pixels of DEM differencing to their 5 % and 95 % quantiles (Hoffmann and Walter, 2006).Subsequently we excluded all elevation differences outside of its two-tailed 1.5 times interquartile range (Pieczonka et al., 2013).The determined horizontal shifts were iteratively reduced until an accuracy of at least 1 m in x and y direction in respect to SRTM-3 was reached for each DEM.
Spatially varying elevation biases were corrected by twodimensional trend surfaces in off-glacier regions (Bolch et al., 2008;Pieczonka et al., 2011).These were calculated from DEM difference images and reduced the mean heightoffset relative to SRTM-3 on stable terrain to zero.Quantile analysis was employed in a similar way as for horizontal alignment to take curvature effects into account.

SAR image co-registration
TerraSAR-X data were processed using GAMMA Remote Sensing software.Data were converted to software-readable SLC format which contains the amplitude (backscatter) and phase information from the signal interaction with the Earth's surface.We defined the image of 1 September 2011 as "slave" and co-registered it at sub-pixel accuracy to the "master" image of 10 August 2011.Orbital offsets and ionospheric shifts were corrected by means of amplitude correlation, determined over well-distributed windows of 128 × 128 pixels.In selected windows the bilinear offset polynomial for image registration was calculated by employing a threshold of 4.0 from the signal-to-noise ratio (Strozzi et al., 2002).

Geodetic glacier mass balances
Geodetic glacier mass balances are based on maps of elevation differences ( h) by subtracting older date elevations (e.g.KH-9 Hexagon) from newer ones (e.g.Pléiades).Such difference images were generated for all possible DEM combinations of KH-9 Hexagon (1973), SRTM-3 (1999), ALOS-PRISM ( 2009), and Pléiades (2013).SRTM voids are particularly observed at steep slopes and mountain ridges, while most of the glacier areas consist of non-interpolated data.We restricted mass-balance calculations to the original SRTM-3 surfaces and excluded gap-filled voids because of high inaccuracies (Kääb et al., 2012).On stable terrain, slight offsets in mean height were induced while differencing DEMs which were both extracted from optical data.The biases were below 1 m and resulted from co-registration of all DEMs exclusively to SRTM-3.Offsets were corrected to keep offglacier elevation differences for all DEM combinations in their mean at zero.We selected 13 larger glaciers of different orientations to calculate their elevation and volume change as well as their mass balance by assuming an ice density of 850 ± 60 kg m −3 (Huss, 2013).Mass change was estimated for the entire glacierized area of Muztagh Ata as well by also taking the mass balances of individual glaciers into account.
The glacier size was defined by the largest extent of the correspondent mass-balance investigation period.

Outlier detection and gap filling
Data gaps smaller than 0.01 km 2 were closed by a mean filter based on surrounding h values.Outlier detection and gap filling of remaining h voids in glacier areas were employed separately for each glacier accumulation and ablation zone.These were separated by equilibrium line altitudes (ELAs), based on estimations of the first Chinese glacier inventory (Shi et al., 2008).ELAs were cross-checked in ALOS-PRISM and Pléiades ortho images and adapted if necessary (see Table 4).ELAs are also based on geometric calculations of Seong et al. (2009b) and in situ measurements at Muztagh Ata Glacier No. 15 (Yao et al., 2012).
For each 25 m elevation band in the ablation zone, we restricted the minimal and maximal allowable elevation differences to its 5 and 95 % quantiles and replaced outliers with its marginal quantile values.In the case of Kematulejia and Kuosikulake Glacier, this restriction was tightened to the 31.7 and 68.3 % quantiles because of higher noise.Remaining no-data gaps were filled by mean elevation differences calculated for each 25 m elevation section in the ablation zone.Poor image contrast in the snow-covered accumulation zone led to high noise of h values and large areas without elevation estimates.Since no plausible statistical replacement values could be derived, we set missing h pixels in the accumulation zone to zero by assuming only minor elevation changes for these areas (Schwitter and Raymond, 1993).Elevations outside the range of the 31.7 and 68.3 % h quantiles of each glacier accumulation zone were considered as outliers and also set to zero.Statistical outlier detection and gap filling employed to individual glaciers were not possible for the entire glacierized area due to diverging glacier elevation changes at similar altitudes.For the remaining glacierized area, we subsequently defined a h threshold of ±100 m for the ablation area and set all h pixels of the accumulation area to zero.

SRTM-3 C-band radar penetration
SRTM-3 C-band penetrations strongly depend on the condition of the upper snow/firn/ice layers below the glacier surface.Landsat 7 ETM+ imagery recorded on 7 February 2000 shows slight snow coverage with mostly snow-free glacier tongues at Muztagh Ata.In this study, we assume that SRTM-3 approximately detects the surface at the end of the melt season in 1999.C-band penetrations were corrected separately for glacier accumulation and ablation zones based on estimates by Kääb et al. (2012).In doing so, we averaged penetration depth estimates of the three nearby and southwards-situated Hindu Kush, Karakoram, and Jammu-Kashmir study sites.This results in penetration assumptions of 4.3 ± 0.9 m for firn and snow (accumulation zone) and 1.5 ± 0.9 m for clean-ice ablation zones, by assuming no penetration in the case of supraglacial debris.Added offsets are positive if SRTM is representing glacier surfaces of older date (1999-2009/2013) and negative if it is of newer date .There is no need for such corrections if the DEM differencing is solely based on optical data.

Glacier surface velocities
Surface velocities of Kekesayi Glacier were determined by amplitude tracking that also performed well in most parts of the upper glacier area (e.g.Strozzi et al., 2002;Floricioiu et al., 2010;Rankl et al., 2014).It was not possible to retain the interferometric phase due to temporal decorrelation over 22 days.Phase-based methods such as DInSAR (Goldstein et al., 1993), GInSAR (Sharov et al., 2002) or double difference InSAR (Floricioiu et al., 2010) subsequently failed due to low coherence.The normalized cross-correlation function was estimated in the co-registered "master-slave" images using 64 × 64 pixel windows.Motion in azimuth and range direction was yielded by the peak location of this function.Glacier surface dynamics were determined as where d range and d azimuth are the motions of the glacier surface in range and azimuth directions, respectively.d absolute represents the magnitude of surface velocity and δ flow depicts the direction of glacier flow (Strozzi et al., 2002).

Uncertainties
Mapping precision of clean-ice glaciers can be roughly estimated by a 1-pixel variability of glacier outlines based on the spatial resolution of its reference imagery (Bolch et al., 2010;Frey et al., 2012;Paul et al., 2013).We enlarged buffers to consider the difficult visual interpretation of debris cover and to take the high spatial resolution of some images into account.The glacier reference outlines from Pléiades (2013) were buffered with ±1 m and from ALOS-PRISM ( 2009) with ±2 m.For Hexagon KH-9 and Landsat 7 ETM+ we followed the buffer sizes proposed by Bolch et al. (2010) and used a glacier size variability of ±10 m for 1973 and of ±7.5 m for 2000.Uncertainties of glacier area and length changes are defined by the root sum squares of each error term and dominated by higher mapping inaccuracies of older data sets.We calculated the normalized median absolute deviation (NMAD), the 68.3 and the 95 % quantiles to measure the vertical DEM precision of all difference images from the multi-temporal DEMs (Höhle and Höhle, 2009).Similar to DEM co-registration, this calculation is based on DEM differencing by excluding non-stable terrain and by considering outliers and curvature effects (Table 2).Density of glacier ice is assumed to deviate in the range of ±60 kg m −3 (Huss, 2013).Another influence onto DEM differencing with SRTM-3 is its penetration-depth uncertainty.This was estimated to be ±0.9 m as the highest uncertainty of the averaged penetration depth corrections of Kääb et al. (2012).The final mass-balance and volume-change uncertainties are the root of the sum of each squared error term and consist of the NMAD as well as the uncertainties of ice-density assumption and of C-band penetration depth correction if applicable.For annual mass-balance rates this is converted into water equivalent and divided by the observational years.The uncertainties of volume change are multiplied by the glacier area and converted to ice equivalent.The uncertainty in surface velocities results from the imprecise matching of the glacier surface features within the search windows.We measured residual velocities at a stable and plain surface below the glacier terminus, where the channels carry the water discharge from the glacier.The RMSE was estimated over non-moving terrain of ∼ 5 km 2 to be ±0.58cm per day.

Results
Investigated glaciers were named according to their ID in the GLIMS Database (GLIMS and NSIDC, 2005, updated 2014) (see Fig. 1).

Glacier area and length changes
The glaciers at Muztagh Ata showed heterogeneous variations with some fluctuating (or advancing), but mostly stable or continuously retreating glacier tongues during the last four decades.Area and length changes are highly variable from one glacier to another, even if they are located adjacently.Several glaciers such as Kekesayi (G075225E38255N) are heavily covered by debris and did not indicate any change at their frontal position.Average glacier retreat was observed to be −1.0 ± 0.3 m a −1 from 1973 to 2013.Overall glacier length changes are becoming less negative and show even positive values for later periods (Table 3).The determined overall area change of −0.6 ± 3.9 % (−0.02 ± 0.1 % a −1 ) is low and not significant.This corresponds to a glacier area reduction of −1.6 ± 10.6 km 2 from 274.3 ± 10.6 km 2 in 1973 to 272.7 ± 1.0 km 2 in 2013 (Table 4).The maximal extent of glaciation was observed to be at ∼ 5000 m (Fig. 3a).

Geodetic glacier mass balances
Glacier elevation change is determined by difference images from DEMs of Pléiades and ALOS-PRISM to SRTM-3 for 1999 to 2009/2013 (Fig. 8), Pléiades and ALOS-PRISM to KH-9 Hexagon for 1973 to 2009/2013 (Fig. 6) as well as SRTM-3 to KH-9 Hexagon for 1973 to 1999 (Fig. 7a).Difference images of multiple time periods show clear temporal variations of ice thickness change and movement, which is particularly evident for advancing or fluctuating glaciers.Kuokuosele Glacier (G075156E38175N) showed a strong mass gain at its downstream part from 1973 to 1999, which, however, was limited to its continuously advancing tongue after 1999.Glacier G075075E38189N revealed down-wasting at its retreating tongue from 1973 to 1999, while surface elevation gain was observed in its middle part.This led to subsequent glacier advance with mass gain at its toe and loss in its middle part.Despite its more or less stable tongue, down-wasting was observed for Kuosikulake (G075092E38214N) Glacier from 1973 to 1999.
Following mass gain at its lower part might explain the sud-den advance after 2009.Clear down-wasting despite stable frontal positions was observed for debris-covered Kekesayi (G075225E38255N) and G075171E38163N glaciers during the entire study time period.
Average mass budgets at Muztagh Ata in the range of −0.03 ± 0.33 m w.e. a −1 (1973-2009) to −0.01 ± 0.30 m w.e. a −1 (1973-2013) have been nearly balanced for more than 40 years.For different periods of the investigated time span, however, mass changes strongly vary from one glacier to another.Kekesayi (G075225E38255N), the largest glacier of the Muztagh Ata Massif, shows ice mass loss during all investigated time periods.There are indications that most glaciers had more positive budgets in the last decade as compared to the period before 1999 (Tables A1).The time period between the ALOS-PRISM and Pléiades data takes is only 4 years and should be considered as too short for reliable results.Its difference image (Fig. 7b), however, shows mostly low noise, and the characteristics of surface elevation changes continue well in line with our results from other periods.ELAs were adjusted from the first Chinese glacier inventory (Shi et al., 2008) and other sources.

Glacier surface velocities of Kekesayi Glacier
Surface velocities of Kekesayi (G075225E38255N) Glacier reached up to 20 cm per day in August 2011.This corresponds to a maximal flow of ∼ 70 m per year if a similar flow throughout the year would be assumed (Fig. 5).Ice flow at more than 15 cm per day (∼ 55 m a −1 ) is maximal at its middle part, downstream of the confluence of the tributaries T2 and T3.Surface movements at the glacier margin, independent of the location, are slow due to retarding friction.Surface velocities slow down consistently with the glacier stream and become almost insignificant where stronger surface lowering occurs.Hence, we conclude that the glacier is stagnant to about 3 km upstream from the terminus.

Discussion
7.1 Glacier area and length changes Yao et al. (2012) found in the eastern Pamir the least glacier shrinkage (−0.07 % a −1 ) and retreat (−0.9 m a −1 ) compared to the Tibetan Plateau and the Himalaya.More than 60 surging glaciers were identified in the central Pamir for the time period from 1972 to 2006 (Kotlyakov et al., 2008).Contrary to this trend high shrinkage was observed in the Zulumart range south of Pamir Alay, where glaciers shrank −7.8 % (−0.65 % a −1 ) from 1978 to 1990 which accelerated to −11.6 % (−1.05 % a −1 ) until 2001 (Khromova et al., 2006).Shrinkage was also reported by Shangguan et al. (2006) at the Muztagh Ata and Kongur massifs and was measured to be at −7.9 % (−0.21 % a −1 ) from 1962 to 1999.This was determined to come along with increasing retreat from −6.0 m a −1 (1962/1966-1990) to −11.2 m a −1 after 1990 (Shangguan et al., 2006).Our determined shrinkage (−0.02 ± 0.1 % a −1 ) and retreat (−1.0 ± 0.3 m a −1 ) is not as negative compared to several rates calculated per-glacier.We attribute the more negative rates of the aforementioned studies to the different sites and investigation periods.However, the differences can also stem from uncertain glacier boundaries in the Chinese topographic maps (Bolch et al., 2010) and as a result of the more difficult glacier interpretation in Landsat imagery with a coarser resolution.In total, we would also expect less glacier shrinkage and retreat at Muztagh Ata as in other areas of the eastern Pamir study region of Yao et al. (2012) by reason of, on average, nearly balanced observed mass budgets in this study.

Glacier mass balances
DEM differencing of multiple time periods confirms spatially as well as temporally inhomogeneous glacier mass changes at Muztagh Ata, but on average nearly balanced budgets.These were determined to be −0.the last decade.This deviation may be attributed to the different time periods of the studies and the uncertain penetration of SRTMs C-band radar into ice and snow.
An overall mass loss in the western and central Pamir seems to be more likely when considering the measured continuous glacier shrinkage (Khromova et al., 2006) as well as the negative mass budget of Abramov Glacier in Pamir Alay (measured years 1968(measured years -1997(measured years and 2011(measured years /2012(measured years , WGMS, 2013) ) and the volume loss of Fedchenko Glacier, the by far largest and debris-covered glacier in the central Pamir (Lambrecht et al., 2014).A region of positive anomaly seems to start in the Karakoram (Hewitt, 2005;Gardelle et al., 2012b) and continues over the eastern Pamir (Yao et al., 2012; this study) to western Kunlun (Gardner et al., 2013;Neckel et al., 2014;Kääb et al., 2015) and central Tibet (Neckel et al., 2014).

Down-wasting, surface dynamics and area changes of debris-covered glaciers
Glacier tongues at Muztagh Ata which reach below 4700 m a.s.l. are usually covered by debris, with increasing debris thicknesses of up to several metres at lower altitudes (Yang et al., 2013).Most of these glaciers do not show visual indications of retreat, and Shangguan et al. (2006) could not detect significant area changes at ∼ 90 glaciers at their Muztagh Ata and Kongur study site, possibly due to debris cover.However, our results of DEM differencing exhibit clear surface lowering at the downstream debris-covered glacier parts.This demonstrates that glaciers may have negative mass balances despite thick debris cover and stable terminus positions.Decoupling of area from volume loss can be provoked by supraglacial debris, which can reduce glacier melt rates if debris coverage is exceeding a few centimetres of thickness.Stagnant debris-covered terminus positions must, hence, not indicate balanced glacier conditions (Bolch et al., 2011;Scherler et al., 2011;Lambrecht et al., 2014;Pellicciotti et al., 2015).In this regard, Fedchenko, as the Pamir's by far largest glacier, lost more than −5 km 3 of volume during the last eight decades (∼ −6.0 %), but it shrank by only −1.4 % at its debris-covered tongue (Lambrecht et al., 2014).Similar results were found by Pieczonka and Bolch (2015) for the central Tian Shan.
The largest glacier at Muztagh Ata, Kekesayi (G075225E38255N), appeared, by visual indication, to be stagnant from 1973 to 2013.DEM differencing, however, clearly indicates increasing ice mass loss at its downstream part during all investigated study periods.Surface lowering at the heavily debris-covered tongue reached up to 40 m in sum for the last four decadal measurements.Down-wasting becomes highest where surface velocities decrease to almost insignificant values, particularly about 3 km upward from its terminus.A profile along the central glacier flow line supports an obvious relationship between surface velocity and down-wasting (Fig. 5).This was previously identified with similar methods by Pellicciotti et al. (2015) for debriscovered glaciers in central Himalaya.Yang et al. (2013) set a polynomial fit through multi-annual surface movements of Kekesayi Glacier, which were measured between 1998 and 2010 from Landsat imagery.The average upstream velocity of up to 50 m per year (∼ 14 cm per day) is in the range of our measurements, while Zhou et al. (2014) presents winter velocities that did not exceed ∼ 11 m per year from 2008 to 2010.Glacier flow in 2009 shows lower rates in winter months (∼ 9 cm per day) as compared to summer rates at ∼ 15 cm per day (Yan et al., 2013a).These studies confirm seasonal and annual glacier flow variability at the central part of Kekesayi Glacier, with little or no fluctuations at the terminus.Its tongue is widely covered by supraglacial ponds that absorb large amounts of energy and thus contribute to down-wasting.The insulation effect of thick debris coverage, however, causes such glaciers to melt at lower rates, which might indicate retarded climate response.Down-wasting associated with negligible or little retreat in case of debris cover is also confirmed by studies of Bolch et al. (2008Bolch et al. ( , 2011) ) and Pellicciotti et al. (2015).This underlines the importance of volume change investigations as more reliable indicators for climate-related glacier responses.

Glacier response to climate change
From the late glacial period to the Holocene, the glaciation style of the massif has changed from an expanded ice cap to deeply entrenched valley and cirque glaciers (Seong et al., 2009b, a).Recent glacier variations might be a response to changing atmospheric circulation patterns.Yao et al. (2012) observed strengthening westerlies coming along with an increase of precipitation, whereas the Indian monsoon is quite likely to weaken.Both temperature and precipitation in the study region have been increasing (Shi et al., 2007;Qiu, 2014).In Xinjiang province (north-west China), from 1961 to 2008 both mean annual temperature and precipitation increased per decade by +0.3 • C and +7.4 mm (Zhang et al., 2012).The summer temperature (June-August) at the close by Taxkorgan meteorological station rose by +0.7 • C from 1957 and 2000, while annual precipitation slightly increased at the same time (Shangguan et al., 2006;Tian et al., 2006;Yao et al., 2012).In summer 2003, a 41.6 m deep ice core was drilled at 7010 m a.s.l. at Muztagh Ata (38 • 17 N, 75 • 06 E, see Fig. 1) (Tian et al., 2006;Duan et al., 2007).Its isotope variations were found to be in good agreement with annual air temperature changes measured at Taxkorgan.However, starting in the 1990s, a more rapid warming trend of +2.0 to +2.4 • C per decade was observed, compared to Taxkorgan station measures at + 0.18 • C per decade (Tian et al., 2006).Reconstructed mass-balance rates of the glaciers in the Muztagh Ata region do not agree with our results and show much higher wastage after 1990 (−0.42 m w.e. a −1 ) as compared to  the determined mean at −0.12 m w.e. a −1 for 1960 to 2003 (Duan et al., 2007).
Accordingly, Khromova et al. (2006) and Shangguan et al. (2006) assume that further glacier shrinkage might have been accelerated by rising summer temperatures, particularly since the 1990s.Ablation is reported from June to August because of positive expected mean summer air temperatures beyond the glacier terminus (Shangguan et al., 2006;Yang et al., 2013).It is, however, suggested that glaciers in this region are more sensitive to a change in precipitation than to temperature (Seong et al., 2009a, b).Glaciers at Muztagh Ata are situated at relatively high altitudes, where despite warming the air temperature still remains far below freezing during winter.Increasing precipitation from strengthening westerlies can, hence, lead to higher snow accumulation.This relativizes the negative effects of climate change regarding warming and might be one of the reasons why average shrinkage and ice mass loss at Muztagh Ata is low and insignificant.Under current climate conditions, and by reason of increasing precipitation, Yao et al. (2012) would expect an advance of glaciers in the eastern Pamir.The observed advance in this study might also be a response to three cooling periods with increasing annual precipitation measured from 1961-1968, 1973-1977, and 1985-1993 at Taxkorgan station (Shangguan et al., 2006).

Uncertainties of geodetic mass balances from optical data
Low contrast alterations and oversaturation hampers terrain extraction from optical stereo imagery, particularly at snowcovered accumulation zones.Substitution of low-quality h values by zero in these zones is a consequence of lacking statistical alternatives.Potentially induced biases in volume change are therefore difficult to quantify and would be rather speculative.The impact in this study is less critical for Pléiades as compared to ALOS-PRISM (2.5 m, 8-bit), since its higher geometric (0.5 m) and radiometric (12-bit) resolution led to a lower rate of poor elevation estimates.KH-9 Hexagon (6-9 m, 8-bit), moreover, shows high noise in its DEM at low-contrast terrain.Debris-covered and crevassed glacier surfaces, however, are of much better quality, and to that effect, the relatively high NMAD of the KH-9 DEM is possibly overestimated.Seasonal variations in glacier elevation change are taken into account by the conservative NMAD uncertainty estimate.Multi-temporal results prove to be in line despite such uncertainties.Median values close to zero (the mean is by construction zero) give confidence of a safe, almost Gaussian distribution (Table 2).The vertical precision in this study is in the range of the SRTM-3 accuracy specifications.These are stated to be ±6 m relative and ±16 m absolute (Rabus et al., 2003).We would expect a higher precision in case of a more accurate reference than the SRTM-3 DEM.This assumption is supported by similar NMAD values from optical stereo data in the study of Pieczonka et al. (2013), and by a much lower NMAD of 2.5 m from high-resolution DEM differencing of ALOS-PRISM to Pléiades in this study (Fig. 7b).

SRTM C-band penetration depth correction
Specific C-band penetrations into snow and ice must be corrected for SRTM due to different weather and the subsequent snow-cover conditions during the acquisition in February 2000.The C-band radar waves penetrate into clean glacier ice and particularly through newly fallen layers of fresh snow.Penetrations reach up to 10 m in dry cold firn and 2 m in exposed ice (Rignot et al., 2001).Gardelle et al. (2013) measured mean penetrations of 1.8 ± 1.5 m in glaciers of the Pamir, but admits that this value might be underestiwww.the-cryosphere.net/9/2071/2015/The Cryosphere, 9, 2071-2088, 2015 mated.This is supported by Kääb et al. (2015) who found larger SRTM C-band penetrations of 5-6 m in the Pamir.We therefore referred to larger estimates which were determined for the three nearby Hindu Kush, Karakoram, and Jammu-Kashmir study sites of Kääb et al. (2012).Its westerly influenced glaciers are situated further south but at about the same longitude, and we, hence, suppose similar penetrations for Muztagh Ata despite a higher degree of continentality.Penetration depths in these regions are 5.1 ± 0.7 m, 5.5 ± 0.3 m, and 2.3 ± 0.9 m for firn and snow, as well as 1.7 ± 0.6 m, 1.1 ± 0.5 m, and 1.7 ± 0.7 m for clean ice, respectively (Kääb et al., 2012).We averaged these estimates in consideration of their wide geographic spreading.Slightly negative mass budgets observed from 1973 to 1999, compared to more positive values after 1999, might still indicate underestimated corrections.Possibly biased trends prior and after 1999 are, however, insignificantly low and the derived mass-balance results are well in line.DEM differencing solely based on optical data is not subject to such eventual biases.

Conclusions
Glaciers at Muztagh Ata, situated in eastern Pamir, are of high importance for seasonal freshwater supply and act as a valuable climate indicator.Detailed knowledge of glacier changes in this remote and high mountain region is, however, scarce.We used remote sensing data sets of Hexagon KH-9 (1973), ALOS-PRISM (2009), Pléiades (2013), and Landsat 7 ETM+ in conjunction with SRTM-3 (2000) to investigate four decades of glacier variations at Muztagh Ata.These are heterogeneous and differ spatially as well as temporally.
The debris-covered glaciers show no or only negligible visual changes at their frontal position.Differencing of multitemporal Digital Elevation Models (DEMs), however, reveals clear down-wasting at their tongues, despite mostly thick debris coverage.Some south-west-exposed glacier tongues fluctuated or advanced, with infrequent variations in ice thickness.The total glacier shrinkage of −0.02 ± 0.1 % a −1 , from 274.3 ± 10.6 km 2 in 1973 to 272.7 ± 1.0 km 2 in 2013, is low and not significant.Averaged mass budgets based on geodetic measurements are slightly but insignificantly negative before 1999 (−0.04 ± 0.42 m w.e. a −1 ) and positive afterwards (+0.04 ± 0.27 m w.e. a −1 ).This might still result from a potentially underestimated SRTM-3 C-band penetration into snow and ice.Slightly positive observed budgets after 1999, however, could possibly reflect a regional-wide positive anomaly with increasing snow accumulation from strengthening westerlies.Mass gain for glacier G075058E38248N (so called Muztagh Ata Glacier No. Our study emphasizes the importance of volume change investigations, particularly for debris-covered glaciers.Largely untouched stereo photographs of the US Keyhole-9 spy program provide long-term information of historic glacier situations and were highly valuable for volume change investigations prior to the acquisition date of SRTM.This study presents, in combination with the recently recorded high-resolution Pléiades imagery, the longest time series of geodetic mass balances for the eastern Pamir at the glacier scale.

Table 1 .
Overview of optical stereo imagery used for DEM extraction and subsequent geodetic mass-balance measurements.Optical sensor (stereo) Acquisition date Stereo mode (b/ h-ratio) Spatial/radiometric res.

Figure 2 .
Figure 2. Distortion vectors of réseau crosses from their measured to their initial reference positions in a KH-9 photograph segment covering Muztagh Ata (frame 17a of mission 1206-5).
The vertical accuracy of SRTM-3 decreases in case of steep terrain, and we thus only considered flat areas up to a slope angle of 10 •(Falorni et al., 2005).Offsets were usually apparent as tilts and therefore corrected by linear surfaces.The DEM of ALOS-PRISM indicated a slight second-order polynomial trend offset.All three extracted DEMs of Pléiades were mosaicked to one single file after co-registration.
Assessment of glacier variations4.1 Glacier area and length changes KH-9 Hexagon, ALOS-PRISM, and Pléiades imagery were orthorectified by use of DEMs generated from their own stereo data.Glaciers and ice divides were manually mapped using the very high-resolution (0.5 m) Pléiades ortho-image mosaic and finally cross-checked with Google Earth TM .Glacier mapping, particularly in case of debris coverage, was also based on a hillshade from the Pléiades DEM at 1 m resolution and derived morphometric parameters.The generated inventory representing the glacier situation in 2013 was afterwards manually adjusted to the extents of the years1973,  2000 and 2009.This is based on the orthorectified KH-9 Hexagon (6.0-9.0 m), pan-sharpened Landsat ETM+ (15 m) and ALOS-PRISM (2.5 m) data sets.Changes in glacier length were measured along their central flow line.

Glacier
Figure 3. (a) Hypsometry of the glacierized area at Muztagh Ata, (b) cumulative glaciological mass balance of Muztagh Ata Glacier No. 15 (G075058E38248N) vs. cumulative geodetic mass balance of the entire Muztagh Ata Massif.

Figure 4 .
Figure 4. (a) Hillshade of the Pléiades DEM at 1 m resolution with advancing tongue of Kuokuosele Glacier (G075156E38175N, left) and stable tongue of glacier G075171E38163N (right), (b) compared to 10 m DEM of ALOS-PRISM.

Figure 5 .
Figure 5. Surface velocities and flow directions of Kekesayi (G075225E38255N) Glacier in August 2011.The profile shows the surface velocities and the corresponding down-wasting (1973-2013) along the central glacier flow line, upstream from A to B.

Figure 6 .
Figure 6.Co-registered map of elevation differences during 1973-2013 (a) and 1973-2009 (b) after outlier and gap-filling processing for glacier mass-balance and vertical uncertainty calculation.

Figure 7 .
Figure 7. Co-registered map of elevation differences during 1973-1999 (a) and 2009-2013 (b) after outlier and gap-filling processing for glacier mass-balance and vertical uncertainty calculation.

Figure 8 .
Figure 8. Co-registered map of elevation differences during 1999-2013 (a) and 1999-2009 (b) after outlier and gap-filling processing for glacier mass-balance and vertical uncertainty calculation.
15) is confirmed by in situ measurements for the period 2001-2014.Differencing based on only optical DEMs is not subject to penetration depth uncertainties.Optical approaches indicate nearly balanced budgets for the last four decades (−0.01 ± 0.30 to −0.03 ± 0.33 m w.e. a −1 ).Keke-sayi (G075225E38255N), the largest glacier at Muztagh Ata, shows more negative trends in the range of −0.08 ± 0.30 to −0.11 ± 0.33 m w.e. a −1 from 1973 to 2009/2013.Amplitude tracking of TerraSAR-X data from summer 2011 indicates a stagnant glacier tongue where down-wasting occurs.Upstream velocities fluctuate in its tributaries and are in the range of ∼ 10-15 cm per day.
STD: standard deviation.The mean is equal to zero (RMSE = STD).

Table 3 .
Glacier length changes ( L) at Muztagh Ata for selected glaciers that have mass-balance estimates.

Table 4 .
Equilibrium line altitude (ELA), glacier area (A), and changes ( A) from 1973-2013 at Muztagh Ata for selected glaciers that have mass-balance estimates and for all glaciers of the study site.