Analysis of the melting glaciers in Southeast Tibet by ALOS-PALSAR data

Global warming has resulted in the melting of glaciers in the southeastern region of Tibet. This study used the InSAR time series obtained from ALOS PALSAR data to investigate the melting glacier over Southeast Tibet during 2007 2010 with the small baseline subset (SBAS) technique. Additionally, GRACE (Gravity Recovery and Climate Experiment) gravity field model issued by the Center for Space Research (CSR) was used to extract the equivalent water height (EWH) with the deduction of glacial isostatic adjustment (GIA). The results revealed that the monitoring results by InSAR were basically consistent with the EWHs from GRACE. The glacier deformation in the study area presents a downward trend overall. From the InSAR measurements, a 900km2 area within the belt subsided -1.6 cm yr-1 during 2007 2010, and the subsidence exceeded -24 cm yr-1 in some regions. On the other hand, as an auxiliary data source, the monitoring result of GRACE is large space scale and comprehensive, the rate is about -0.01 cm yr-1 from GRACE. With the variations of glacier shape obtained from MODIS (moderate-resolution imaging spectroradiometer) data, the changing rate of surface temperature was about 0.014°C yr-1. The surface temperature change is negatively correlated with the rate of the glacial subsidence with the correlation coefficient of -0.237, which reflects that the melting glacier is influenced by the temperature rising to a certain extent. Article history: Received 28 January 2018


INTRODUCTION
The polar ice sheets and continental glaciers are melting with the global warming, especially in the Greenland and Antarctic ice sheets (GIS and AIS) (Allison et al. 2009).In the last decades, the melting rate of most of recorded glaciers and ice sheets is increasing (Gardner et al. 2013;Zemp et al. 2015).Continental glaciers play an important role in the natural environment and human activities.The ice sheet and glacier meltwater flows into depression of land and then to the ocean which contributes more directly to the sea-level rise, the seawater intrusion and the land subsidence in the coastal area (Guo et al. 2016a;Ye et al. 2017).
Routine monitoring of glacier melting mainly refers to on-the-spot investigation by field operations including field stereo photogrammetry (Sun and Chen 1980;Jing et al. 2003), differential global positioning system (GPS) technique (Shangguan et al. 2008) and ground penetrating radar (Tian et al. 2014).These techniques can measure the ice-sheet surface elevations of monitoring points directly.However, the area of temperate glacier under routine monitoring is less than 1% because of complicated topography and extreme weather (Trouve et al. 2007).
Comparing with the traditional monitoring, the remote sensing technique can monitor glacial deformation in the large scale.Remote sensing monitoring of glacier melting includes mainly ICESat (Ice, Cloud, and land Elevation Satellite) with laser altimetry (Gardner et al. 2013) and interferometric radar measurements (Hwang et al. 2019).The former is limited by the geographic distribution of ICESat at ground tracks, which is failure to monitor small glaciers without any satellite pass.Due to the all-weather and alltime dynamic monitoring characteristics of microwave, the Terr.Atmos.Ocean. Sci., Vol. 30, No. 1, 7-19, February 2019 interferometric synthetic aperture radar (InSAR) has an irreplaceable role on the extraction and assessment of vertical deformation information and melting areas of continental glaciers.Differential interferometric synthetic aperture radar technique (DInSAR) (Chang et al. 2011) and advanced multiphase interferometry techniques [e.g., permanent scatter interferometry technique (Ferretti et al. 2001), small baseline subset (SBAS) interferometry technique (Berardino et al. 2002)] can provide the high spatial resolution (up to 25 cm) and high accuracy (sub-cm to sub-millimeter accuracy) products for monitoring the deformation in temporal-spatial scale.InSAR can also identify the geological structure in the glacier deformation region (Chang et al. 2010) and detect signals of geological activity (earthquakes, glacial erosion, etc.).In recent years, the mass balance detected by GRACE can also reflect the vertical glacier surface displacements (Guo et al. 2014), but the main problem is its instability and non-exclusiveness of gravitational field inversion, especially in a small spatial scale.
The southeast region of Tibet is one of the central distribution regions of monsoonal temperate glacier in China.However, the historical record of glacier change in this region on a long-term scale is insufficient.This area is located on the strong geological structure belt, and the collision and extrusion among the plates lead to irregular deformation, these bring great challenges to the traditional means of monitoring glacier deformation (e.g., on-the-spot investigation).
Glacial area in the Tibetan Plateau and surroundings is still an issue to be precisely inventoried.According to the latest research and work, more data are published in the southern Tibetan Plateau and Himalayas using new technologies and methods (Zhu et al. 2017), and results indicate that, in the past three decades, the glacial area reduction in Southeast Tibet has the rate of -0.57% yr -1 (Zhang et al. 2015).As regards the analysis on precipitation and glacial change, results show that the rapid glacier retreat could be attributed to decreasing precipitation.However, more detailed analysis and mechanism researches still have not been discussed (Wiltshire et al. 2014).The situation emphasizes the necessity for a comprehensive study of glacial area based on satellite images and in situ glacial thickness and glacial mass balance.In this paper, we will study the glacial melting in the region by the ALOS images and compute the deformation rate in the vertical direction of the glacier.In addition, the coupling relationship between glacial melting and surface temperature change is also discussed.

