Seasonal Dynamics of a Temperate Tibetan Glacier Revealed by High-Resolution UAV Photogrammetry and In Situ Measurements

Abstract

The seasonal dynamic changes of Tibetan glaciers have seen little prior investigation, despite the increase in geodetic studies of multi-year changes. This study compares seasonal glacier dynamics (“cold” and “warm” seasons) in the ablation zone of Parlung No. 4 Glacier, a temperate glacier in the monsoon-influenced southeastern Tibetan Plateau, by using repeat unpiloted aerial vehicle (UAV) surveys combined with Structure-from-Motion (SfM) photogrammetry and ground stake measurements. Our results showed that the surveyed ablation zone had a mean change of −2.7 m of ice surface elevation during the period of September 2018 to October 2019 but is characterized by significant seasonal cyclic variations with ice surface elevation lifting (+2.0 m) in the cold season (September 2018 to June 2019) but lowering (−4.7 m) in the warm season (June 2019 to October 2019). Over an annual timescale, surface lowering was greatly suppressed by the resupply of ice from the glacier’s accumulation area—the annual emergence velocity compensates for about 55% of surface ablation in our study area. Cold season emergence velocities (3.0 ± 1.2 m) were ~5-times larger than those observed in the warm season (0.6 ± 1.0 m). Distinct spring precipitation patterns may contribute to these distinct seasonal signals. Such seasonal dynamic conditions are possibly critical for different glacier responses to climate change in this region of the Tibetan Plateau, and perhaps further afield.

1. Introduction

Recent ground-based and large-scale remote sensing studies have revealed that monsoonal glaciers in the southeastern Tibetan Plateau have experienced significant mass loss during the last two decades, exceeding the average for High Mountain Asia [1,2]. Comprehensive studies of glacier dynamics and their controlling factors are essential for understanding glacier response to climate change and for predicting glacier-related disasters, such as glacial lake outburst floods [3,4], ice avalanches [5] and glacial debris flows [6], which are common in this monsoon-dominated glacierized region [7].

Traditional in situ glaciological measurements rely mainly on ablation stakes drilled into the ice at key locations. The resulting point measurements suffer from a relatively low degree of spatial representativeness and are typically labor-intensive to emplace. Satellite remote sensing products can be used to study cryospheric change across a range of spatial (glacier to global scale) and temporal (days to decades) scales, but these data have historically been limited by their relatively coarse spatial resolution as well as extensive data gaps and restricted image quality. This is especially the case for the southeastern Tibetan Plateau, for which the quality of satellite images is generally poor due to pervasive monsoonal cloud cover.

The use of unpiloted aerial vehicles (UAVs) can overcome many of the shortcomings associated with both satellite remote sensing and in situ measurements. UAV surveys are being increasingly used in glaciological studies worldwide [8,9,10,11,12,13]. High-resolution UAV imagery can be used for semi-automated multi-temporal mapping of surface geomorphology and the generation of high-resolution digital elevation models (DEMs) using low-cost Structure-from-Motion (SfM) with multi-view stereo photogrammetry (hereafter “UAV-SfM”) [14,15]. This method is particularly useful for debris-covered glaciers or rugged debris-free glaciers, where the surface texture and topography are suitable for image tie point extraction and matching for photogrammetric reconstruction. However, to date most applications of UAVs for glaciological surveying have relied on the use of GPS-located ground control points (GCP) for georeferencing and to improve the robustness of 3D camera positioning and landscape reconstruction. Using this method, the accuracy of the final product can be significantly influenced by the density and spatial distribution of GCPs [16,17]. The Real Time Kinematic (RTK) and Post-Processed Kinematic (PPK) positioning methods are a promising means of reducing, or entirely removing the requirement for dedicated survey ground control (e.g., in [18,19,20]). The RTK model relies on maintaining a reliable real-time radio link between the UAV and a GNSS base station; however, radio communication problems in mountainous regions limit the feasibility of using the RTK model. In contrast, the PPK method, whereby UAV positional fixes are corrected using base station data after a flight, has been completed and is potentially a more feasible solution for accurate mapping tasks in mountainous regions.

Across High Mountain Asia, high-resolution UAV surveys have been performed predominantly on a few debris-covered glaciers in the Nepalese Himalayas (e.g., in [8,9,21,22]). All of these studies have focused on the UAV data acquired at the beginning (~June) and end (~September) of the ablation season. In addition, to the best of our knowledge, no studies have sought to quantify seasonal glacier dynamic changes, including both “cold” (accumulation) and “warm” (ablation) seasons for Tibetan and Himalayan glaciers. In addition, no studies, to our knowledge, focus on the estimation of ice flux divergence by combining ground-based stake ablation measurements and UAV-based elevation surveys. This is important because the elevation change of the glacier surface is not only controlled by the surface climatic mass balance (ablation and accumulation), but also the ice flux divergence, which results from differences in ice flux into and out of a specific area [23,24]. Examination of seasonal glacier dynamics is critical for predicting future glacier change and understanding their relevant drivers and mechanisms in high-altitude regions.

