Three-Dimensional Time Series Movement of the Cuolangma Glaciers, Southern Tibet with Sentinel-1 Imagery

Abstract

Many debris-covered glaciers are broadly distributed across High Mountain Asia and have made a number of contributions to water circulation for Qinghai-Tibet Plateau (QTP). The formation of large supraglacial lakes poses risks for glacier lake outburst floods (GLOFs). Therefore, it is important to monitor the movement of glaciers and to analyze their spatiotemporal characteristics. In this study we take Cuolangma glaciers in the central Himalayas as study targets, where glacier No.1 is a lake-terminating debris-covered glacier and glacier No.2 is a land-terminating debris-covered glacier. The 3D deformation time series is firstly estimated by using the Pixel Offset-Small Baseline Subsets (PO-SBAS) based on the ascending and descending Sentinel-1 datasets spanning from January to December 2018. Then the horizontal and vertical time series displacements are obtained to show their spatiotemporal features. The velocities of glacier No.1 in horizontal and vertical direction were up to 16.0$$\pm$$0.04 m/year and 3.4$$\pm$$0.42 m/year, respectively, and the ones of the glacier No.2 were 12.0$$\pm$$0.07 m/year and 2.0$$\pm$$0.27 m/year, respectively. Next, the correlation between the precipitation and the surface velocity suggests that the glacier velocity does not show a clear association with daily precipitation alone. Finally, the debris-covered glaciers evolution is evaluated which shows that the tongue of the glacier No.1 is wasting away and the transition of glacier No.2 from land-terminating to lake-terminating is a probable scenario in the later period of glacier wastage. This research can significantly serve for glacier multidimensional monitoring and the mitigation of hazardous disaster caused by debris-covered glaciers in the central Himalayas.

1. Introduction

Mountain glaciers are one of the most sensitive natural indicators for climate changes in the Himalayas [1]. The melting water of mountain glaciers is one main live resource for downstream people, especially for the arid and inland areas. Furthermore, debris-covered glaciers play a significant role in water-cycles in High Mountain Asia (HMA), altering glacier melting and their spatial patterns [2]. Shrinking of Himalayan debris-covered glaciers poses challenges to societies, such as the increasing changes of seasonal runoff caused by the glacier melting and the increasing risk of glacier lake outburst floods (GLOFs) caused by the expansion of unstable glacial lakes [3]. The percentage of debris-covered glaciers are roughly 10% in the Himalayas [2] and their tongues are generally situated in the lowest elevation with the increasing number of glacial lakes at the terminate [4]. Additionally, previous studies show that mass loss of debris-covered tongues were at the same rates as debris-free glaciers [5,6,7], but the lake-terminating glaciers mass loss were more negative than land-terminating glaciers over the period of 2001–2015. [3]. The principal reason may be the significant increased rate of supraglacial lakes melting [2].

The estimation of debris-covered glacier flow velocity provides better insight into glacier spatiotemporal evolution. The variation of glacier flow velocity and the mass balance are largely determined by their susceptibility to climate change. Recent studies show that flow velocities are 100–200 m/year near Bhutan on the north side of Himalayas whereas the maximum flow velocities on the south side is around tens of meters, emphasizing a great spatial variation of glacier velocity in the Himalayas [8,9,10]. However, the flow velocities for debris-covered tongues varied from 0 to 37 m/year on the south side of Everest region [10]. Additionally, the dynamic behaviors of lake-terminating glaciers are in contrast to the land-terminating glaciers in the Himalayas [4]. Since 2000, the Himalayas have experienced considerable reduction of flow velocities due to thinning and driving stress reduction [3]. Therefore, it is crucial to conduct multidimensional movement monitoring of debris-covered glacier.

Most of debris-covered glaciers are difficult to approach and study with regular in-situ measurements. Satellite remote sensing techniques including optical feature-tracking and interferometric synthetic aperture radar (InSAR) have been abroad used for mapping comprehensive, nonuniform, and frequent glacier velocities [11,12,13,14,15]. However, optical remote sensing is often hampered by the inclement weather conditions [16]. Phase-based InSAR and amplitude-based pixel-tracking methods have become significant instrument to study glacier movement in a tradition, low-priced manner [17]. Phase-based InSAR methods including differential InSAR(D-InSAR) and multiple aperture interferometry (MAI) are always affected by the loss of coherence caused by the large time interval of SAR image pairs and large gradient glacier movement [18,19]. In contrast, the amplitude-based pixel-tracking method can work very well to estimate glaciers displacements with rates ranging from decimeter per year to tens meter per year in both line-of-sight (LOS) and azimuth directions [20]. With accessible to multi-track SAR data, it is possible to acquire ascending and descending SAR data over the same region during the specific period, which can then be used to construct 3D deformation models of glaciers [21,22]. More importantly, the time series of 3D motion can also be derived with the pixel-offset small-baseline subset (PO-SBAS) technique [23]. In recent years, numerous geological disasters have occurred in the northern of Himalayas, Tibet, China [24,25]. Generally, it is fundamental to comprehend the full 3D flow pattern for interpreting of debris-covered glacier movement and particularly recognizing their evolution process [26]. However, what was the 3D debris-covered glacier deformation pattern and how did it change with time? And what was the relationship between the glacier movement and precipitation? How does a moraine-lake form in a debris-covered glacier? We take the Cuolangma glaciers as example to answer these questions.