STUDY AREA AND DATA
The study area is located on the southeastern of Lhunze County, Tibet of China, between 27°55'N -28°31'N and 92°29'E -93°15'E, with a total area of ~3740 km 2 .The area is bordered by Bhutan to the west and the territory of Portal to the south (Fig. 1).The Yarlung Zangbo River flows through the eastern part of the region, bringing abundant water vapor here, which breeds rich glaciers, the snow-capped mountains and forests.The study area is located on the southeast of the Karakoram Kunlun-Jiali fault zone, which is the south boundary of the Qinghai-Xizang plateau with a strong right-slip feature (Armijo et al. 1989).The southern part of the Brahmaputra fault is the tectonic boundary of the southern Qinghai-Xizang plateau, which is mainly formed by the subduction and collision of the Indian plate and the Eurasia plate (Taylor and Yin 2009).Glacier-covered peaks are mainly located in the big arc of Himalaya Mountains with the average altitude of 5000 m, where the moraine and snow remain all the year round.
The study area is characterized as the plateau temperate continental monsoon, with the average annual temperature of 5.5°C and the mean annual rainfall of about 290 mm (Lobsang et al. 2010).
The main data set used in this study consists of a stack of 23 ALOS images obtained from the Japan Aerospace Exploration Agency (JAXA) (https://www.restec.or.jp/) from January 2007 to December 2010.The twelve images are the fine dual polarization (FBD) and the rest are the fine single polarization (FBS).The mean angle of incidence and the course angle of SAR sensor are 38.769° and 83.44°, respectively.The ALOS satellite equipped with an L band radar system with the wavelength of 236 mm, and the study area is located in the staggered mountains.Due to the influences of the shadow, ghost and shrink of radar images, the vegetation covered by snow and other factors, there are a lot of noises in the acquired radar echo.The PALSAR (Phased Array type L-band Synthetic Aperture Radar) sensor is more sensitive to the vertical displacement because of the long L-band and the large incidence angle.At the same time, it is less sensitive to vegetation cover in the optical spectrum.Therefore, the sensor can detect the deformation of mountain glaciers with the high accuracy.
The auxiliary data is the spectral image of Landsat 8, and the sensor TM (Thematic Mapper) includes 9 bands.The width of image is 185 × 185 km, and the spatial resolution is 30 m including a 15-m panchromatic band.These images provide high-accuracy two-dimensional information of glacier, hence the glacier coverage area can be obtained effectively.
The CSR RL05 data used to reflect the changes in the reserves of glacier and other water storages are obtained from the Gravity Recovery and Climate Explorer (GRACE), which was jointly launched by the National Aeronautics and Space Administration (NASA) and the Deutsches Zentrum für Luft-und Raumfahrt (DLR).The GRACE is an effective technique to monitor the mass redistribution of the Earth's system, which contains twin satellites flying around the near polar orbit with the altitude of 450 km (Tapley et al. 2004;Chen et al. 2006;Velicogna and Wahr 2006).The system can provide monthly global time-varying gravitational field models.The spherical harmonics are truncated to degree and order 60 with the spatial resolution of 400 km.
In order to analyse the relationship between glacier melting rate and temperature change rate, the synthetic products of surface temperature in China obtained from the MODLT1T products are used.The calculation method is to take the monthly mean value with the temporal and spatial resolutions of 1 day and 1 km, respectively.The data are released by the International Scientific Data Mirroring of CAS Computer Network Information Center (http://www.gscloud.cn/).

InSAR Processing
The main monitoring technology for this study is In-SAR.It is difficult to obtain high-density coherent points, considering that the monitoring target is mountain glaciers.Hence, the small interference vertical baselines are relatively suitable because they have a particularly small effect due to the reduced sensitivity caused by the associated noise.This paper processed the images with the SBAS time-series approach (Berardino et al. 2002) to extract mountain glacial deformation.The main idea is to obtain the inversion result of the low-resolution deformation and the residual elevation error after the multi-look processing.Firstly, all available baseline combinations are calculated.Then the interference image pairs whose time-space baselines are less than a certain threshold selected.The optimal super master image is selected to generate a multi-look differential interference figure.These interference images are divided into several sets in terms of the baseline condition.Finally, the deformation time series of each small set can be computed by the least squares method, and the time series of the whole time period is obtained by using singular value decomposition (SVD) to calculate the multiple small baseline sets.However, if this method is used to directly solve the phase, the settled phase value will have a large jump in time and will not conform to the physical laws of deformation sequence.In order to obtain a physical deformation sequence, the average phase velocity between two acquisition time can be solved as an unknown.
To derive the temporal and spatial changes of the melting glacier in southeast Tibet, 23 ALOS PALSAR images from an ascending track during January 2007 to December 2010 are used in this study.The outline of the SAR frames is shown in Fig. 1a.All vertical baseline groups can be calculated with Next ESA SAR Toolbox software (Engdahl et al. 2012), and the results are listed in Table 1.The mountain glaciers are mixed with the alpine vegetation and exposed rock.Taking the computational efficiency and redundancy into account, the main image strategy in SBAS is the combined analysis of interference graphs used in QPS (Quasi Permanent Scatterers) technique (Perissin and Wang 2012) to generate interference graph in combination with the maximum connection coherence.In this study, the coherence of the interference graph is used as the weight to measure the quality of the connection graph to judge the pros and cons of the different interference combination graphs.The coherent value is expressed as where , p i j c is the coherence of target point p of interferometric image pair (i, j) which calculated from the mean of the neighbourhood window space [win(p)], s i represents the i-th complex SAR image (with i = 1…N).Then the interferogram between the images i and j can thus be expressed as Ii, j = s i •s j * .Np represents the number of target points in the scene, and , i j c is a coherent value obtained by averaging all the absolute values of , p i j c .The interference coefficient on the coherent combination is calculated with the algorithm of minimum spanning tree (Li et al. 2005), see Fig. 2a.The corresponding radar data acquisition time are listed in Table 1.The glaciers and bare soil areas signify the high to moderate coherence in the case of medium and small baseline.The coherent of the plateau vegetation is generally low (< 0.2) even in the case of small baseline which is shown in Fig. 3.The interference of the same season (summer and winter) can obtain high coherence (> 0.8).This situation is also reflected in the year, but the coherence in the spring is generally low (Fig. 2).To obtain the ideal combination of interference, these baselines of poor coherent interference should be deleted (Fig. 2).The interference pairs with the coherence coefficient greater than 0.2 are selected according to the above results.
In the SBAS processing, the differential interferograms of the SAR data are processed by the DORIS software (http://doris.tudelft.nl/)after selecting interference pairs, and the ALOS interferograms were processed by the repeat-pass method in the DORIS (Kampes et al. 2003).Because the ALOS satellite is well controlled in a squint angle and the Doppler centroids are generally less than 100 Hz, the zero-Doppler refocusing can be used to calculate the local offset so that a precise coregistration of the SAR data can be implemented effectively.The precise coregistration is necessary for SAR interferograms caused by along-track differences in Doppler-centroids (Prats-Iraola et al. 2012;Fattahi et al. 2017).Furthermore, the Goldstein algorithm (Goldstein and Werner 1998) was used for each interference image to reduce phase noise for the next phase unwrapping.Then, we selected the high coherence points in the unwrapped differential interferograms by utilizing the coherence information and calculated the phase of each high coherence point.Then, we used spatial high-pass filtering and time low-pass filtering to remove atmospheric phase and noise respectively.The 90-m digital elevation model (DEM) derived from the NASA Space Shuttle Radar Topography Mission (SRTM) was used as the reference topography model for the topography-related phase correction and geocoding.Then, the study performed the time-series analysis of differential interferograms by the StaMPS/MTI method (Hooper 2008).
In order to extract the deformation of the glacier, it is necessary to distinguish the glacier from kinds of ground objects as much as possible.In the study, the texture characteristics of glacier in SAR image are considered.The texture is one of the most important characteristics of image classification.It relies on the local information and is not affected by the single pixel value.The texture information of the same object is stable.Based on the second-order gray level co-occurrence matrix (GLCM) principle, the texture characteristics of SAR images can be effectively confirmed (Haralick et al. 1979).The method can distinguish different objects such as glaciers, buildings, vegetation and so on.Gray level co-occurrence matrix describes the relationship between the two pixels in a local area or the whole area, or a certain distance, in fact, the texture feature information for identifying different objects is contained in the correlation between domain pixels.It is used to calculate the probability of a pair of pixels (i, j) in which the gray value is i and j with the distance of d pixels in the θ direction.The θ signifies the location and orientation of the space, mainly including horizontal, vertical and 45-degree angles.The grayscale symbiotic matrix generates 14 features around it

/ /
where p (i, j, d, θ) is the number of occurrences of i and j (the joint conditional probability density between gray levels) for the pixel pairs whose distance in the θ direction is d.It characterizes the value of element in a grayscale matrix.
Every ALOS image is grayscale-quantized, and then the size, distance and direction of the window are determined.Firstly, with speckles of a SAR image retrained by the enhanced Frost filter method, then we started with the window size of 9 × 9 in the upper left corner and calculated the texture features with the sliding step length of 1.Each texture feature was taken as a single band for color band synthesis.Taking into account the small spatial scale of the image range, the decision tree classification algorithm was used for the fast and robust classification calculation.The classification results obtained from the RF decision tree algorithm (Fang et al. 2013) are shown in Fig. 4a, and the area of glacier is about 943.4 km 2 .
In order to verify the extraction result, two images of Landsat 7 in the same period are processed by striping, radiometric calibration, registration and geometric correction.The optical images are used to supervise the classification and the indexing of glacier area by NDSI.The contour and classification results are shown in Fig. 4b, and the area is about 901.2 km 2 with the relative precision of 95.5%.By comparing the classification results (Fig. 4), we found that the scope of the glacier extracted from the texture of the SAR image is somewhat overestimated due to the fuzzy edge information.The advantages of obtaining the glacier information near the equator in the cloudy situation cannot be ignored.
The results of smoothing and filtering of two images are combined, and the confusion matrix generated from the classification is compared for the quantitative evaluation of the final result (Cohen 1960;Congalton 1991).The results are shown in Table 2.And the glacier footprint and the average LOS velocity on the actual map are shown in Fig. 5.
Finally, the outline of glacier coverage is vectorized to clip the entire image file, which contains the results of SBAS time analysis.Then, we calculate the mean of the grid values in the glacier coverage, draw the time series and calculate the rate of annual ablation over the entire glacier area by the least squares algorithm.The entire processing is shown in Fig. 6.

GRACE Processing
In order to prove the rationality of the InSAR-derived results, the CSR GRACE gravitational sphere harmonic model data (RL05) from 2007 -2010 are downloaded to extract the equivalent water height with the GIA (Glacial isostatic adjustment) correction (Jacob et al. 2012).Firstly, the model is truncated to degree and order 60 of Stokes coefficients, and the coefficient of C 20 is replaced by the observed value of SLR (Cheng and Ries 2014).Taking into account the high-degree spherical harmonic coefficient of high correlation system error, these items should be de-correlated (Swenson and Wahr 2006).The smooth filtering is used to eliminate the noise effects of the high-order spherical harmonic coefficients.Secondly, GIA is converted to the same spherical harmonic coefficient as the GRACE time-varying gravitational field model and deducted.The ICE-5G model (Peltier 2004) where a is the mean radius of the earth, ave t and w t are the mean earth's density (5517 kg m -3 ) and the water's density (1000 kg m -3 ), P lm (cosθ) is the fully normalized Legendre function, N max is the maximum order of the monthly gravity field model, k l is the l-order load Lough, θ and φ refer to the latitude and longitude respectively, and C lm d and S lm d are the changes of GRACE time-varying gravitational field model at two different times.The Gaussian filter and fan-shaped filter are used to remove the high-degree errors.The weight factor w l should be added into the Gaussian filter when the order is 1.Similarly, the sector filter needs to add the weight factors w l and w m of Gaussian filter with the order l and the number of times m, and the filtering radius is 340 km.
Whereas the GRACE time-varying gravitational field model contains some high-frequency errors, the researchers proposed filtering techniques, such as Gaussian filtering, to weaken the effects of high-frequency errors in its spherical harmonic coefficients (Wahr et al. 1998;Zhang et al. 2009).The larger the filter radius of the Gaussian filter, the higher the signal-to-noise ratio of the filtered result.But it also leads to the signal attenuation and signal leakage between different regions (Chen et al. 2013;Guo et al. 2016b).In this study, the main loam, snow, and vegetation water provided by GLDAS Noah simulate the change of water reserves in the karst area, and the grid scale factor in the study area is calculated to recover the filtering and the signal leakage caused by the spherical harmonic coefficient truncation.
To assess the accuracy of the results, we download the GRACE-Tellus (GRCTellus) gridded data (https://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons/) to compare with the results of our results.

MODIS Processing
In this study, the monthly temperature products of MODIS were used to extract and calculate the change rate of temperature.We acquired raster data of the Chinese surface temperature between 2007 and 2010 and extracted them by mask of a SAR image.The resolution of temperature data is consistent with that of SAR data.The sampled raster cell values are derived and the corresponding annual change rate is computed.After geocoding and displaying, we can investigate the relationship between glacier deformation and temperature change.

Glacier Melting Velocity Analysis
Considering that glacier movement may be very obvious in the time range, we assumed that the deformation in this area is consist of seasonal nonlinear deformations and prolonged linear subsidence trends (Kampes et al. 2004).After a computational study, we observe that the deformation results of the InSAR shown in Fig. 7 do not show a good periodicity over the past four years.Our iterative solution to the linear-triangular compound function model fails to converge.For this reason, linear fitting was performed using the principle of least squares.
The figure indicates that the elevation of glacier in the study area generally decreased with the average rate of 1.6 cm yr -1 .In 2008, the fluctuation was much larger, which was caused by the heavy snowfall and freezing disaster (Zeng 2009).Overall, the melting rate of mountain glacier is rapid.
Figure 8 shows the equivalent water high (EWH) change retrieved from GRACE RL05 time-varying gravitational field model with the GIA correction during 2007 -2010.The results serve as an auxiliary reference for monitoring results of SAR glacier deformation inversion.Considering the small spatial scale of the interesting region, the obtained result is much weaker than the sensitivity of the SAR results.But the results basically reflect the consistency in a linear trend.The results of GRACE show that the total water storage in the area was in a downward trend during the period but with an obvious seasonal cycle.

Surface Temperature Change
The temperature change rate obtained from MODIS products is shown in Fig. 9. Generally, the average monthly temperature increased during the period with the rate of 0.014°C yr -1 , the maximum rate exceeded 0.07°C yr -1 , which is consistent with the published 0.014 -0.08°C yr -1 of the relevant study (Du et al. 2016).

Correlation Analysis
Firstly, the GRACE-based EWH shown in Fig. 8 is computed as the mean of the study area.Then the correlation analysis between surface temperature and the EWH shows that the variation of EWH is consistent with that of surface temperature in the same month.Figure 10a shows this corresponding correlation scatter plot, the average rate of change of surface temperature in the study area is the Xaxis, and the result of EWH is the Y-axis.The correlation coefficient is 0.673 (Table 3).In order to calculate the correlation between the temperature change and the melting velocity of glacier, t, the raster temperature data in the glacier coverage area is selected to obtain the average temperature change rate, this result is shown in Fig. 10b which the result of InSAR is the Y-axis.It can be seen that the correlation coefficient is -0.239 (Table 4), which manifests that the rising temperature could accelerate the melt of glacier.From the analysis of the GRACE results, it is true that the increase in temperature does contribute to the extent of glacial ablation in the region, but the results of interferometric SAR also indicate that the increase in temperature is not the most important cause of glacial ablation in the study area.

CONCLUSION
This study used the ALOS-PALSAR images to investigate the deformation of glacier in southeastern Tibet from 2007 -2010.Moreover, multiple data sets were used to study the potential factors of deformation.After analysing the coherence characteristics of glaciers and other ground objects, it was found that the study area was so complex that it was difficult to distinguish the glacier.In the inverse process of the glacier ground objects, the QPS method was introduced to solve the difficult problem to extract the interesting information in the mountainous area.The scene object classification and the glacier coverage area were added, and the SAR image texture feature and spectral band calculation method were used to solve and optimize the classification and extraction of the ground objects.
The results indicated that the glaciers in the southeast Tibet were suffering from strong melting, and the glaciers in upper and lower two deformed regions divided by a fault zone decreased with different rates in a quasi linear manner in the four years.The maximum subsidence rate reached about -2.4 cm yr -1 between 2007 and 2010.Some subsidence points around the glacier were measured and found that the movements of most of them were non-linear.In addition, InSAR observations contribute to the model construction of glacier displacement and explain the subsidence of ice in different regions.Compared with the change of equivalent water high, including GIA and other factors, GRACE plays an important role in extracting the change of water reserves in the large scale.SAR is more sensitive than GRACE in extracting the vertical deformation of glacier in the small scale.
At the same time, the results of MODIS monthly temperature show that there was a correlation between the rates of temperature and glacier deformation in a seasonal trend.It is noted that the deformation time series of GRACE had a positive correlation with the surface temperature change rate.The phenomenon shows that the land water reserves also show an upward trend when the summer temperature rises.The reason is that the results of GRACE monitoring represent the change of water storage in the whole region, not only related to glacier melting water caused by land surface temperature change, but also to local precipitation.Therefore, the monitoring results of InSAR can be used as a reference for the separation of glaciers and other water elements by GRACE satellite.Because of the weak correlation between the rates of temperature change and glacier ablation in the study area, it is clear that the temperature change is not the most critical factor for the glacier deformation in the southeastern region of Tibet, which may be related to other factors such as local precipitation and geological activity.The problems and suggestions in this study are as follows: (1) Tellus (http://grace.jpl.nasa.gov/data/get-Data/jpl_glob-al_mascons/), the monthly mass grids provided here contain global water storage anomalies relative to a time-mean as derived from GRACE time-variable gravity data.The results obtained from our inversion have a high degree of consistency with it.However, as an auxiliary reference for the InSAR derived results, the signal leakage from the Indian aspect is still not fully suppressed, but the overall trend of deformation tends to a long-term linear settlement, reflecting the rationality of the InSAR linear trend results.(3) Due to the lack of reliable high resolution temperature data from the ground point of the study area, the variation of vertical displacement of the ice glaciers obtained by InSAR inversion can not be well correlated with it.Therefore, the correlation analysis is carried out on the average of the total calculation results of the study area, and the corresponding data are given.In further work, small-scale glacial deformation researches will contain the analysis of other climatic factors such as precipitation.

Fig. 1 .
Fig. 1.Geographical location and terrain cartography of the study area.(a) Landsat 8 image of the study area.The inset at the top left signifies the location of Tibet in China.The inset at the top middle signifies the study area in Tibet.The red rectangle shows the frames of radar sensor (ALOS) used in this study.The white area is the glacier.(b) 90m-DEM (provided by SRTM).The black curves are main faults.

Fig. 2 .
Fig. 2. Multiple combines of different baselines based on minimum spanning tree algorithm (MST).(a) The different baselines combination referring to coherency.(b) The filtered interferometric pairs.

Fig. 4 .Fig. 5 .
Fig. 4. The two-dimensional information of the glacier.(a) The extraction based on texture (SAR).(b) The extraction based on NDSI index.

Fig. 8 .
Fig. 8.The EWH change from GRACE.The linear trend in red (this study) is about -0.011 cm yr -1 .The linear trend in blue (GRCTulles) is about -0.006 cm yr -1 .The error bars represent the standard deviation of the two results.

Table 3 .
The accuracy verification of InSAR results is limited by the lack of existing data.ALOS data transmission cycle is up to 46 days, data density is low and time is undersampling.Further work is to add data and provide details, such as ASAR.(2)Comparing with the results published on the GRACE Correlation analysis of EWH and STC (surface temperature change).

Table 4 .
Correlation analysis of Glacier velocity and STC (surface temperature change).