In this study, high-resolution UAV surveys and ground stake measurements were performed three times between September 2018 and October 2019 on Parlung No. 4 Glacier, a temperate, debris-free glacier in the southeastern Tibetan Plateau (Figure 1). The objectives of the study were three-fold: (i) to evaluate the accuracy of UAV-SfM DEMs derived using PPK-enabled direct georeferencing; (ii) to quantify and compare distributed surface elevation changes and surface velocities in cold and warm seasons; (iii) to quantify the relative contributions of surface mass balance and emergence velocity to the observed surface elevation changes in the ablation zone.

2. Study Site

Parlung No. 4 Glacier is located in the southeastern Tibetan Plateau (Figure 1), which is strongly influenced by the Bay of Bengal vortex and the Indian monsoon system. In the spring season (March to May), warm and humid air from the Bay of Bengal vortex penetrates along the corridor of the Brahmaputra River [25], while the Indian summer monsoon influences the region during June to September [26]. Thus, the monthly precipitation distribution exhibits a double-peak with a substantial amount occurring in both spring and summer along the southern side of Brahmaputra River [25]. The combination of the complex topography and abundant precipitation results in an area of ~10,000 km² of “monsoonal-type” glaciers at high elevations, characterized by high accumulation and high ablation rates [27]. Annual precipitation at the equilibrium line altitude (ELA) reaches about 2.5–3.0 m, and the mean summer air temperature is usually above 1 °C [28]. There are 165 glaciers with areas > 10 km² and the largest glacier has a total area of 204 km². Both in situ measurements and remote sensing studies have indicated that the magnitude of recent ice loss in this region exceeds the average for High Mountain Asia [1,29,30,31,32]. Geodetic measurements indicate a glacier mass losses of ~−0.19 m w.e. a⁻¹ across High Mountain Asia, but with higher mass losses of ~−0.60 m w.e. a⁻¹ in the southeastern Tibetan Plateau during the period from 2000−2018 (e.g., in [32]). The decreased precipitation due to the weakened Indian summer monsoon system and increased temperature have had reduced mass gain and enhanced mass loss in recent decades [1,33].

Parlung No. 4 Glacier is a typical debris-free valley glacier in the southeastern Tibetan Plateau. It has a total glacierized area of 11.7 km² and extends from 5937 m above sea level (asl) at the head of the accumulation zone, to 4610 m asl at the glacier terminus (Figure 1). The glacier has a wide (~2.8 km) accumulation zone and narrow (~0.8 km) ablation zone, and heavy crevassing characterizes the glacier surface from 5000 m asl to 5100 m asl in a ~1 km-long zone which is 4–5 km up-glacier of the terminus (Figure 1). This zone marks the transition between a relatively high (~15°) mean surface slope up-glacier and low (~7°) surface slope down-glacier. The mean ELA was located in ~5400 m asl in recent years [34]. An automatic weather station (AWS) is located at ~4800 m asl. AWS observations have previously revealed that net incoming shortwave radiation controls the surface melting in the ablation zone [35].

3. Data and Methods

3.1. UAV Flights and Data Processing

Optical images of the ablation zone were acquired by UAV surveys conducted on September 23rd 2018, June 7th 2019 and October 7th 2019 (

Table 1. UAV photogrammetric flights used for acquisition of glacier images.

Survey Time	Number of Images	Flight Altitude	Flight Time	Ground Resolution (cm)	Coverage Area (km2)	Flight Time (Beijing Time)
		(m)	(min)
23 September 2018	460	324	83	9	6.3	11:30–12:52
7 June 2019	504	428	84	12	6.4	11:08–12:32
7 October 2019	652	328	149	9	5.4	09:15–10:44

). Both optical images and in situ stake measurements (described below) show that the surveyed area was snow-free in September 2018 and October 2019, but there were a few patches of snow cover around the lateral moraines and upper glacier in June 2019 (Figure A1). We used a senseFly eBee Plus fixed-wing UAV, which has a wingspan of 110 cm and a take-off payload capacity of 1.1 kg. The UAV fuselage contains a GPS receiver, altimeter, wind meter and a senseFly S.O.D.A. (Sensor Optimized for Drone Applications) 20 megapixel RGB digital compact camera that is electronically triggered by the autopilot system to acquire images at specific positions, which are determined by a user-specified forward and side overlap. The UAV has in-built GNSS RTK/PPK functionality with a manufacturer-stated horizontal accuracy of 3 cm and vertical accuracy of 5 cm for RTK positioning. SenseFly’s eMotion 3^® (Cheseaux-sur-Lausanne, Switzerland) flight and data management software was used to plan the UAV flight path prior to deployment, which comprised nadir imagery captured along parallel flight lines with pre-programmed longitudinal and lateral image overlaps of 70% and 65%, respectively. The total surveyed glacierized area was ~2.2 km² (Figure 2).

