Recent mass balance of Purogangri ice cap , central Tibetan Plateau , by means of di ff erential X-band SAR interferometry

Introduction Conclusions References


Introduction
The Tibetan Plateau (TP), also known as the third pole, is characterized by many glaciers and ice caps.Information of whether that ice is losing or gaining mass is a valuable indicator to understand the climate variability on the TP (Yao et al., 2012).This is of high importance as direct measurements of climatic parameters are sparse, especially in the central part of the TP (Thompson et al., 2006).In this study, we take a closer look at the recent mass balance of the Purogangri Ice Cap (PIC), Tibet's largest ice field (Shi et al., 2009).The PIC is located at 33 • 05 N, 89 • 10 E in the continental westerly-dominated north-central part of the TP (Thompson et al., 2006).Due to its remoteness and high altitude of 5800 m a.s.l, on average, field work at the PIC would involve high costs, large logistical efforts and great physical strain.Therefore, remote sensing is a promising alternative to conduct mass balance measurements of this ice cap.
In this study, we estimate the mass balance by geodetic means.In general, the geodetic mass balance is calculated by subtracting data sets of glacier surface elevation acquired at different times (e.g., Rignot et al., 2003;Kääb, 2008;Haug et al., 2009;Bolch et al., 2011).For the PIC two such data sets are available.The ice cap was almost completely mapped by X-band SAR Interferometry (InSAR) in February 2000 during the Shuttle Radar Topography Mission (SRTM) resulting in a Digital Elevation Model (DEM).Almost exactly 12 years later, on January 26, 2012, the ice cap was mapped by single-pass X-band InSAR from TerraSAR-X and its addon for digital elevation measurements (TanDEM-X) (Krieger et al., 2007).
Utilizing these interferometric data sets, we followed two different approaches to estimate the geodetic mass balance of the PIC between 2000 and 2012.First, we employed Differential Synthetic Aperture Radar Interferometry (DInSAR).In this approach we subtracted a simulated SRTM-X interferogram from a single-pass TSX/TDX (TerraSAR-X/TanDEM-X) interferogram.In order to compare the DInSAR derived estimate with a more common method, we constructed a DEM from the TSX/TDX acquisition and calculated the surface elevation differences between the two InSAR derived Published by Copernicus Publications on behalf of the European Geosciences Union.

N. Neckel et al.: Recent mass balance of the Purogangri Ice Cap
DEMs.Finally, the results of both approaches were converted to mass changes showing the enormous potential of the new TSX/TDX data for geodetic mass balance estimates in remote regions.

Shuttle Radar Topography Mission
The SRTM was conducted from 11 to 22 February 2000 in order to derive a near-global DEM (Rabus et al., 2003;Farr et al., 2007).For this purpose, data from C-and X-band SAR were acquired.In this study we used the X-band DEM, which was processed at the German Aerospace Center (DLR).In comparison to the C-band ScanSAR system with a swath width of 225 km, the X-band SAR system was operated with a swath width of only 45 km leaving large data gaps in the resulting X-band DEM (Rabus et al., 2003).Fortunately, 90 % of the PIC is covered by the data set (Fig. 1).The DEM is sampled to a grid posting of one arc second and is referenced to the WGS84 ellipsoid.

TerraSAR-X and its add-on for digital elevation measurements (TanDEM-X)
TSX was launched in June 2007 followed by its twin satellite TDX in June 2010.The two satellites are flying in a unique helical formation acting as a single-pass InSAR system with a flexible baseline configuration (Krieger et al., 2007).The main goal of the TanDEM-X mission is the generation of a global consistent DEM with a 12 m × 12 m grid posting and a vertical accuracy of < 2 m (Moreira et al., 2004).
In this study, we employed the experimental Co-registered Single look Slant range Complex (CoSSC) product acquired in bistatic InSAR stripmap mode on 26 January 2012 (Table 1).The data were acquired on an ascending satellite pass with an incident angle of 44 • and VV polarization.The CoSSC product has already been focused and co-registered at the TanDEM-X Processing and Archiving Facility (PAF) via the integrated TanDEM processor (Duque, 2012).For the interferometric processing of the CoSSC product we employed the GAMMA SAR and interferometric processing software (e.g., Werner et al., 2000).
Besides the tandem acquisition, we employed two TSX stripmap scenes acquired on 2 and 13 August 2011 (Table 1).The data were acquired on a descending satellite pass with an incident angle of 31 • and HH polarization.From the phases of these scenes a coherence image was calculated, which was used to support and validate the Landsat derived glacier outlines shown in Fig. 1.Further details concerning the production of the DEM and the generation of the glacier outlines are described in the Methods section.