In this study, we derive 3D time series displacement of two different terminus types of Cuolangma glaciers in the central Himalayas by applying PO-SBAS methods with Sentinel-1 imagery. We examine the 3D time series movement of Cuolangma glaciers in detail, where the one is a lake-terminating glacier and the other is a land-terminating glacier. Then the deformation features in terms of precipitation and seasonality are analyzed. A general explanation for evolution of debris-covered glaciers will trace the development of a glacial lake, verifying the role of lake-terminating debris-covered glacier in catastrophic disaster.

2. Study Area and Datasets

2.1. Study Area

Cuolangma glaciers are located at the Mid-Himalayas, which consists of a series of parallel ranges from the Ganges Basin to the QTP [16,27]. It should be stressed that the status of debris-covered glacier and its surroundings in the Himalayas have been paid much attention as the increased risks of GLOFs [16]. Since the 1980s reduced height of ice dam and accelerated accumulation of glacial lake water occurred with the rising temperatures [4]. In addition, the Mid-Himalayas provide waters for major rivers in Asia, including the Ganges, Indus, and Yangtze Rivers [27]. The prevailing westerlies and the Indian summer monsoon have a profound contribution to complicated climate pattern [28], as shown in Figure 1a. In summer, the Indian summer monsoon brings significant rainwater to the south side of Himalayas, resulting in annual precipitation of about 2500–3000 mm in the eastern Nyainqentanglha mountains [29]. In winter, the westerlies have played a determined role in precipitation in the west of Himalayas [6,28]. The two glaciers have debris covering in their ablation areas, which provides facilities for supraglacial ponds’ growth [16]. We take the two different type glaciers of Cuolangma glaciers as the study object. The location and topographic features can be seen in Figure 1, where Figure 1b and Figure 1c show the terrain of glacier No.1 and No.2, respectively. The area of glacier No.1 is 24 ${\mathrm{km}}^2$ and the elevation varies from 5500 to 6100 m, while the area of glacier No.2 is 26 ${\mathrm{km}}^2$ within the elevation from 4900 to 5800 m. The glacier No.1 is a lake-terminating glacier; the supraglacial ponds and lakes are typically studded on the tongues of the large debris-covered glacier [2]. However, the glacier No.2 is a land-terminating debris-covered glacier with the likely consequence of transiting into a lake-terminating glacier.

2.2. Datasets

To generate 3D movement maps of the Cuolangma glaciers, we processed Sentinel-1 images acquired from both ascending and descending tracks with Interferometric Wide Swath (IWS) mode. The pixel spacing in the LOS and azimuth directions are 2.33 m and 13.96 m, respectively. Table 1 presents the principal parameters for SAR imagery, whilst Figure 1a shows their coverages superimposed on the SRTM DEM shaded relief image. In total, 30 ascending and 30 descending SAR images were acquired from tracks 85 and 121, respectively, spanning from January 2018 to December 2018. From these, a total of 173 high-quality offset pairs which meet the co-registration accuracy with 0.02 pixels, were selected as shown in Figure 2 including 85 ascending pairs and 88 descending pairs. Additionally, the daily rainfall products from the Global Precipitation Measurement (GPM, https://gpm.nasa.gov) were involved for triggering factor analysis.

3. Methods

Figure 3 illustrates the flowchart used to retrieve the time series of 3D glacier movement, including three steps: offset-tracking method, estimation of 3D displacement, and the time series of 3D displacement.

3.1. Offset-Tracking

Offset-tracking, including amplitude tracking and coherence tracking, is based on the calculation of the offsets in both range and azimuth directions by maximizing the cross-correlation of the SAR image patches, which have been used to measure surface displacement caused by seismic activity, coal exploration, and glacier movements. [18,20,30,31,32,33,34,35,36]. Amplitude tracking correlates two images using the amplitude information and is suitable for situations with clear ground features or poor coherence of the image pair. Coherence tracking correlates two images by calculating coherence within the given patch if no clear features are available [20]. The amplitude-based method is used in this study because it is more sensitive to large gradient displacement measurement as seen in the research area and it has a high computation efficiency [20,21,22,23,24].

Normally, the calculated total offsets contain four elements, that is, the offsets caused by (1) glacier movement, (2) differences in satellite orbit position, (3) ionospheric delays, and (4) the topographic relief [20]. The offsets caused by ionospheric delays are generally ignored for C-band Sentine-1 SAR images in mid and low-latitude areas [20,23,37]. Previous work [38,39] confirmed that topography errors can be calculated by this mathematical model, as shown in Equations (1) and (2):

(1)${\mathrm{Topo}}_{\mathrm{az}}=\sin \mathsf{\alpha }\cot \mathsf{\theta }\frac{\Delta \mathrm{h}}{{\mathrm{R}}^{\mathrm{az}}}$

(2)${\mathrm{Topo}}_{\mathrm{rg}}=-\frac{\mathrm{B}}{\mathsf{\rho }\sin \mathsf{\theta }}\frac{\Delta \mathrm{h}}{{\mathrm{R}}^{\mathrm{rg}}}$

where $\mathsf{\alpha }$ is the azimuth angle, $\mathsf{\theta }$ is the look angle of the SAR images, $\mathsf{\rho }$ is the slant range, ${\mathrm{R}}^{\mathrm{az}}$ is the azimuth pixel spacing, $\mathrm{B}$ is the perpendicular baseline,${\mathrm{R}}^{\mathrm{rg}}$ is the slant-range pixel spacing, and $\Delta \mathrm{h}$ is the topographic error. In this study, we estimate and remove the offsets caused by height difference based on Equations (1) and (2).

Three main steps are employed to estimate LOS and azimuth offsets by GAMMA [23,40]: (1) precise registration of SAR images, (2) the calculation of offsets, (3) the separation of the deformation signals from offset measurements. As the offset pairs with short temporal intervals ($\le 48$ days) and spatial baselines ($\le 200$ m) were used to obtain initial offset fields, we first chose the window size of 256$\times$256 pixels in azimuth and range directions. Then, the window sizes of 128$\times$128 pixels and 64$\times$64 pixels were used to refine the results adaptively [37]. We set the correlation coefficient (CRC) threshold at 0.2, and any patches below this value were masked. Lower CRC values usually appear for the offset pairs with large spatial baseline, large-gradient movement, and temporal changes in the surface scattering characteristics [23,37].

3.2. Estimation of 3D Deformation Fields

As the difference of acquisition dates between adjacent ascending and descending Sentinel-1 SAR data is only 3 days, it is reasonable to derive 3D time series movement by fusing ascending and descending SAR measurements. The range and azimuth offsets derived from SAR offset-tracking are the function of the real 3D surface deformations [23]. Taking the ascending images as an example, the ground-range (GR) deformation ${\mathrm{D}}_{\mathrm{H}}^{\mathrm{A}}$ can be expressed as the function of north-south deformation ${\mathrm{D}}_{\mathrm{N}}$ and east-west deformation ${\mathrm{D}}_{\mathrm{E}}$, whilst the LOS deformation ${\mathrm{D}}_{\mathrm{LOS}}^{\mathrm{A}}$ is the function of vertical deformation ${\mathrm{D}}_{\mathrm{U}}$ and ground-range deformation ${\mathrm{D}}_{\mathrm{H}}^{\mathrm{A}}$, as shown in Equations (3) and (4) [23]:

(3)${\mathrm{D}}_{\mathrm{H}}^{\mathrm{A}}={\mathrm{D}}_{\mathrm{E}}\sin \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)+{\mathrm{D}}_{\mathrm{N}}\cos \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)$