UAV flights were performed in the morning (Table 1) to minimize aircraft destabilization by katabatic winds. UAV flight lines were perpendicular to the central longitudinal axis of the glacier and the UAV maintained a constant flying height above the glacier surface, which in turn resulted in a constant ground sampling distance (GSD) resolution for each survey: for the September 2018 and October 2019 surveys, the GSD was 9 cm pixel⁻¹, and for June 7th 2019, the GSD was 12 cm pixel⁻¹. A Huaxing A10 GNSS GPS was used as the static base reference station and was deployed on the same stable, off-glacier boulder during the three surveys (Figure 2). Both the drone flight log file (.bb3) and the static GNSS of the reference station were imported into eMotion^® software to re-georeference the images using the PPK method. The software provides a final image geolocation by using the GPS time of the camera trigger marker in the RINEX data of the reference station and linearly interpolating between the two closest points of the 1 Hz record. Post-processed camera positions were attached to EXIF metadata of every geotagged image. The images taken in 2018 and in 2019 were used as inputs to SfM reconstruction using pix4Dmapper software (v. 4.3.31) (Prilly, Switzerland). Post-processed image geotags were incorporated into the software’s bundle adjustment algorithm, where they served the dual purpose of improving the bundle adjustment and georeferencing the 3D point cloud and DEM outputs. The primary outputs for each survey were an orthomosaic and a DEM, which were exported at grid cell resolutions of 0.09 m and 0.12 m, respectively.

DEM differencing was applied to the DEMs in ArcGIS Desktop software (v.10.2) (Redlands, CA, USA) to determine the overall ice height change between successive surveys. The open-source image georectification and feature tracking toolbox ImGRAFT with orientation correlation mode [36] was used to derive maps of distributed surface displacement between successive surveys, using resampled (1 m) greyscale versions of the UAV-SfM orthomosaics, and a moving search window of 10 × 10 pixels as the primary inputs. This resulted in a spatially distributed map of the magnitude and direction of horizontal ice surface displacement vectors, which were used for investigating the seasonal surface movement.

3.2. In-Situ Measurements and Estimation of Emergence Velocity

Five 12 m long PVC stakes (named S1–S5) were inserted into the ice using a Heucke ice drill to measure surface mass balance at different elevations within the surveying area (Figure 2). The stake heights were measured after each UAV survey. The uncertainty was estimated to be within ± 0.1 m due to the manual stake measurements.

The continuity equation for the conservation of mass for ice flow [23], assuming constant density, expresses that ice thickening or thinning ($\mathsf{\delta }\mathrm{H}/\mathsf{\delta }\mathrm{t}$) at any given section on a glacier is the arithmetic product of changes as the following:

$$\frac{\mathsf{\delta }\mathrm{H}}{\mathsf{\delta }\mathrm{t}}=\frac{\mathrm{b}}{\mathsf{\rho }}+\frac{\left({\mathrm{Q}}_{\mathrm{i}\mathrm{n}}-{\mathrm{Q}}_{\mathrm{o}\mathrm{u}\mathrm{t}}\right)}{\mathrm{A}}\approx \Delta \mathrm{h}+\Delta {\mathrm{V}}^{\uparrow }$$

(1)

where b is the surface mass balance (m w.e.);$\mathsf{\rho }$ is the ratio of densities of ice to water;${\mathrm{Q}}_{\mathrm{in}}$ and ${\mathrm{Q}}_{\mathrm{out}}$ are the ice flux entering (in) and leaving (out) the total area (A) of interest. In the ablation zone, the first term in equation ($\mathrm{b}/\mathsf{\rho }$) is the amount of surface lowering ($\Delta \mathrm{h}$) recorded by stakes and the second term in equation $\frac{\left({\mathrm{Q}}_{\mathrm{in}}-{\mathrm{Q}}_{\mathrm{out}}\right)}{\mathrm{A}}$ is also called the emergence velocity when $\Delta {\mathrm{V}}^{\uparrow }$ is positive. Each term in Equation (1) is defined as positive when it is directed upward. If two of these above terms are known, the remaining factor can be calculated using the continuity equation. Under this framework, we subtracted the mean UAV-SfM orthomosaic-derived ice elevation changes ($\mathsf{\delta }\mathrm{H}/\mathsf{\delta }\mathrm{t}$) in 20 m elevation bands from the corresponding ice surface lowing ($\Delta \mathrm{h}$) measured from the ablation stakes to quantify the magnitude of emergence velocity ($\Delta {\mathrm{V}}^{\uparrow }$) in the ablation zone. This method is helpful to elucidate the relative contribution of glacier surface mass balance and emergence velocity to surface elevation change and clarify the differences in these metrics between cold and warm seasons on the southeastern Tibetan Plateau.