Landsat
In order to compare the interferometrically derived results with changes in glacier extent, we utilized six Landsat Enhanced Thematic Mapper Plus (ETM+) and two Landsat Thematic Mapper (TM) scenes (Table 1).We used all bands of the orthorectified level T1 products provided by the United States Geological Survey (USGS).No horizontal shift was observed by visual comparison amongst the Landsat imagery, the co-registered TerraSAR-X coherence image and the SRTM-X DEM.In order to reduce data gaps induced by scanline errors in the 2012 ETM+ data, we combined two ETM+ scenes acquired in September 2012 (e.g., Chen et al., 2011).To enhance the spatial resolution of the Landsat ETM+ scenes to 15 m, pan-sharpening employing principal component analysis was performed.

Delineation of ice
In order to measure any glaciological parameter, the delineation of the ice body is an essential step.Glacier outlines are available globally through the Randolph Glacier Inventory (Arendt et al., 2012).However, for the PIC these glacier outlines are spatially inaccurate when compared to the Landsat scenes listed in Table 1.Due to the high spectral contrast between glacier area and non-glacier area in our study region, we conducted an unsupervised 2-class classification to delineate the ice body for the years 2000 and 2012.This classification was performed on the pan-sharpened Landsat ETM+ scenes of 2000 and 2012, which incorporate the principal components of all  Landsat bands.We also tested a method based on band ratios and the use of a specific threshold as suggested by Paul et al. (2004) and Bolch et al. (2010b).For our study region, a visual comparison showed a good match between both methods.
Choosing only two classes revealed a general determination of on-glacier and off-glacier areas.Due to a few snow covered mountain ridges and remains of snow avalanches, the derived glacier masks were adjusted manually according to the Global Land Ice Measurements from Space standards (Rau et al., 2005).In consequence of occasional cloud cover in the Landsat ETM+ data of 2012, we employed an additional coherence image of an 11 day repeat-pass of TerraSAR-X acquired on 2 and 13 August 2011 (Table 1).Due to the short wavelength of TerraSAR-X (3.1 cm) and the changing surface characteristics of the ice cap (e.g., precipitation, wind drift, melting, freezing), a low coherence was calculated for ice covered areas whereas the coherence outside the glacier was distinguishably high (e.g., Frey et al., 2012).Two small cloudy patches in the 2012 glacier mask were fixed manually on the basis of the coherence image.Ice divides were estimated on the basis of watersheds extracted from the SRTM-X DEM utilizing the r.watershed tool implemented in the GRASS GIS software (GRASS Development Team, 2012).
An error estimation of the resulting glacier masks was performed based on a modified approach from Granshaw and Fountain (2006) and Bolch et al. (2010a) who compared their semi-automatically generated results to independently digitized glacier outlines from high resolution aerial imagery at random locations.In this study, we used the additional TerraSAR-X coherence image as a reference base.Randomly selected parts of the TerraSAR-X coherence image were digitized manually and compared to the 2012 Landsat derived glacier extent.On the basis of these differences we calculated a mean relative error of ± 2.3 % for the estimated glacier area.It should be noted that this error estimate is rather conservative as it also includes one year of glacier change.As the 2000 glacier mask was derived with the same method we assume a similar error for the 2000 glacier extent.