(4)${\mathrm{D}}_{\mathrm{LOS}}^{\mathrm{A}}=-{\mathrm{D}}_{\mathrm{U}}\cos {\mathsf{\theta }}^{\mathrm{A}}+{\mathrm{D}}_{\mathrm{H}}^{\mathrm{A}}\sin {\mathsf{\theta }}^{\mathrm{A}}$

where ${\mathsf{\alpha }}^{\mathrm{A}}$ is azimuth angle of the SAR images, and ${\mathsf{\theta }}^{\mathrm{A}}$ is the look angle. Based on Equations (3) and (4), the relationship between LOS deformation and real 3D surface deformation can be deduced as follows [23,35]:

(5)${\mathrm{D}}_{\mathrm{LOS}}^{\mathrm{A}}=-{\mathrm{D}}_{\mathrm{U}}\cos {\mathsf{\theta }}^{\mathrm{A}}+ {\mathrm{D}}_{\mathrm{E}}\sin \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)\sin {\mathsf{\theta }}^{\mathrm{A}}+{\mathrm{D}}_{\mathrm{N}}\cos \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)\sin {\mathsf{\theta }}^{\mathrm{A}}$

Meanwhile, the azimuth deformation ${\mathrm{D}}_{\mathrm{AZ}}^{\mathrm{A}}$ can be expressed as the function of east-west ${\mathrm{D}}_{\mathrm{E}}$ and north-south ${\mathrm{D}}_{\mathrm{N}}$ deformations [23]:

(6)${\mathrm{D}}_{\mathrm{AZ}}^{\mathrm{A}}=-{\mathrm{D}}_{\mathrm{E}}\cos \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)+{\mathrm{D}}_{\mathrm{N}}\sin \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)$

Therefore, the real 3D deformation of the glacier can be retrieved using the azimuth and LOS offsets derived from ascending and descending images. For simplicity, the matrix form is shown below.

(7)$\mathrm{BX}=\mathrm{L}$

where

$$\begin{gathered} \mathrm{B}=[\begin{array}{ccc}-\cos {\mathsf{\theta }}^{\mathrm{A}}& \sin \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)·\sin {\mathsf{\theta }}^{\mathrm{A}}& -\cos \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)·\sin {\mathsf{\theta }}^{\mathrm{A}}\\ 0& -\cos \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)& \sin \left({\mathsf{\alpha }}^{\mathrm{A}}-3\mathsf{\pi }/2\right)\\ -\cos {\mathsf{\theta }}^{\mathrm{D}}& \sin \left({\mathsf{\alpha }}^{\mathrm{D}}-3\mathsf{\pi }/2\right)·\sin {\mathsf{\theta }}^{\mathrm{D}}& \cos \left({\mathsf{\alpha }}^{\mathrm{D}}-3\mathsf{\pi }/2\right)·\sin {\mathsf{\theta }}^{\mathrm{D}}\\ 0& -\cos \left({\mathsf{\alpha }}^{\mathrm{D}}-3\mathsf{\pi }/2\right)& \sin \left({\mathsf{\alpha }}^{\mathrm{D}}-3\mathsf{\pi }/2\right)\end{array}] \\ \mathrm{X}={\left[\begin{array}{ccc}{\mathrm{D}}^{\mathrm{U}}& {\mathrm{D}}^{\mathrm{E}}& {\mathrm{D}}^{\mathrm{N}}\end{array}\right]}^{\mathrm{T}} \\ \mathrm{L}={\left[\begin{array}{cccc}{\mathrm{D}}_{\mathrm{LOS}}^{\mathrm{A}}& {\mathrm{D}}_{\mathrm{AZ}}^{\mathrm{A}}& {\mathrm{D}}_{\mathrm{LOS}}^{\mathrm{D}}& {\mathrm{D}}_{\mathrm{AZ}}^{\mathrm{D}}\end{array}\right]}^{\mathrm{T}} \end{gathered}$$