4. Results

4.1. Accuracy of DEMs and Orthophotos by UAV Measurements

The relative horizontal accuracy between the DEMs was determined by manually comparing the displacements of a total of 43 visible large stable boulders in the proglacial area and along the lateral moraines of the glacier from the three UAV-SfM orthophotos (Figure 2). The statistical results show that the mean horizontal displacements (+1σ) are 0.10 ± 0.04 m for the cold season, 0.12 ± 0.05 m for the warm season and 0.13 ± 0.06 m at the annual scale, respectively (Figure 3). The horizontal accuracy is, therefore, estimated to correspond to less than ~1.5 GSD.

The relative vertical accuracy between the DEMs was firstly assessed by comparing the elevation difference in rectangular domains in the proglacial region (Region I) where we assume the topography has remained stable over an annual timescale (Figure 2). We also selected another rectangular region (Region II) along the true left moraine for testing the vertical accuracy in ice-marginal regions, although this region has less long-term stability due to the influence of the freezing and thawing process and steep moraine slopes, which can destabilize.

The histogram of height differences between three DEMs in Region I shows that the most of the vertical difference ranges from −0.2 m to 0.2 m (Figure 4). The mean DEM difference was −0.03 ± 0.11 m for the cold season, 0.02 ± 0.09 m for the warm season, and −0.01 ± 0.1 m for the annual scale, respectively. The median values were −0.03 m for the cold season, 0.02 m for the warm season and −0.01 m at the annual scale. For Region II, the accuracy is slightly decreased but still with decimeter accuracy, with the mean DEM difference of −0.22 ± 0.32 m for the cold season, −0.10 ± 0.28 m for the warm season and 0.13 ± 0.16 m for the annual scale, respectively. Such examination on two different regions revealed that the vertical uncertainties derived from the comparison of DEMs should be limited to decimeter-accuracy.

Such decimeter-scale horizontal and vertical uncertainties between the three orthophotos and DEMs confirm the ability of UAV-SfM to accurately capture elevation changes and surface velocity without any GCPs in highly glacierized mountainous regions, particularly where the magnitude of ice surface elevation change exceeds the estimated errors by an order of magnitude or more. These findings support the use of UAV-SfM photogrammetry using PPK-enabled direct georeferencing to facilitate the measurement of seasonal-annual, glacier-scale mass loss and velocity measurements.

4.2. Surface Elevation Changes at Different Seasons

The spatial distributions of ice surface elevation changes for different seasons are shown in Figure 5 and summarized in Figure 6. In the cold season (October 2018 to June 2019, a total of 258 days) we observed isolated regions of surface lowering (max. 4 m) near the glacier terminus, as well as high variability of surface lowing and uplift in the upper zone of the heavy crevassing region, which are likely the product of horizontal ice flow and serac collapse. Aggregated over the entire surveying area, the mean cold season elevation change was 2.0 m (median of 1.6 m and standard deviation of 2.2 m), with 86% of the elevation change falling in the range 0 m to 5 m (Figure 6a). The mean elevation changes within 20 m elevation bands show an altitudinal gradient of ~0.88 m 100 m⁻¹, with a minimum elevation change of 0.55 ± 0.82 m at the glacier terminus to a maximum of 4.20 ± 5.7 m near the icefall zone at ~5000 m (Figure 7a). Above ~5000 m asl, elevation change shows a decreasing trend (−0.94 m 100 m⁻¹) and the elevation change reduces to ~3.0 m in the uppermost region (5090 m asl). In the cold season, the total net volume increase in all pixels is + 4.4 × 10⁶ m³, corresponding to 8 mm day⁻¹ of ice volume increase.