Glacier elevation changes
In this study, glacier elevation changes were calculated using two different methods.The first method is based on differential SAR interferometry while the second method uses www.the-cryosphere.net/7/1623/2013/The Cryosphere, 7, 1623-1633, 2013 common DEM differencing.Both approaches are using the TSX/TDX acquisition from January 2012 and the SRTM-X DEM acquired in February 2000 (Table 1).The main focus of this section is on the DInSAR approach followed by a brief description of the DEM differencing.
The interferometric phase of the single-pass TSX/TDX interferogram may be described by where φ TSX/TDX is the difference of phases φ TSX and φ TDX simultaneously acquired by TSX and TDX in January 2012 (Fig. 2, upper left).As the baseline of this satellite pass is relatively small (155 m) and the data were acquired simultaneously, the same atmospheric conditions are assumed for both SAR antennas, which sets the atmospheric contribution φ atm in Eq. ( 1) to zero.The same applies for φ scat , which describes the phase difference due to different scattering on the ground.This only leaves φ orbit , which is the phase difference induced by the different acquisition geometry of the SAR sensors, and φ topo , which describes the phase difference induced by topography.
The DInSAR approach applied in this study can be described by where φ SRTM−X is the interferometric phase of the February 2000 SRTM-X acquisition.However, as the raw interferometric data of the SRTM-X acquisition is not available, φ SRTM−X was simulated from SRTM-X DEM data using the satellite geometry and baseline model of the TSX/TDX pass from January 2012 (Fig. 2, lower part).Since the same satellite geometry and baseline model was used for both interferograms, the differential phase φ diff is solely based on changes in φ topo between data acquisitions (Fig. 2, bottom).In order to solve the 2-π ambiguity of the difference interferogram, a phase unwrapping step was performed using GAMMA's minimum cost flow (MCF) algorithm (Costantini, 1998).In the next step 20 ground control points (GCPs) were selected randomly in off-glacier regions across the entire scene.The relative differential phase was set to zero at these points as no changes in topography were assumed in these regions between 2000 and 2012.
Before calculating the difference interferogram and prior to the simulation of φ SRTM−X , precise horizontal offset registration and fitting between the TSX/TDX and the SRTM-X data set is mandatory.Therefore, we calculated a DEM from the TSX/TDX interferogram.For this we removed the orbital contribution from the interferogram ( φ orbit ) by subtracting a simulated flat-earth phase trend (e.g., Rosen et al., 2000).The resulting flattened interferogram was unwrapped using GAMMA's MCF algorithm and converted to a DEM (Fig. 2, upper part).For the calculation of absolute heights, GCPs were obtained from the respective off-glacier pixel locations in the SRTM-X DEM.In the interferometric processing of the TSX/TDX DEM we observed an artificial linear  phase ramp presumably related to an inaccurate flat-earth estimate.This linear phase ramp was minimized by an additional baseline refinement based on off-glacier phase values from the differential interferogram.In order to calculate horizontal offsets between both data sets, two SAR images were simulated from the DEMs using the orbital parameters of the TSX/TDX pass (e.g., Kropáček et al., 2012).Offsets were calculated using cross correlation optimization of the simulated SAR images employing GAMMA's offset_pwrm module.These offsets were used to refine the horizontal data registration via a refined geocoding lookup table, which was used to translate the SRTM-X DEM from geographic coordinates into SAR coordinates (Fig. 2, central part) and, conversely, the final difference map from SAR coordinates to geographic coordinates (Fig. 5a).In order to compare the DInSAR derived elevation changes, we also constructed a difference map by common DEM differencing (Fig. 5b).Here we used the DEM that was created to estimate the horizontal offsets between the TSX/TDX data set and the SRTM-X DEM (Fig. 2, upper part).This DEM was geocoded to geographic coordinates with a grid posting of one arc second using the refined geocoding lookup table created above (Fig. 2, central part).This way both DEMs feature the same grid posting and are horizontally aligned.Figure 3 shows a subset of the map derived by common DEM differencing before (Fig. 3a) and after applying the refined geocoding lookup table to the data (Fig. 3b).
In the next step, a constant vertical offset and a linear trend were removed from both difference maps.The latter was estimated by a two dimensional first order polynomial fit in off-glacier regions and is probably a residual not covered by the baseline refinement mentioned above.The same linear trend was also removed from the TSX/TDX DEM.Finally, the data sets were translated to a metric cartographic coordinate system with a grid spacing of 25 m × 25 m employing bilinear interpolation.