We use the least squares criterion (${\mathrm{V}}^{\mathrm{T}}\mathrm{PV}=\min$) to obtain the optimal estimation of 3D displacements and error as follows,

(8)$\mathrm{X}={\left({\mathrm{B}}^{\mathrm{T}}\mathrm{PB}\right)}^{+}{\mathrm{B}}^{\mathrm{T}}\mathrm{PL}$

3.3. Time Series of 3D Displacement

When we establish the adjustment model for the 3D (east-west, north-south, up-down) displacement, the 3D time series displacement can be estimated using the singular value decomposition (SVD). Taking the north-south direction as an example, we first generate M offset pairs with n+1 SAR images:

(9)$\frac{\mathrm{n}+1}{2}\le \mathrm{M}\le \mathrm{n}\left(\frac{\mathrm{n}+1}{2}\right)$

Second the north-south displacement between each offset pair is independently calculated using the offset-tracking method and the estimation of 3D displacement.

(10)${\mathsf{\delta }\mathrm{N}}_{\mathrm{j}}={\left[\begin{array}{ccc}{\mathsf{\delta }\mathrm{N}}_1& \cdots & {\mathsf{\delta }\mathrm{N}}_{\mathrm{M}}\end{array}\right]}^{\mathrm{T}}$

where ${\mathsf{\delta }\mathrm{N}}_{\mathrm{j}}\left(\mathrm{j}=1,\cdots ,\mathrm{M}\right)$ are the north-south displacement between each offset pair.

Next, in the M offset pairs, the acquisition orders of primary and secondary images are ${\mathrm{m}}_{\mathrm{j}}$ and ${\mathrm{s}}_{\mathrm{j}}$($\mathrm{i}=1,\cdots ,\mathrm{M}$), respectively. If the primary and secondary images are ranged in order, thus,

(11)${\mathsf{\delta }\mathrm{N}}_{\mathrm{j}}=\mathrm{dN}\left({\mathrm{t}}_{{\mathrm{m}}_{\mathrm{j}}}\right)-\mathrm{dN}\left({\mathrm{t}}_{{\mathrm{s}}_{\mathrm{j}}}\right)\left(\mathrm{j}=1,\cdots ,\mathrm{M}\right)$

Equation (11) includes M equations with N unknown parameters and can be rewritten in matrix form as:

(12)$\mathsf{\delta }\mathrm{N}=\mathrm{A}\left(\mathrm{dN}\right)$

where A = $\left[\begin{array}{ccccc}0& -1& 0& +1& \cdots \\ 0& 0& +1& 0& \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \end{array}\right]$, each row of matrix A [$\mathrm{M}\times \mathrm{N}$] corresponds to an offset pair.

Finally, SVD is used to solve Equation (12) to derive the time series deformation of the glacier in north-south direction. Similarly, the time series of other two directions can be achieved in a similar way.

4. Results and Analyses

4.1. 3D Glacier Velocities Fields

In this study, the 3D glacier velocities were derived by PO-SBAS method to better explain their complicated movement. Figure 4 shows the horizontal and vertical velocity of these two glaciers. The map indicates that the most active part of glacier No.1 in horizontal velocity is at the downstream part. While the most active part of the glacier No.2 is at the upper part. Identically, as for the vertical component, the downward velocity happens in the western of the two profiles (AA′, BB′), while the upward velocity occurs in the eastern part.

Figure 5 illustrates the relationship between the elevation and the glacier surface velocity, vertical velocity. The four profiles are shown in Figure 4, in which AA′ and CC′ are the flow direction of glacier No.1 and glacier No.2 respectively; BB′ and DD′ intersect the flow profiles of glacier No.1 and glacier No.2 at their widest places. The error bar represents the standard deviation (STD) of the 3D velocities in the range parallel to the selected profiles.

The horizontal velocity is increasing gradually from the upper part to the middle part, with a maximum annual velocity of 15.77$\pm$0.04 m/year. Immediately, the glacier flows towards the lower section. At a distance of 5 km from the initial point along profile AA′, the horizontal velocity is up to 14.46$\pm$0.12 m/year, which is the second maximum velocity, behind the middle part due to the narrow proglacial tongue in Figure 5a. The vertical velocity of mainstream direction is around zero as a whole except a significant abrupt change within 0.5 km. From the perspective of the STD of glacier velocity, the errors in horizontal and vertical directions are mainly distributed in two locations. Within 2 to 3 km, there is an obviously increased velocity in horizontal direction. Immediately afterwards, in the tongue part of 4.0 to 6.0 km, the horizontal velocity gradually slows down. At the same time, the velocity in vertical direction suddenly makes change. The upward velocity of 1.88 ± 0.23 m/year occurred and then the glacier experienced downward velocity of −0.25 ± 0.14 m/year. The unstable flow velocity at these locations of glacier No.1 may be due to the complicated mechanism between the tongue and the glacial lake.