In the warm/ablation season (122 days) we observed widespread surface lowering (Figure 5b). The mean elevation difference was −4.7 m (median −4.8 m and standard deviation of 1.9 m) between June and October 2019 (Figure 6b), corresponding to 39 mm day⁻¹ of ice surface lowering, or about 0.9 m³ s⁻¹ of melted water in the warm season (122 days), assuming an ice density of 900 kg m⁻³. Ice surface elevation change for 70% of the survey area falls in the range −4.0 m and −6.0 m (Figure 6b). Only a small percentage of pixels (1.5%) near the icefall region shows the elevation increase due to the horizontal displacement of ice cliffs and ogives. Disregarding the icefall region, where loss-gain elevation change pairings complicate analysis, the mean ice surface elevation change exhibits a weak decreasing trend (−0.2 m 100 m⁻¹) with increasing altitude (Figure 7b). Similar to the data in the cold season, there is less spatial variability below 4900 m asl, but an increased degree of variability in the icefall region and then back to reduced variability above the ice cliff region.

From September 2018 to October 2019 the mean ice surface elevation change in the survey area was −2.7 m (median −3.0 m and standard deviation of 2.4 m; Figure 6c), which corresponds to a total net volumetric loss of 5.9 × 10⁶ m³. The mean elevation change exhibits a steep trend from the mean value of −5.0 m at the glacier terminus to about −0.8 m at the upper part (Figure 7c). The vertical gradient of elevation change became as steep as −1.0 m/100 m. At the annual scale the altitudinal ice surface elevation change gradient (Figure 7c) is the sum of seasonal gradients (Figure 7a,b). Overall, the annual gradient of elevation change is mainly controlled by the gradient in the cold season.

4.3. Surface Velocity

We derived distributed surface velocity fields for the glacier surface fields for areas below 4850 m asl. Above this elevation, complex surface deformation, mostly in the heavily crevassed region, hindered feature tracking. Velocity and flow direction vectors are shown in Figure 8. Overall, the horizontal flow displacement of the glacier is characterized by down-glacier deceleration in both the cold and warm seasons, in line with observations of velocity fields on other mountain glaciers.

During the cold season (258 days) the mean horizontal displacement rate is ~0.06 ± 0.03 m day⁻¹, with standard deviation, maximum and minimum surface velocities of 0.03 m day⁻¹, 0.13 m day⁻¹ and 0.02 m day⁻¹, respectively. For the warm season (122 days), the mean surface velocity is faster than that in the cold season (+40%); the mean daily velocity is 0.09 m day⁻¹ with a standard deviation, maximum and minimum of 0.03 m day⁻¹, 0.18 m/day and 0.04 m day⁻¹, respectively. The speed-up during the warm season agrees with other studies in the Himalayas [9,37]. This speed-up is likely caused by enhanced basal sliding due to increased lubrication of the bed and the development of a distributed subglacial drainage system early in the ablation seasons. Over our total study period (380 days) the mean surface movement of all pixels is 25.9 ± 10.3 m. Adjusted to an annual velocity, this equates to 23.9 m a⁻¹ or a daily velocity of 0.07 m day⁻¹. Annual displacements reach a maximum of 52.8 m a⁻¹ at 4850 m asl and a minimum of 10.6 m a⁻¹ at the glacier terminus. By manually tracking the horizontal displacement of a cluster of surface debris in the ice cliff region using the UAV-SfM orthomosaics, we estimate a flow velocity of ~78 m a⁻¹ in this region. Given that this debris cluster was offset from the glacier’s central flow line, where surface flow velocities are fastest, the surface velocity in the icefall region may approach or exceed 100 m a⁻¹.

4.4. Surface Ablation and Emergence Velocity

As mentioned in Section 3.2, ice thickening or thinning ($\Delta$H) at any given section on the ablation zone is the arithmetic sum of height changes by surface lowering ($\Delta \mathrm{h}$) and the emergence velocity ($\Delta {\mathrm{V}}^{\uparrow })$. In situ measurements of $\Delta \mathrm{h}$ were available by using five stakes along the central flowline of Parlung No. 4 Glacier. Assuming the same amount of ablation in the same elevation bands of stakes (±10 m), the sectional emergence velocity ($\Delta {\mathrm{V}}^{\uparrow })$ could be roughly estimated by subtracting $\Delta$H of 20 m-band DEM change from the in situ measured stake $\Delta \mathrm{h}$.

Table 2. In-situ measured surface lowering ($\Delta \mathrm{h}$), UAV-SfM orthomosaic-derived surface elevation change ($\Delta$H) and emergence velocity ($\Delta {\mathrm{V}}^{\uparrow }$) near the five stakes on the ablation zone of Parlung No. 4 Glacier over seasonal and annual timescales.