Estimation of mass changes and error computation
In order to translate the derived surface elevation changes into volume changes, an estimate of the glacier area is needed.To respect changes in glacier extent between the dates of data acquisition, we employed the geometric union of the 2000 and 2012 glacier masks (Li et al., 2012).We then multiplied the mean surface elevation changes of the 13 separated glaciers (Fig. 1) and of the entire ice cap by the respective area estimate.Finally, we applied an average ice density of 900 ± 17 kg m −3 to convert volume changes to mass changes (Gardner et al., 2013).This was done separately for the DInSAR results and for the results of the DEM differencing.For a first accuracy assessment of the SRTM-X DEM we utilized data from the Geoscience Laser Altime-ter System (GLAS) carried on-board the Ice Cloud and Elevation Satellite (ICESat).We employed the GLA 14 data product from all ICESat campaigns provided by the National Snow and Ice Data Center (NSIDC).SRTM-X surface elevations were extracted by bilinear interpolation at each ICESat footprint location.ICESat measurements were excluded from the analysis if the difference between GLA 14 and SRTM-X elevation exceeded 150 m, which can be attributed to cloud cover during the time of data acquisition.Compared to the ICESat data, we found a mean and standard deviation of 3.93 ± 2.07 m for the SRTM-X DEM (SRTM-X DEM higher).We also compared the InSAR derived TSX/TDX DEM with the ICESat data sample available in our study region.Applying the same procedure as for the SRTM-X DEM, we found a mean and standard deviation of 3.24 ± 1.11 m against the spatially limited ICESat sample (TSX/TDX DEM higher).As off-glacier SRTM-X values were used to translate the unwrapped TSX/TDX interferogram into absolute heights, a similar vertical bias is found for both DEMs.However, it should be noted that ICE-Sat measurements are only available in a relatively flat offglacier region (Fig. 1) making the ICESat sample distribution not fully representative for our study region.
For an error estimate of the derived surface elevation changes we used the Normalized Median Absolute Deviation (NMAD) in off-glacier regions, which is considered to be less sensitive to outliers than the standard deviation (Höhle and Höhle, 2009).As a certain degree of autocorrelation can be assumed in DEMs, we analyzed semivariograms and detected a decorrelation distance of ∼ 100 m (e.g., Nuth and Kääb, 2011).To avoid the effect of autocorrelation we reduced the grid spacing of non-glacier grid cells to 200 m × 200 m.In order to calculate the NMAD for a reasonable number of grid cells close to the glacier, we only employed grid cells located in a 1 km buffer around the ice cap.We estimated a NMAD of 2.38 m for the DInSAR approach and of 3.70 m for the DEM differencing.In this study, the overall error of the derived surface elevation changes is given by where the NMAD is the random part of the error while the systematic part of the error is estimated by a weighted root mean squared error (wRMSE) of the off-glacier trend residuals (Fig. 4c).From the residuals that represent the center points of the histogram bins shown in Fig. 4a we calculated a RMSE which was weighted by the histogram frequency.In this way, we estimated a systematic error of 0.32 m for the DInSAR approach and of 0.27 m for the DEM differencing.
For the overall error of mass changes we used the root of sum of squares of the estimated errors of glacier area and surface elevation differences and an error of ± 17 kg m −3 for the density assumption.

Results
Similar to Lei et al. (2012), who estimated surface elevation changes of the PIC between 1974 and 2000, we found predominant surface lowering on outlet glaciers and glacier thickening in the interior of the ice cap (Fig. 5).However, on average the mass budget of the ice cap was close to equilibrium between 2000 and 2012.Table 2 lists the surface elevation changes next to the area changes and the annual mass budget between 2000 and 2012 for the observed part of the PIC.For the observed time period we estimated an annual change rate of −0.15 ± 0.01 km 2 a −1 for the entire glacier area, suggesting a general but slow retreat of the ice cap in the last decade.It is evident in Figs.4b and 5 that rapid thinning occurred in the lower regions of the ice cap while simultaneously a modest thickening occurred in the upper parts where the major glacier area is located (Fig. 4a).
The slightly negative mass budget of the PIC is in agreement with Yao et al. (2012), Gardner et al. (2013) and Neckel et al. (2013), who estimated balanced mass changes in the northwestern part of the TP for a similar time period, contrary to the remaining parts of the TP.Also, Lei et al. (2012) found that the PIC retreats at a much slower rate than other glaciers on the TP.
Prominent positive elevation changes were found for one glacier tongue in the eastern part of the Ice Cap (Fig. 6, World Glacier Monitoring Service id: 5Z213E0012).Between August 1999 and September 2011, the glacier terminus advanced by 515 m at an average rate of 43 m per year.