The horizontal velocity along profile CC′ in Figure 5c firstly increased and gradually decreased, where the maximum velocity is 11.0$\pm$0.07 m/year at 2 km. The velocity in vertical direction for glacier No.2 is 2.0$\pm$0.27 m/year. The unstable flow velocity of profile CC′ mainly occurs within 1 km. The maximum velocity errors are 3.45 ± 0.34 m/year at 0.95 km and 1.40 ± 0.18 m/year at 0.9 km in horizontal and vertical directions, respectively. The reason may be that there is abundant ice at the upper part that is conducive to movement, while the terminus type is land with less ice mass exchange.

Figure 5b and d present that the 3D glacier velocities are the fastest in the mainstream and slower at the two edges along profiles BB′ and DD′ due to increasing frictional drag between the edge and the valley wall. The maximum horizontal velocities are 15.33$\pm$0.11 m/year increasing by 98% and 4.13$\pm$0.01 m/year increasing by 92% from edge to the mainstream of glacier No.1 and glacier No.2, respectively. The ice motion usually increases to maximum value in the mainstream direction and immediately decreases down to zero at the other edge [32]. The same characteristic can also be observed in vertical velocities. On the whole, the errors in the mainstream directions of profile AA′ and CC′ are obvious, but the changes of 3D velocities for profile BB′ and DD′ that are perpendicular to the mainstream directions are relatively stable. The horizontal velocity of glacier No.1 and No.2 clearly and intuitively reflect the spatial distribution of glacier surface, while the change of vertical velocity is more complicated. This is because the vertical velocity not only reflects the thinning or thickening of a glacier but includes the deformation of the glacial moraine.

4.2. Time Series of 3D Glacier Movement

Under the method described above, the time series of 3D glacier movement from January to December 2018 are retrieved. Figure 6 depicts the 3D displacement time series of glacier kinematics at No.1, and Figure 7 for glacier No.2. Color indicates the vertical displacement while the arrow represents the horizontal displacement of the two glaciers. It can be seen that the horizontal displacement of glacier No.1 gradually increased from January to December, showing the longer length arrows in Figure 6. The time series of horizontal motion of a glacier, shown by displacement vectors, is located at the center line of the middle of glacier No.1. The vertical motion shows that the lower part is increasingly downward. However, the upward occurs in the middle part by 77% from January to December due to the mass accumulation at the intersection of two glacier branches.

The horizontal movement of glacier No.2 is mainly distributed in the upper part with the increased length of displacement vectors in Figure 7. The upward movement of glacier No.2 occurs at the eastern part, while the downward movement extensively focusses on the western part and the lower part of glacier in Figure 7.

In order to analyze the temporal deformation characteristics at different positions, we select six points (P1–P6 as shown in Figure 6a and Figure 7a) at front, middle, and the terminus of two glaciers. Figure 8 shows the time series displacement value; to avoid the influence of inconsistent precision, standard deviation (STD) of the 3$\times$3 pixel around each point was estimated for horizontal and vertical direction which is indicated by error bars.

At the front of the glacier No.1 and No.2, the time series of P1 and P4 are close to a linear tendency in horizontal and vertical directions. The horizontal displacement is 6.90$\pm$0.16 m on 25 December while the vertical displacement is –0.58$\pm$0.08 m (see Figure 8a). For the glacier No.2, the horizontal motion at point P4 increased with the maximum cumulative displacement of 11.32$\pm$0.43 m. The vertical deformation is around zero (see Figure 8d). For glacier No.1, the cumulative deformation of point P2 is 9.74$\pm$0.35 m in horizontal direction on 25 December. However, the deformation of the glacier with time is minimal in vertical directions (Figure 8b). Similarly, the dramatic displacement occurs in horizontal direction with 4.50$\pm$0.06 m at point P5 while the downward displacement is –0.51$\pm$0.12 m on 25 December for glacier No.2 (Figure 8e). There are some common characteristics at the termini of the two glaciers, i.e., the vertical change is negative, suggesting an upward phenomenon at the end of the glacier due to mass loss (Figure 8c,f).

4.3. Uncertainty Analysis

Errors in the PO-SBAS method results are mainly caused by offset-tracking and the 3D inversion method. Firstly, there are three main sources of error by means of offset-tracking: SAR image registration, geocoding, and topography [41,42]. In the process of offset-tracking method, we set the spatiotemporal threshold and the correlation coefficient threshold and used different window size to continuously refine the result of offset-tracking. The DEM errors were estimated and removed through Equations (1) and (2). So, the accuracy of co-registration was guaranteed within 0.02 pixels, which is approximately 0.3 m and 0.05 m in azimuth and LOS direction, respectively. The study area is located at towing mountains with hostile environments and blocked roads in winter, which makes it difficult to obtain the real observations by field survey, such as GPS, leveling data, etc. Based on the assumption that the movement of the non-ice region is zero theoretically [43], the STD of velocity was estimated in a stable region in Figure 4c, showing the result were 0.6 m/year and 0.2 m/year in horizontal and vertical direction, respectively. The values in a stable region are all much lower than the average velocities as shown in Figure 4, which indicates the reliability of the result.