Stake	Elevation (m asl)	Cold Season (258 days)			Warm Season (122 days)			Annual (380 days)
		$\mathbf{\Delta }\mathit{H}\left(\mathbf{m}\right)$	$\mathbf{\Delta }\mathit{h}\left(\mathbf{m}\right)$	$\mathbf{\Delta }{\mathit{V}}^{\mathbf{\uparrow }} \left(\mathbf{m}\right)$	$\mathbf{\Delta }\mathit{H}\left(\mathbf{m}\right)$	$\mathbf{\Delta }\mathit{h}\left(\mathbf{m}\right)$	$\mathbf{\Delta }{\mathit{V}}^{\mathbf{\uparrow }}\left(\mathbf{m}\right)$	$\mathbf{\Delta }\mathit{H}\left(\mathbf{m}\right)$	$\mathbf{\Delta }\mathit{h}\left(\mathbf{m}\right)$	$\mathbf{\Delta }{\mathit{V}}^{\mathbf{\uparrow }} \left(\mathbf{m}\right)$
S1	4660	1.2 ± 0.9	−2.2 ± 0.1	3.4 ± 1.0	−5.1 ± 0.5	−6.5 ± 0.1	1.4 ± 0.6	−4.5 ± 0.7	−8.7 ± 0.1	4.2 ± 0.8
S2	4720	1.3 ± 0.9	−1.3 ± 0.1	2.7 ± 1.0	−4.7 ± 0.5	−5.4 ± 0.1	0.8 ± 0.6	−3.3 ± 0.8	−6.8 ± 0.1	3.5 ± 0.9
S3	4760	1.5 ± 1.0	−1.1 ± 0.1	2.6 ± 1.1	−4.8 ± 0.8	−5.1 ± 0.1	0.3 ± 0.9	−2.9 ± 0.8	−6.2 ± 0.1	3.2 ± 0.9
S4	4780	1.8 ± 1.1	−1.1 ± 0.1	2.8 ± 1.2	−4.9 ± 0.9	−5.0 ± 0.1	0.1 ±1.0	−3.0 ± 0.9	−6.0 ± 0.1	3.0 ± 1.0
S5	4870	2.7 ± 1.7	−1.0 ± 0.1	3.7 ± 1.8	−4.3 ± 1.6	−4.9 ± 0.1	0.5 ± 1.7	−1.3 ± 1.7	−5.9 ± 0.1	4.5 ± 1.8

shows the magnitude of each component across seasonal and annual scales at the elevational bands of each five stakes.

Our in situ stake measurements for the cold season show that a certain amount of surface ablation has already occurred before the onset of the Indian summer monsoon with the maximum of −2.2 ± 0.1 m near the glacier terminus and the minimum of −1.0 ± 0.1 m at the uppermost stake (S5). The UAV-SfM orthomosaic-derived $\Delta$H at the same elevation bands within the stakes ranges from 1.2 ± 0.9 m to 2.7 ± 1.7 m. The derived sectional $\Delta {\mathrm{V}}^{\uparrow }$ ranges from 2.6 ± 1.1 m to 3.7 ± 1.8 m. In the warm season, the surface ablation decreases from the glacier terminus to the upper glacierized region, with the maximum value at glacier terminus of −6.5 ± 0.1 m. The $\Delta {\mathrm{V}}^{\uparrow }$ ranges from 0.1 ±1.0 m to 1.4 ± 0.6 m, which slightly compensates the surface ablation. At the annual scale, the magnitude of surface ablation can reach between approximately −6 and −8 m and the $\Delta {\mathrm{V}}^{\uparrow }$ ranges from 3~4.5 m, which contributes to a surface lowering of between approximately −1.3 and −1.5 m near the stakes on the ablation zone.

5. Discussion

5.1. The Contributions of Surface Ablation and Emergence Velocity to Surface Elevation Change

In order to clarify the overall contribution of each component in the ablation zone, the values of each component at five stake locations was averaged to represent the mean condition in different seasons (Figure 9). Regarding the uniform distribution of stakes along the central flowline (Figure 2), we assume that such a simple spatial average can be taken to represent conditions across the wider ablation zone. Using the UAV-SfM orthomosaic-derived elevation change data as a result, the averaged δH/δt among the five stakes (Table 2) agree well with the values averaged by all regional pixels below the elevation of 4850 m (S5), with 1.7 ± 1.1 m vs. 1.5 ± 1.0 m in the cold season, −4.8 ± 0.9 m vs. −4.8 ± 0.9 m in the warm season and −3.0 ± 0.9 m vs. −3.3 ± 1.2 m across an annual timeframe, respectively. The mean emergence for five stakes was therefore estimated to be 3.0 ± 1.2 m. These data clearly demonstrate that considerable ice supply from the upper glacierized region contributes to the downstream elevation increase for the vast majority of the survey area in the cold season (Figure 5a). A larger gradient serves to suppress ice surface lowering at higher altitudes because of strong ice fluxes from the upper glacier, whilst allowing high rates of surface lowering to be maintained in the ablation season at lower elevations (Figure 7).

