Dataset of the May 2015 Karayaylak Glacier surge in eastern Pamir

Ice, Snow and Environment Over High Asia Zone II • Versions EN1 Vol 3 (4) 2018

Zhang Zhen, Liu Shiyin, Wei Junfeng

DOI: 10.11922/csdata.2018.0009.zh
PID: 21.86101/csdata.2018.0009.zh
: 10.11922/sciencedb.572

： 2018 - 04 - 11

： 2018 - 04 - 19

： 2018 - 10 - 31

1695 12 0

Abstract & Keywords

Abstract: Karayaylak Glacier (38°35′6″N–38°44′48″N, 75°7′47″E–75°22′29″E) is located in the northern slope of Jiubie Peak of Kongur Mountains, eastern Pamir, China. Its surge in May 2015 caused property losses to local residents. Monitoring and analysis of the surged glacier has important practical significance for exploring the glacier surge mechanism and for establishing a disaster early-warning mechanism. Based on Landsat images, this study obtained the surface velocities of Karayaylak Glacier from April 13, 2015 to July 11, 2015 by using feature-tracking method. With a spatial resolution of 240 m, the obtained data were stored in 32-bit floating point GeoTiff format. A statistical analysis indicates that the results have an overall accuracy of less than ± 0.42 m d^-1. Its elevation change between 2013 and 2015 was also obtained from ASTER-derived DEM. With a spatial resolution of 240 m, the data were stored in 32-bit floating point GeoTiff format. A statistical analysis of the residual error in non-glacial areas indicates that the results have an overall accuracy of about ± 0.7 m. The dataset can serve as the basis for investigating surged glaciers in this region, or as input data for building glacier dynamic models. It provides statistical support for glacier surge research. In addition, the dataset not only supports further research in risk assessment and the early warning of glacier disasters, but also ensures the economic and social development of the China–Pakistan Economic Corridor.

Keywords: Pamir, glacier surge; glacier surface velocity; glacier elevation change; Landsat; ASTER; Karayaylak Glacier

Dataset Profile

Chinese title	2015年东帕米尔高原克拉牙依拉克冰川跃动数据集
English title	Dataset of the May 2015 Karayaylak Glacier surge in eastern Pamir
Data corresponding authors	Zhang Zhen (zhangzhen@aust.edu.cn), Liu Shiyin (shiyin.liu@ynu.edu.cn)
Data authors	Zhang Zhen, Liu Shiyin,Wei Junfeng
Time range	2013 – 2015
Geographical scope	38°35′6″–38°44′48″N ,75°7′47″–75°22′29″E
Spatial resolution	30 m, 240 m	Data volume	0.4 MB
Data format	*.tif
Data service system	<http://www.sciencedb.cn/dataSet/handle/572>
Sources of funding	National Natural Science Foundation of China (41701087,41671066); Fundamental Work Program of the Ministry of Science and Technology of China (MOST) (2013FY111400); New Talents Research Program of Yunnan University (YJRC3201702); and International Collaboration Project (131C11KYSB20160061-4)
Dataset composition	The dataset consists of two subsets titled “Velocity” and “ElevationChange”, which record glacier surface velocity and glacier elevation change, respectively.

1. Introduction

Glacier surging can be defined as intermittent and rapid cyclic movements. During the period of glacier surging, the movement speed is accelerated suddenly, which leads to rapid glacier mass transfer and new distribution of the ice body. Glacier surges may not only flood the forests, roads, villages and meadows in the lower reaches, but also cause obstruction of river courses, glacial lake outburst, debris flow and ice collapse, and seriously threaten the safety of life and property.¹ Surged glaciers have been reported in all major glacierized regions of the world, and glacier surged relatively frequent in Pamir and Karakoram regions.^2-5 In May 2015, the west branch of the Karayaylak glacier (Figure 1) on the north slope of Jiubie Mountain in eastern Pamir Plateau surged suddenly, resulting in the damage of 61 families’ houses, the disappearance of hundreds of livestock and the flooding of some meadows. This event has not only attracted extensive attention from the media and researchers, but as it happened in an area along the China-Pakistan Economic Corridor which is also an important traffic node of the Silk Road Economic Belt, it gives special significance to this study. In order to understand the control mechanism of normal and surged glaciers and accurately identify the process of glaciers surging, a dynamic model needs to be established for the active glaciers. In the process of modeling, the data of glacier velocity and glacier elevation change is indispensable. However, it is difficult to conduct field observation of the glacier because of its harsh environment. Therefore, remote sensing becomes the main means of monitoring the surged glacier.^6-8 The glacier surface velocity based on remote sensing image refers to the displacement change of the eponymous feature point of two or more periods in unit time. Both the glacier surface velocity and glacier elevation change derived from remote sensing data has been widely used in glacier research.^9-12