Discussion
In this study, elevation changes were calculated in two different ways.The first approach makes use of the interferometric phase of a differential interferogram, while the second approach starts from two co-registered DEMs.The results of both methods show a very similar pattern of surface elevation changes with a slightly smaller error estimate for the DInSAR approach (Table 2).The different noise level of the results is evident from the histograms shown in Fig. 5 and     in Fig. 4, which shows a higher amount of noise for the result from the DEM differencing.This may be attributed to the fact that no filter was applied to the SRTM-X DEM in the case of the DEM differencing.As shown in the flowchart (Fig. 2), both approaches are using adaptive filtering prior to phase unwrapping (Goldstein and Werner, 1998).However, in the case of the DInSAR approach this filter is applied to both data sets, i.e., the filtering takes place after the differencing (Fig. 2, lower part), while for the DEM differencing the filtering is solely applied to the TSX/TDX interferogram (Fig. 2, upper part).For this reason, more noise is present in the result of the DEM differencing.The difference between both approaches (i.e., the difference between Fig. 5a and b) shows a normal distribution with a mean and standard deviation of 0.14 ± 2.18 m, suggesting that the noise, which we attribute to the SRTM-X DEM, is random.
The systematic error of the DInSAR approach, on the other hand, is estimated to be slightly higher than that of the DEM differencing.This is probably due to residual inaccuracies of the baseline refinement not covered by the additional first order polynomial fit.
A crucial point in both methods is an accurate horizontal co-registration between the input data sets.An accurate horizontal co-registration was achieved with the help of a refined geocoding lookup table (Fig. 3).As we used the same lookup table for both methods, we can assume that the differences between the two methods are not affected by a different horizontal data registration, but rather by vertical differences.
Even though we estimated a lower random error for the DInSAR approach, we cannot conclude that the DInSAR approach is more reliable.The DInSAR approach, however, provides a more homogeneous result because the random noise of both input data sets is minimized simultaneously during the adaptive phase filtering of the differential interferogram.
In general, both methods agree within their error bars (Table 2).However, the mean elevation differences of the listed individual glaciers show rather large variations for glaciers 2, 5, 6 and 7 when comparing the two methods.We attribute these variations to the random noise, as the estimated surface elevation changes are rather low for these specific glaciers.
Penetration effects of SAR signals in snow and ice are widely discussed in the literature (e.g., Rignot et al., 2001;Gardelle et al., 2012a).The uncertainty induced by these effects are difficult to quantify when comparing InSAR derived DEMs with other data sets (e.g., Sauber et al., 2005;Berthier et al., 2006;Nuth and Kääb, 2011).The great advantage of this study is that the data sets used were acquired at the same wavelength at nearly the same time of the year.However, a certain bias introduced by the X-band penetration depth may have affected our results as the snow pack properties in 2000 and 2012 were probably not identical.Another bias can be expected from snow depth variations between 2000 and 2012 for which no measurements are available.
In 2000 two ice cores were drilled near the main ice divide of the PIC by the Byrd Polar Research Center in a joint US-Chinese project (Thompson et al., 2006).The location of the drill sites is shown in Figure 1.Personal communication with Davis (2013) and Thompson et al. (2006) suggest that the PIC has no obvious firn layer, therefore we applied an ice density of 900 ± 17 kg m −3 for the conversion from volume to mass changes for the entire ice cap.In order to test the sensitivity to a change in ice density, we also applied an ice density of 800 kg m −3 , which resulted in a 10 % decrease of the total mass balance estimates shown in Table 2.
Our results suggest an almost balanced glacier regime for the PIC between 2000 and 2012, which is in agreement with Neckel et al. (2013), who found a slightly positive trend in mass balance between 2003 and 2009 for the Zangser Kangri ice cap, located ∼300 km to the northwest.However, this behavior is contrary to the remaining parts of the TP where glacier mass balances are negative (Yao et al., 2012;Gardner et al., 2013;Neckel et al., 2013).The slightly negative mass budget of the PIC may be related to the same mechanisms as the frequently mentioned Pamir-Karakoram anomaly (e.g., Hewitt, 2005;Gardelle et al., 2012bGardelle et al., , 2013)).A possible mechanism could be a compensation of the temperature driven melt-off due to an increase of precipitation at high altitudes.This is in agreement with Li et al. (2011), who showed a significant increase in annual temperature and precipitation during the period 1961-2008 from meteorological station measurements on the TP.Also Dyurgerov and Dwyer (2000) found an increase in vertical gradient due to an annual ablation increase and a simultaneous increase in accumulation at higher altitudes on several glaciers in the Northern Hemisphere.According to Raper and Braithwaite (2009), such a behavior is typical for glaciers in wet and warm climate conditions rather than in a dry and cold continental climate.
Overall, we found negative elevation changes in glacier tongue regions except for one glacier in the eastern part of the ice cap.This glacier shows thickening at the terminus while negative values are found further up the glacier (Fig. 6).
These areas could be interpreted as reservoir and receiving areas of a surging glacier (Paterson, 1994).