Finally, in order to further confirm the quality of our 3D velocities results, we made comparison of glaciers velocity in horizontal direction between our results and the ones from published reference [44]. As the absence for time series velocities in published reference, the velocity in 2018 was used to make comparison. The published results from the reference were obtained by optical feature-tracking method with the resolution of 250 × 250 m. The surface velocity datasets could be downloaded from the website (https://its-live.jpl.nasa.gov/), which contain coordinates and annual surface velocity. Therefore, we used the six points that were discussed in the previous analysis due to the large difference in resolution. The six points can be seen in Figure 6a and Figure 7a and located at the upper, middle, and terminus of the two glaciers, respectively. The detailed process is to sample our surface velocity results to the same resolution as the reference dataset in 2018, comparing the differences and making statistics as shown in Table 2. The maximum difference is 1.73 m/year, and most of them are less than 1.27 m/year, which verifies the feasibility and reliability of our results.

5. Discussion

5.1. The Correlation between Precipitation and Glacier Surface Velocity

It is critical to investigate glacier movement caused by climatic conditions, including precipitation and temperature [44,45]. Fujita et al. [46] revealed that glaciers located in a summer-precipitation climate show higher sensitivity than those in a winter-precipitation climate. Both Shi et al. [47] and Liu et al. [48] proposed a glacier classification system for the Himalayas based upon their continentality. Maussion et al. [49] indicated that glacier velocities are impacted by the uncertainty of snowfall occurrence and magnitude in spring or summer with the Weather Research and Forecasting (WRF) model. They further proposed a new classification using an objective clustering approach based on precipitation seasonality. Precipitation seasonality importantly affected the climatic sensitivities of glacier different types [49]. However, because of the lack of meteorological documentation in the Himalayan climate, it is challenging to understand the precipitation factors.

To further analyze the relationship between precipitation and the glacier flow velocity, we conducted comparison for the glaciers No.1 and No.2. The daily precipitation is obtained from the Global Precipitation Measurement (GPM). The relationship between the flow velocities of two glaciers and the precipitation variation were discussed as shown in Figure 9. Firstly, as glacier No.1 and glacier No.2 are located on the leeward and windward sides of the Indian summer monsoon respectively, the precipitation patterns are significantly different. The peak of rainfall at glacier No.1 was 53 mm on 15 July. June to August was a period with abundant precipitation. The maximum daily rainfall at glacier No.2 was 39 mm on 21 August, and the period of abundant rainfall was in July. In summary, there are different precipitation patterns for these two glaciers. Secondly, the glacier No.1 had undergone apparent two acceleration processes: (I) during the period of 24 March to 29 April, the velocity increased from 5.7$\pm$0.6 m/year to 11.8$\pm$0.4 m/year; (II) during the period of 28 June to 27 August, the velocity increased from 10.7$\pm$0.5 m/year to 12.6$\pm$0.7 m/year, where the maximum rainfall is 53.1 mm on 15 July. As for glacier No.2, the conspicuous acceleration process can be found from 23 May to 10 July with the velocity increased from 5.7$\pm$0.3 m/year to 9.6$\pm$0.5 m/year. On the whole, the surface velocities of the two glaciers have experienced cyclical changes of an acceleration-deceleration-acceleration-deceleration stage. The daily rainfall during the acceleration period is abundant and the rainfall during the deceleration period is poor for glacier No.1. The reason may be that a large amount of surface meltwater and rainfall penetrate into the glacier bedrock in abundant rainfall period, resulting in the decrease of basal frication and increase of flow velocity. However, the horizontal velocity for glacier No.2 does not show a clear association with daily precipitation alone. Therefore, the glacier velocity cannot merely be attributed to daily precipitation. The casual factors of glacier movement are complex and interactional in the central Himalayas.

5.2. 3D Glacier Movement with Seasonal Characteristics

Besides the variations in precipitation, different seasons also exert effects on glacier motion. As the glaciers in our study area are both summer-accumulation glaciers [49], they are more susceptible to increasing temperature. Therefore, the different responses to the same climate circumstance are generated by the different glacier accumulation types.

To make seasonal analysis, we chose four representative classifications, i.e., December to February (DJF), March to May (MAM), June to August (JJA), and September to November (SON). Note that in order to ensure cumulative deformation, “December” here refers to December in 2017. We used PO-SBAS method to connect December of 2017 with January and February of 2018. Figure 10 shows the contributions to the 3D glacier deformation in four different periods, where the seasonal deformation in horizontal and vertical directions is firstly analyzed. It is mostly dry in DJF because the precipitation has a spatially variability with 20–25% [49]. Furthermore, the results illustrate that the displacements are 0.82$\pm$0.10 m and 1.12$\pm$0.05 m in horizontal direction for glacier No.1 and No.2 respectively. In vertical direction, the displacement is –0.07$\pm$0.01 m in downward movement for glacier No.1 and 0.05$\pm$0.01 m in upward motion for glacier No.2. The precipitation in MAM is impacted by the variability of the Indian summer monsoon. It can be detected that the contribution for MAM to glacier No.1 3D movement is greater than that of the glacier No.2. For horizontal and upward direction, the displacements for glacier No.1 are 2.82$\pm$0.10 m and 0.30$\pm$0.07 m, respectively, while the displacements for glacier No.2 are 0.85$\pm$0.08 m and 0.07$\pm$0.04 m respectively. In JJA, the cycle patterns change severely with the highest precipitation in a year, both in India and central Himalayas [49]. The period of JJA for glacier No.1 is impacted by 3.60$\pm$0.24 m and –0.13$\pm$0.02 m in horizontal and downward directions, respectively. In this same period, glacier No.2 is affected by 2.38$\pm$0.30 m and –0.50$\pm$0.05 m in horizontal and downward directions, respectively. It is worth noting that the period in SON is limited by low precipitation. The displacements are 1.51$\pm$0.16 m and 0.28$\pm$0.02 m for glacier No.1 in horizontal and upward respectively, while the displacements are 1.34$\pm$0.08 m and –0.15$\pm$0.02 m for glacier No.2 in horizontal and downward, respectively. Overall, June to August is the most important period for 3D movement of glacier No.2, while for glacier No.1, the four seasonal periods all have a propound influence, especially in vertical direction.

We attribute the difference in two glaciers to their different glacier types, where glacier No.1 is a lake-terminating debris-covered glacier and glacier No.2 is a land-terminating debris-covered glacier. The abundant debris of the glacier No.1 not only provides facilities for the ice cliff and supraglacial ponds development but also creates complicated supraglacial, englacial, and subglacial drainage system components. Further, a supraglacial lake is developed at the lowest part of glacier No.1, which is combinedly produced by wind-driven surface currents and density-driven undercurrents to deliver power to the ice front and lake floor via melting ice [50]. Moreover, lake changes, an important source of GLOFs, had a temporally and spatially diversity effect on glacier movement.

5.3. The Evolution of Debris-Covered Glaciers

Debris-covered glaciers affect glacier length, ice motion, and mass balance in high-mountain regions [2,10,51,52]. The emergence of extensive debris cover means that the linkage of cause and effect is more complicated than on debris-free glaciers [2,4,53,54]. The recession of tongues on large debris-covered glaciers promotes the growth of supraglacial lakes, which brings about an increased risk of GLOFs. Therefore, it is meaningful and important to consider the evolution of debris-covered glaciers. The general conceptual model of the debris-covered glaciers’ evolution was advanced by Been et al. [2]. Three process regimes can be detected in Figure 11.

In regime 1, the surfaces of glaciers are not extensively crevassed and the meltwater can be transferred to their terminus. The ablation phenomenon typically occurs in lower areas and the highest melt rates are located in the mid-ablation areas [2]. The mechanism for glacier No.1 in this study differ from the regime 1 in that there are evident crevasses in the tongue. With the increasing climatic warming, the balance between ice accumulation and ablation rates are altered on debris-covered glacier tongues, causing a transition to some characteristic cliffs. So, a new stage was formed called the regime 2. Consequently, mid-ablation zones are experiencing the increased mass loss and reduced ice accumulation. In turn, this result encourages glacier slowdown and stagnation [10]. Accordingly, the optical image in Figure 4c shows a large moraine-dammed lake and many small concave-up lakes that are facilitated by a reduction of glacier gradient. A lake can form immediately at part of the glacier surface and can expand rapidly by calving and melting above or below the waterline [2]. Therefore, the existence of a supraglacial lake is a dominant sign for determining whether a glacier is in regime 2 or regime 3. These large lakes provide sufficient convenience for potential glacier lake outburst floods. Therefore, the evolution of glacier No.1 or other lake-terminating debris-covered glaciers need to be carefully supervised henceforward.

6. Conclusions

In this study, the time series of the 3D movement of the Cuolangma glaciers in 2018 were estimated from multitrack SAR images by PO-SBAS method, where the one is a lake-terminating debris-covered glacier and the other is a land-terminating debris-covered glacier. Firstly, the multidimensional spatial pattern and temporal evolution of the two glaciers were characterized, which indicated the maximum cumulative displacement of lake-terminating glacier occurs at their middle part, while the maximum cumulative displacement of land-terminating glacier occurs at their upper section.

Then, the correlations analysis between glacier surface velocity and daily precipitation revealed that the surface velocity does not display a clear relationship with daily precipitation alone. In addition, the seasonal displacement shows that June to August exert the biggest contribution on 3D movement of land-terminating glaciers while the seasonal change of lake-terminating glaciers is not evident for their 3D movement. Moreover, the evolution of debris-covered glaciers was evaluated to better understand the formation for a glacial lake, providing reference for the prevention and management of glacial lake outburst floods in the central Himalayas.

More importantly, regarding the absence for vertical component by conventional offset-tracking method, the accuracy and temporal resolution of the glacier motion are limited. Time series of 3D movement for glaciers provide more comprehensive and accurate information for clarifying their movement process. PO-SBAS method for interpreting the 3D movement will be widely applied to more scenarios for different glaciers, such as the evolution of surge of glaciers, the mechanism of glacial lake outburst floods, and the work for disaster deduced by glacier motion. However, the PO-SBAS model could be further optimized with the help of in-situ measurements of glaciers in High Mountain Asia.

Author Contributions

Conceptualization, L.Y., C.Z., and Z.L.; methodology, L.Y. and C.Z.; software, L.Y.; validation, L.Y., C.Z., and C.Y.; data analysis, L.Y., C.Z., and Q.Z.; writing—original draft preparation, L.Y.; writing—review and editing, all the author.; funding acquisition, C.Z., Z.L., C.Y and Q.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work is jointly supported by the National Natural Science Foundation of China NSFC (Grant No. 41874005), China Geological Disaster Investigation Project (DD20190637) and the Fundamental Research Funds for the Central Universities (CHD, 300102269204, 300102269303), Shaanxi Provincial Natural Science Basic Research Project (2019JM-245).

Acknowledgments

The authors would like to thank the European Copernicus for providing the Sentinel-1 SAR data freely. The precipitation data was download from the website https://gpm.nasa.gov/. One-arc-second SRTM DEM is freely downloaded from https://e4ftl01.cr.usgs.gov/MODV6_Dal_D/SRTM/SRTMGL1.003/2000.02.11/. We are also indebted with editor Xiao Becky and three anonymous reviewers for their valuable remarks and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Study area and SAR dataset coverages: (a) Spatial coverage of ascending and descending Sentinel-1 datasets represented on SRTM DEM; (b) The terrain of glacier No.1, (c) The terrain of glacier No.2.

Figure 2. The distributions of temporal and perpendicular baselines of ascending and descending offset pairs used in this study: (a) shows the offset pairs generated by ascending Sentinel-1 data; (b) shows the offset pairs generated by descending Sentinel-1 data.

Figure 3. Flowchart of 3D time series displacement estimation.

Figure 4. The 3D velocities maps and optical images of the two glaciers. Color represents the vertical velocity and arrows represents the horizontal velocity. (a) the 3D velocities for glacier No.1; (b) the 3D velocities for glacier No.2; (c) Landsat-8 image of glacier No.1 in 2018; the red rectangle indicates the stable region; (d) Landsat-8 image of glacier No.2 in 2018.

Figure 5. 3D velocities and elevations of glacier No.1 and glacier No.2 in 2018, where red and green lines represent the horizontal and vertical annual velocity, respectively. Profiles of AA′ and BB′ are shown on Figure 4a and profiles of CC′ and DD′ are shown on Figure 4b. The error bar indicated the standard deviation of the velocity in the area near the selected profiles. (a) Profile A-A′; (b) Profile B-B′; (c) Profile C-C′; (d) Profile D-D′.

Figure 6. Time series of 3D displacement for glacier No.1 from 23 January to 25 December, 2018. Color represents the vertical displacement and arrows represent the horizontal displacement. (a) 20180123–20180429; (b) 20180123–20180604; (c) 20180123–20180803; (d) 20180123–20181002; (e) 20180123–20181225.

Figure 7. Time series of 3D displacement for glacier No.2 from 23 January to 25 December, 2018. Color represents the vertical displacement and arrows represent the horizontal displacement. (a) 20180123–20180429; (b) 20180123–20180604; (c) 20180123–20180803; (d) 20180123–20181002; (e) 20180123–20181225.

Figure 8. 3D deformation time series for points P1-P6. Red square and green triangle represent the horizontal and vertical time series movement, respectively. Locations of P1-P3 are shown on Figure 6a and P4-P6 are shown on Figure 7a. Negative values refer to downward movement. The red error bars indicated the STD of the average 3 × 3 pixels in horizontal directions and the green error bars indicated the errors for vertical direction. (a)–(f) represents the 3D time series displacement of points P1-P6, respectively.

Figure 9. Correlation between glaciers horizontal velocities and precipitation in 2018. (a) point P2 at glacier No.1; (b) point P4 at glacier No.2. Locations of points P2 and P4 are shown on Figure 6 and Figure 7, respectively. The error bars indicated by the STD of the average 3 × 3 pixels in horizontal direction.

Figure 10. Graph of the 3D glacier displacement in four different periods. (a) the horizontal displacement for glacier No.1; (b) the vertical displacement for glacier No.1; (c) the horizontal displacement for glacier No.2; (d) the vertical displacement for glacier No.2; DJF: December, January, February; MAM: March, April, May; JJA: June, July, August; SON: September, October, November. The positive values indicate upward movement while the negative values show downward movement, respectively.

Figure 11. Diagrammatic representation of debris-covered glacier regimes. (a) Regime 1: active ice flow; (b) Regime 2: downwasting ice; (c) Regime 3: calving retreat. Red arrows represent basal flowing, red dotted lines represent the equilibrium line altitudes (ELAS), blue circles represent small lakes, blue shadow represents the basal surface, light orange shadow represents the terrain, orange shadow represents the surface at the terminal of a glacier, the white strips at the terminal in (a,b) represent the flow of glacial melting, gray shading in (c) indicates a glacial lake, and green shadow indicates the surface at the terminal of a lake. (modified from Been et al. [2]).

Table 1

The principal parameters of the synthetic aperture radar (SAR) datasets.

SAR Sensor	Sentinel-1	Sentinel-1
Track No.	85	121
Orbit direction	Ascending	Descending
Incidence angle at scene center (°)	39.327	39.275
Azimuth angle (°)	−12.570	−169.822
Pixel spacing (azimuth $\times \mathrm{range}$)	13.97$\times$2.330	13.97$\times$2.330
Number of scenes	30	30
Acquisition period(yyyymmdd)	20180123–20181225	20180102–20181228

Table 2

Comparison of glaciers flow velocity in horizontal direction between our results and published ones (Dehecq et al. [44]).

Points	Our Results (m/year)	Published Results(m/year)	Difference(m/year)
P1	1.50	2.03	0.53
P2	15.78	15.40	0.38
P3	8.83	9.78	0.95
P4	6.40	8.13	1.73
P5	6.32	6.34	0.02
P6	4.01	5.28	1.27