In contrast, the whole ablation zone suffers dramatic surface ablation in the warm season, and this leads to the surface of the glacier lowering considerably (Figure 9). The mean UAV-SfM orthomosaic-derived elevation change over the survey area (−4.8 ± 0.9 m) was dominated by surface ice ablation (−5.4 ± 0.1 m) in this period, whilst the relative contribution of the emergence velocity was only limited to be 0.6 ± 1.0 m, which was one fifth of the value recorded for the cold season. At the annual scale, the surface ablation lowering signal (−6.7 ± 0.1 m) was significantly suppressed by the dynamic thickening (3.7 ± 0.9 m), which equated to approximately 55% of annual surface ablation for our domain, and which mainly occurred during the cold season. These findings shed light on the distinct seasonal glacier dynamics and surface ablation conditions in both the cold and warm season for this temperate glacier and such phenomena could be happening elsewhere for similar monsoon-dominated temperate glaciers in the southeastern Tibetan Plateau.

The annual emergence velocity for the ablation zone of Parlung No. 4 Glacier (~3.7 ± 0.9 m a⁻¹) is obviously larger than that in the Qiyi glacier (0.3–0.4 m a⁻¹) on the northeastern Tibetan Plateau [38] and Urumqi Glacier No. 1 (~0.2~0.3 m a⁻¹, Xu Chunhai, pers. comm.) in the Tien Shan [39] and debris-covered Lirung Glacier (0–0.35 m a⁻¹) in the Nepalese Himalayas [40], but is similar to the values for the Mer de Glacier (~4.0 m a⁻¹) in the French Alps [24] and for the central zone of the Khumbu Glacier (5–6 m a⁻¹) in the Nepalese Himalayas [41]. The temperate Parlung No. 4 Glacier is located in the monsoon-dominated southeastern Tibetan Plateau region, characterized with large mass accumulation and strong surface ablation. Our results indicate that seasonal dynamic thickening, particularly the change in the cold season, should be carefully considered when predicting the mass balance response of Parlung No. 4 Glacier, or similar large glaciers in the southeastern Tibetan Plateau, to future climate change.

5.2. The Possible Mechanism of Dynamic Thickening in Cold Season and the Impact on Regional Glacier Changes

Our above analysis indicated that seasonal ice dynamics lead to considerable glacier thickness change and are a significant components of mass change for the ablation zone. We attribute the comparatively high rates of dynamic thickening in the cold season for Parlung No. 4 Glacier to several possible factors. Firstly, temperate glaciers in this region receive substantial snowfall in the winter-spring season—an ice core drilled at 5500 m asl on the study glacier revealed that the annual snow accumulation is 2.45 m w.e. a⁻¹ during the period 1998 to 2005 [42] and about 2.0 m a⁻¹ for neighboring Zuoqiupu Glacier (5600 m asl) between 1956 and 2006 [43]. Both ground-based meteorological records and atmospheric modelling data from the High Asia Refined (HAR) analysis show that monthly precipitation exhibits a double-peak pattern with a considerable amount falling in both the spring and warm seasons [25,44]. Energy-based mass balance models have revealed that temperate glaciers in this region receive the majority (~50%) of accumulative snowfall in the spring season (March-May), and has thus been classified as a “spring-accumulation type” glacier [25] which is distinct from more typical “summer-accumulation-type” glaciers on the Tibetan Plateau [45]. The abundant accumulation at high elevations in the winter–spring season possibly contributes to high emergence velocity in the ablation zone in the cold season (Figure 9). Secondly, the warm season surface melt is so high that it would actually reduce the supply of ice down-glacier in the warm season, which in turn highlights the importance of the emergence velocity in the cold season. The magnitude of surface ablation at the elevation of 5200 m asl reached ~3 m w.e. [25]. In particular, for the ice fall region (Figure 1), the abundant ice crevasse will greatly increase the exposed area for surface ablation. The total ice fluxes through the section of crevasse region will be greatly reduced, which will possibly influence the downward surface elevation change in the warm season. Thirdly, the glacier geometry also appears to be a potential factor because the Parlung No. 4 Glacier has a wide accumulation zone, providing favorable geometry for high ice emergence velocities in the narrower ablation zone (Figure 2). The mean width of the ablation zone below the elevation of 5000 m is 0.75 km, while the accumulation zone is 2.8 km wide. The glacier width narrows at 5100~5000 m asl, which coincides with the location of the 1.2 km-long crevassed region due to ice flux surge. We acknowledge that these three factors are the possible explanations for the seasonal difference of glacier dynamic change for the Parlung No. 4 Glacier in the southeastern Tibetan Plateau and more detailed process analyses are required to explore the possible mechanism.