The Karayaylak glacier is located on the north slope of Kongur Mountains (7530 m above sea level) in the eastern Pamir Plateau, with a total length of about 20.3 km. According to the second Chinese glacier inventory,¹³ glacier area is 115.6 km², which is the largest glacier in the eastern Pamir Plateau. The glacier covers a large area of debris, whose surface debris coverage rate is about 22.3%, the glacier snow line elevation is 4220 m, and the terminal elevation is 2817 m. In our study, the average surface velocity data of of the Karayaylak glacier during surging was extracted by Landsat 8 images, and the elevation change before and after surging was extracted by ASTER stereo image. The method has been largely verified,^9-12 and the data is reliable, which can provide necessary support for the further study of the glacier.

Fig.1 Overview of Karayaylak glacier

2. Data collection and processing

2.1 Data sources and digital elevation model (DEM)

In this study, Landsat 8 satellite images (Table 1) were collected before and after glacier surging (from April 13, 2015 to July 11, 2015) to extract the glacier surface velocities. The ASTER L1A stereo image pair data obtained on June 11, 2013 and June 17, 2015 were used to extract the digital elevation model (DEM) in ENVI respectively. The plane control coordinates were derived from Landsat images, the elevation coordinates were derived from topographic map DEM, the coordinate system was set as WGS1984 UTM 43N, and the spatial resolution was set at 30 m.

Table 1 Source data used in our study

Data	Time	Optimum resolution (m)	Objective
Terra ASTER	June 11, 2013	15	Extract DEM
Terra ASTER	June 17, 2015	15	Extract DEM
Landsat OLI	April 13, 2015	15	Extract surface velocity
Landsat OLI	April 29, 2015	15	Extract surface velocity
Landsat OLI	May 8, 2015	15	Extract surface velocity
Landsat OLI	May 15, 2015	15	Extract surface velocity
Landsat OLI	July 11, 2015	15	Extract surface velocity
China’s second glacial inventory	2009	/	Glacier boundary reference
Landsat OLI	October 9, 2013	15	Control point plane reference
Topographic map	1971–1976	1:50 000 (measuring scale)	Control point elevation reference

2.2 Data processing

2.2.1 Glacier velocity data processing

The glacier surface velocity is calculated by the displacement of the characteristics of the glacier surface between different images. For optical image, its displacement can be obtained through the COSI-CORR (Co-Registration of Optically Sensed Images and Correlation), a software package developed based on IDL by California Institute of Technology and then integrated in ENVI software which supports the ENVI version 5 classic (ENVI classic interface). We used the version updated in October 2014 (http://www.tectonics.caltech.edu/slip_history/spot_coseis/download_software.html). The displacement was calculated by COSI-CORR based on the feature matching and mutual correlation of the optical remote sensing image through frequency-domain mutual correlation algorithm.¹⁴ All phase correlation methods rely on the Fourier transform theory: the relative displacement between a pair of similar images is extracted from the phase difference of their Fourier transform. i₁ and i₂ represent the two images, x and y are the coordinate values, and △_x and △_y represent the displacement of x and y, respectively, then:

\({i}_{2}\left(x,y\right)={i}_{1}\left(x-{△}_{x},y-{△}_{y}\right)\)(1)

where I₁ and I₂ represent the Fourier transform obtained through the Fourier transform theorem, which can be expressed as:

\({I}_{2}\left({\mathrm{\omega }}_{x},{\omega }_{y}\right)={I}_{1}\left({\omega }_{x},{\omega }_{y}\right){e}^{-j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)}\)(2)

where ω_x and ω_y are the frequency variables in column and row, and e is the natural constant, and j is the imaginary number. Therefore, the normalized cross-spectrum (C) of the images i₁ and i₂ can be expressed as:

\({C}_{{i}_{1}{i}_{2} }\left({\omega }_{x},{\omega }_{y}\right)=\frac{{I}_{1}\left({\omega }_{x},{\omega }_{y}\right){I}_{2}^{*}\left({\omega }_{x},{\omega }_{y}\right)}{\left|{I}_{1}\left({\omega }_{x},{\omega }_{y}\right){I}_{2}^{*}\left({\omega }_{x},{\omega }_{y}\right)\right|}={e}^{j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)} \)(3)

where * denotes the complex conjugate. The Dirac delta function (\(\delta \left(x+{△}_{x},y+{△}_{y}\right)\)) is obtained by inverse Fourier transform (\({\mathcal{F}}^{-1}\)) of the normalized cross-spectrum:

\({\mathcal{F}}^{-1}\left\{{e}^{j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)}\right\}=\delta \left(x+{△}_{x},y+{△}_{y}\right)\)(4)

The relative displacement of images can then alternatively be estimated from the coordinates of the correlation peak. In case of subpixel displacements, this peak is not a Dirac delta function anymore, but a down-sampled version of a Dirichlet kernel. Therefore, the phase correlation method of COSI-CORR can be followed by two steps: (1) the relative migration of the images reconstructed by calculating the linear phase of their cross-power spectrum; (2) the relative displacement of the images determined by the position of the strict correlation peak value. In order to calculate the relative displacement, it is necessary to apply ortho-rectification to one of the images, and then use the orthophoto as the reference image to coregister another image. The Landsat data used in this study were obtained from USGS, and all the data were ortho-rectified using the ground control point of GLS2005 by USGS. Guo et al.¹⁵ verified the correction accuracy, and the results showed that the Landsat image provided by USGS after ortho-rectification had a high correction accuracy: most of the correction precision was around half pixel, and some of the images even reached the accuracy of 1/6 ~ 1/10 pixel. Therefore, it was considered that the accuracy of the ortho-rectified Landsat images satisfied the glacier surface velocity analysis. The Landsat OLI band 8 (panchromatic band, with a spatial resolution of 15 m) was calculated by the COSI-CORR, where frequency domain algorithm was used to calculate the correlation coefficient, with the reference window set to 128 and the search window set to 32. The resulting displacement data consisted of three images that recorded east-west displacement (EW), north-south displacement (NS) and signal-to-noise ratio (SNR), respectively. We removed the results of SNR < 0.9 and cloud or shadow cover, and obtained the glacier surface velocity data.

2.2.2 Elevation change data processing

The glacier elevation change could be obtained by the difference value of ASTER DEM before and after glacier surging. Different ASTER DEM might have spatial matching error, and the error was eliminated according to the procedure in Figure 2 before the difference calculation. The elevation change error (dh) caused by theoretical space matching can be described by Formulas 5 & 6:⁹

\(\frac{dh}{\mathit{tan}\left(\alpha \right)}=a*\mathit{cos}\left(b-\phi \right)+c\) (5)

\(c=\frac{\stackrel{-}{dh}}{tan\left(\stackrel{-}{\alpha }\right)}\) (6)

where α and \(\phi \) represent the slope and aspect of referenced DEM, respectively. \(\stackrel{-}{dh}\) is the overall elevation bias between the two elevation datasets. Parameters a, b and c could be calculated by regression analysis. The migration correction in the X, Y and Z directions between DEM data can be obtained by Formulas 7 & 9:

\(X=a*sin\left(b\right)\)(7)

\(Y=a*cos\left(b\right)\)(8)

\(Z=c*tan\left(\stackrel{-}{\alpha }\right)\)(9)

The stable topography of ice-free regions was selected to iterate over such processes, and the iteration was completed when the standard deviation of dh decreased less than 2% or (X² +Y² ) decreased less than 0.25. And 5% and 95% percentiles of elevation deviation were selected to eliminate the influence of outliers.¹⁶

3. Sample description

Due to the influence of cloud/snow in the accumulation area, the data results were not ideal. Considering glacier surging occurred in the ablation area, we extracted this part of data which had a fairly higher accuracy. The data in this dataset are all stored in 32-bit floating-point GeoTIFF format, and in two folders, namely, “Velocity” and “ElevationChange”. The former folder contains five files named “yyyymmdd_yyyymmdd_velocity.tif”, which represents the glacier surface velocity between two time intervals, in md^-1. The results are presented in Figure 3. There is only one file (elevationchange.tif) in the latter folder, which is the glacier elevation change in m from 2013 to 2015. The data results are presented in Figure 4 (both Figures 3 & 4 are rendered by ArcGIS).

Fig.2 Schematic diagram of DEM correction⁹

Fig.3 Glacier velocity (The background image is hillshade created from SRTM; glacier extent is from the Second Glacier Inventory of China¹³)

Fig.4 Glacier elevation change (The background image is hillshade created from SRTM; glacier extent is from the Second Glacier Inventory of China¹³)

4. Quality control and assessment

The error of glacier velocity can be evaluated by Formulas 10 & 11):

\({\mathrm{E}}_{\mathrm{d}}=\sqrt{{\mathrm{E}}_{1}^{2}+{\mathrm{E}}_{2}^{2}}\)(10)