Conclusions
In this study, we estimated the recent mass budget for large parts of the PIC by geodetic means.We employed SRTM-X DEM data from February 2000 and compared it with TSX/TDX data acquired in January 2012 following two different approaches.Very similar results were achieved by a DInSAR based approach and DEM differencing.We estimated modest positive elevation changes in the accumulation area of the ice cap while rapid negative elevation changes were found in glacier tongue regions.This behavior might be connected to an increase in ablation due to warmer climate conditions with a simultaneous increase in accumulation due to higher precipitation rates.Overall, the mass budget of the PIC is close to equilibrium with a mean annual mass budget of −44 ± 15 and −38 ± 23 mm w.eq. a −1 (millimeter water equivalent) estimated from the DInSAR approach and from the DEM differencing respectively.In the same time period, the ice cap retreated at a relatively slow rate of −0.15 ± 0.01 km 2 a −1 .Contrary to the remaining glacier tongues, we detected one continuously advancing glacier tongue in the eastern part of the ice cap.This study shows the enormous potential of the new TSX/TDX data to derive glacier elevation changes when combined with an older elevation data set.The data situation at the PIC is preferable as the SRTM-X DEM was acquired with the same method as the new TSX/TDX data.

Fig. 1 .
Fig. 1.Landsat ETM+ image, acquired on 7 October 2000, of the Purogangri Ice Cap (PIC) overlayed by SAR image footprints and glacier outlines for the years 2000 (yellow) and 2012 (black).Bands 4, 3, and 2 are combined as red, green, and blue, respectively.ICESat footprints are indicated as small red dots.Ice core drill sites of the Byrd Polar Research Center are shown as yellow triangles.

Fig. 3 .
Fig. 3. Elevation differences estimated from DEM differencing before (a) and after the co-registration via a refined geocoding lookup table (b).Location of the data example is shown in Fig. 5b.

Fig. 4 .
Fig. 4. Glacier hypsometry in 50 m elevation bins (a).Surface elevation changes of the PIC are shown as a function of elevation, a third order polynomial fit is shown as black solid line (b).Surface elevation changes in off-glacier regions are shown as a function of elevation, a linear trend is shown as black solid line (c).

Fig. 5 .
Fig. 5. Elevation changes of the PIC between February 2000 and January 2012.Figure (a) shows the difference map derived by DInSAR while Fig.(b) shows the difference map derived from absolute heights.The glacier outline is based on the geometric union of the 2000 and 2012 Landsat derived glacier extent.Histograms show the on-glacier elevation differences.
Fig. 5. Elevation changes of the PIC between February 2000 and January 2012.Figure (a) shows the difference map derived by DInSAR while Fig.(b) shows the difference map derived from absolute heights.The glacier outline is based on the geometric union of the 2000 and 2012 Landsat derived glacier extent.Histograms show the on-glacier elevation differences. id

Fig. 6 .
Fig. 6.Positive surface elevation changes in glacier tongue region of glacier 5Z213E0012 (World Glacier Monitoring Service id).DIn-SAR derived surface elevation changes are color-coded.In the background is the TSX/TDX DEM.Glacier terminus positions are based on Landsat imagery.Location is shown in Fig. 5a.

Table 1 .
Overview of satellite data and date of data acquisition.Perpendicular baseline, B ⊥ and pass direction are given for SAR acquisitions (Fig.1).For Landsat data path and row number is given next to the spatial resolution.