Such dynamic thickening in the cold season could influence the glacier response to climate change in the southeastern Tibetan Plateau, partly because of their significant role of decelerating the elevation lowing and maintaining the glacier geometry for the ablation zone. As described above, the emergence velocity contributes as much as ~55% compensation for surface ablation over the lower glacier, which will play an important cumulative role in decades or longer time scales. Meanwhile, the comparison of terminus positions between the three orthomosaics showed that the Parlung No. 4 Glacier may undergo a seasonal cycle in terminus position, characterized by weak advance (0 m~10 m) in the cold season but remarkable retreat (between approximately −5 and −26 m) in the warm season.

Recent large-scale geodetic results reveal a clear contrast between moderate elevation decrease in the south and dramatic elevation decease in the north of southeastern Tibetan Plateau [46]. Previous energy-based mass balance modelling has revealed that the surface mass balance of spring-accumulation type glaciers are less sensitive to climate change when compared to summer-accumulation type glaciers [25]. Such contrasting sensitivities of surface mass balance between different types of glaciers were therefore adopted to explain such spatial difference on mass loss pattern [46]. Our findings provide an additional perspective to explain such contrasting spatial variability of glacier elevation changes in this monsoon-dominated region. We suggest that there are possibly large differences in emergence velocity in the cold season between spring–winter accumulation-type glaciers and summer accumulation-type glaciers. To our knowledge, there is no detailed study of seasonal ice dynamics on the summer accumulation-type glaciers on the Tibetan Plateau and Himalayas, and this is an obvious area for future research to focus on. Future detailed studies of seasonal dynamic conditions between different accumulation-type glaciers would be helpful for understanding the drivers, mechanisms and long-term glacier responses to climate change, and implications for glacier-related hazards [4,5].

Temperate glaciers in the monsoon-influenced southeastern Tibetan Plateau have experienced overall dramatic mass loss and unfavorable climatic combinations in the last two decades [1,29,32,33]. Glacier surface velocities have slowed dramatically. The mean rate of slowdown across the wider southeastern Tibetan Plateau is −0.64 ± 0.02 m a⁻¹, representing a ~49% decrease during the period 2000–2017 [47]. If the Parlung No. 4 Glacier exhibits a similar trend, we would expect this to manifest as reduced ice flux into the ablation zone and a weakened relative contribution (percentage) of emergence velocity to ice surface elevation change, resulting in enhanced surface lowering.

6. Conclusions

We have quantified the seasonal glacier dynamics in the ablation zone of the Parlung No. 4 Glacier, a temperate glacier in the southeastern Tibetan Plateau, using repeat high-resolution UAV-SfM photogrammetry and field stake measurements. Our UAV-SfM method employed direct georeferencing using PPK positioning to produce decimeter-accuracy DEM and orthomosaic products, which allowed the robust detection of ice surface elevation change and movements.

DEM differencing revealed an increase in glacier ice surface elevation (mean 2.0 ± 0.15 m) during the cold season (September–June) followed by a decrease in ice surface elevation (mean −4.7 ± 0.15 m) during the warm season (June–September). The results of feature-tracking revealed higher mean surface velocities in the warm season compared to the cold season (~0.09 ± 0.03 m day⁻¹ vs. ~0.06 ± 0.03 m day⁻¹ or +40%), with the highest flow velocities (52.8 m a⁻¹) at the up-glacier limit of our survey area, and the lowest (10.6 m a⁻¹) at the glacier terminus. In addition, corresponding in situ and geodetic measurements of ice surface height change allowed us to help us to isolate the relative contribution of emergence velocity at seasonal and annual timescales. We found that the positive surface height change in the cold season is caused by high emergence velocities (mean 3.0 ± 1.2 m), which exceed ablation during this period, whereas magnitude of ablation exceeds emergence velocity in the warm season (−5.4 ± 0.1 m vs. 0.6 ± 1.0 m). Over an annual timescale, the surface ablation lowering signal in the ablation zone was significantly suppressed by dynamic thickening, and which offsets approximately 55% of the surface lowering due to surface ablation.

Such seasonally distinct glacier dynamics could contribute to the spatial variability of glacier mass loss in the southeastern Tibetan Plateau. Significant accumulation in the cold season can suppress ice surface lowering in the ablation zone to the point where the glacier exhibits very little annual surface lowering despite high summer ablation. These results highlight the importance of careful assessments of ice flux divergence in addition to glacier surface elevation changes. On-going, regular UAV and ground truth measurements are essential for understanding the drivers and mechanisms of glacier response to climate change for spring- and summer accumulation-type glaciers in this and other monsoon-dominated regions.