\({\mathrm{E}}_{\mathrm{v}}={\mathrm{E}}_{\mathrm{d}}/\mathrm{d}\)(11)

where E_d and E_v are, respectively, the errors of displacement and velocity between two time intervals. E₁ and E₂ are the displacement errors in the east-west direction and north-south direction, and d is the number of days between the two dates. In our study, the east-west and north-south displacement of the ice-free stable topography was statistically calculated (Table 2), and it was found that the average value of the east-west and north-south displacement of the ice-free stable topography was close to zero. It is assumed that the stable topography of ice-free area does not shift in theory in the east-west direction and north-south direction, the actual displacement can be taken as the displacement change residue, and the error in glacierized region is similar to the ice-free residual. Therefore, we took the standard deviation of the east-west and south-north displacements of ice-free region as E₁ and E₂ , and estimated the error of the glacier surface velocity (Table 2). The parts of SNR < 0.9 are represented by null value.

Table 2 Glacier speed error

Data pair	Mean value of east-west displacement on ice-free stable topography (m)	Standard deviation of east-west displacement on ice-free stable topography (m)	Average value of north-south displacement on ice-free stable topography (m)	Standard deviation of north-south displacement on ice-free stable topography (m)	Glacier velocity error (m d^-1)
20150413_20150429	0.4	4.2	0.6	4.5	±0.38
20150429_20150508	0.1	3.6	0.8	3.3	±0.54
20150508_20150515	1.8	4.3	1.1	3.7	±0.81
20150515_20150711	1.2	4.8	0.7	4.9	±0.12
20150413_20150711	0.3	5.3	2.1	5.7	±0.09

After co-registration, the DEM difference value of two phases tends to 0 in the ice-free area. In our study, it is assumed that the elevation change value of ice-free area is 0 in theory, and the actual result can be considered as the residual of DEM elevation difference. The elevation change error (σ) can be evaluated according to the mean elevation difference (MED) and standard deviation (STDVno-glacier) of the two periods of DEM after spatial matching and correction:¹⁷

(12)

(13)

where SE is an introduced intermediate variable, namely, standard error of the mean, and N is the number of pixel elements. In order to eliminate the spatial auto-correlation effect, it is sampled at 600 m, and N is the number of pixels after spatial decorrelation. The final results showed that the mean value of residual error (MED) on the ice-free region decreased to 0.39 m, and the error of elevation change (σ) was ±0.70 m. However, there are no data in a little area because of affection by the cloud.

5. Value and significance

The remote sensing monitoring of glacier velocity could be extracted by feature matching method based on optical image or SAR image and radar interferometry. Due to the limitation of satellite revisiting period, it is difficult to obtain high coherent images in glacier area. Many studies on feature matching in mountain glaciers have proved our method reliable, and the results obtained are consistent with field measurements.¹⁸ The most successful method for remote sensing monitoring of glacier elevation change is geodesy method used in this study. The DEM required by geodesy mainly comes from the method of optical stereoscopic image pairs and radar interferometry. Radar interferometry is limited to coherence and possible penetration of snow and ice. In our study, only part of the ablation area was considered, the visual contrast was better, and the results were reliable. Therefore, to some extent, the results of our study provide data support for the study of surged glacier which lacks of field observation.

The glacier velocity data and elevation change data are the basic data for studying glacier mass balance, glacier dynamics model and glacier disaster warning and prediction. In our study, Landsat 8 images were used to extract the glacier velocity data during the glacier surging, reflecting the motion characteristics before and after glacier surged. The elevation change data before and after surging were extracted by using ASTER stereo image, which reflected the characteristics of ice redistribution. Both of them provide basic data for the development of surged glacier dynamics model.

6. Usage notes and recommendations

All data of the 2015 surged glacier in eastern Pamir Plateau were stored in GeoTIFF format, and the spatial coordinate system was WGS1984 UTM 43N. These data can be read and operated by common GIS and remote sensing software such as ArcGIS, SuperMap, ENVI and ERDAS. The daily variation of glacier speed is represented by image pixel value in unit of m d^-1, with a spatial resolution of 240 m. The glacier elevation change is represented by image pixel value in unit of m, and the spatial resolution is 30 m. The void of this data set can be processed by spatial interpolation, so that the dataset can be applied to the glacier dynamics model to study the surged glacier.

Acknowledgments

The authors would like to thank USGS for providing the Landsat and SRTM 1 data, NASA EARTHDATA for providing the ASTER data, and the Cold and Arid Regions Science Data Center for providing the Second Glacier Inventory Dataset of China.

Meier MF & Post AS. What are glacier surges? Canadian Journal of Earth Sciences 6(1969): 807-817．