Characterizing Changes in Geometry and Flow Speeds of Land- and Lake-Terminating Glaciers at the Headwaters of Yarlung Zangbo River, Western Himalaya

1. Introduction

The Himalayas, spanning over 2400 km from west to east, are profoundly influenced by the Indian monsoon and the mid-latitude westerlies [1]. This region stands as the most heavily glacierized area in the low latitudes of the Earth, harboring the largest ice volume outside the polar regions [2]. Himalayan glaciers serve as vital water resources for densely populated regions in South Asia and the Qinghai-Tibet Plateau [3,4], playing a critical role in downstream river runoff and modulating seasonal variability [5,6,7]. Furthermore, meltwater from these glaciers sustains essential ecosystem services crucial for socio-economic development in mountainous areas, including agricultural irrigation, hydropower generation, and domestic water supply [8,9,10].

Himalayan glaciers have been experiencing mass loss over several decades [11,12], with a twofold increase in loss rates observed since the 21st century [13,14,15]. This trend poses challenges to regional water security and increases risks of glacial hazards (e.g., glacial lake outburst floods and ice avalanches) [16,17]. Due to diverse climate patterns, the glacier mass loss rates across the Himalayas exhibit high spatial heterogeneity [18,19]. However, large variations in mass balances among individual glaciers within the region were also documented [13,20,21], suggesting that there are factors partly independent of climate affecting the processes of mass changes. One important factor is the rapid growth of proglacial lakes during recent decades in the Himalayas [22,23,24,25,26]. The proglacial lake expansion can affect the glacier ablation through frontal calving and subaqueous melt [23,26], and it can also influence the longitudinal coupling stress and basal slip, thereby impacting ice dynamics and mass balance [27,28]. Remote sensing studies have revealed the contrasting patterns of mass and velocity changes between lake- and land-terminating glaciers in the Himalayas [29,30,31]. Generally, the lake-terminating glaciers lose more mass and flow faster than land-terminating glaciers [23,30,31]. The ice-contact lakes complicate the physical processes of glacier mass balance and dynamics, and raise modeling challenges for future glacier projection [32]. It is therefore crucial to investigate the multi-decadal changes in both mass and velocity of Himalayan glaciers, which would provide further insights into their nonlinear response to climate change and lake evolution.

Previous studies have revealed statistically robust differences in surface elevation and velocity changes between lake- and land-terminating glaciers [24,25,29,30,31]. However, the majority of these investigations have been concentrated in the eastern Himalayas. A recent investigation by Scoffield et al. suggests that these differences are insignificant in the western Himalaya due to earlier evolutionary stages of ice-contact lakes in this area compared to those in the eastern Himalaya [33]. Nonetheless, the number of ice-contact lakes in the western Himalaya is projected to substantially increase by 2100, indicating that future ice loss may become pronounced due to lake expansion [34]. Observations of glacier changes in the western Himalaya can help us understand their dynamic behaviors, especially those of lake-terminating glaciers at the early development stage. Elucidating the key processes in glacier-lake interactions will be crucial for developing process-based models and improving large-scale projections.

Space-based measurements have been extensively used to monitor Himalayan glacier changes at various geographic scales. However, the accuracy of these satellite-derived results is often limited by the lack of field data. Due to logistical challenges, field measurements in the Himalayas are restricted to only a few glaciers [35], and are particularly scarce on lake-terminating glaciers [30], thereby constraining the calibration and validation processes. While numerical glacier evolution models have been utilized to project future glacier changes [36,37,38], few of these models explicitly incorporate the dynamics of lake-terminating glaciers. Consequently, they are unable to fully capture the impacts of proglacial lakes on glacier evolution. Numerical ice flow modeling is of importance for understanding the physical mechanisms underlying lake-terminating glacier changes and provides valuable insights into stress balance and flow regime changes. Current studies often employ simplified analysis using shallow-ice approximation models [31]. However, a more sophisticated model, such as a higher-order model, is essential for comprehensively understanding the flow regimes of glaciers with different termini types.

This study investigates the long-term changes in land- and lake-terminating glaciers at the headwaters of the Yarlung Zangbo River (YZR) in the western Himalayas. Our aim is to reveal the temporal impacts of proglacial lakes on glacier mass and velocity changes, and to contrast the distinct dynamic behaviors of these two types of glaciers within the same region. To achieve this, we use multi-source digital elevation models (DEMs) to calculate multi-decadal changes in glacier surface elevation. We validate the satellite-derived High Mountain Asia (HMA) glacier velocity product using field measurements from our study region. We then assess glacier surface velocity changes using this satellite-derived product. Additionally, we track the terminus changes of representative glaciers to examine the interaction between glaciers and lakes, thereby enhancing our understanding of glacier mass and velocity changes. Finally, we employ a higher-order ice flow model to investigate changes in glacier flow speed regimes by assimilating satellite-derived velocities.

2. Study Area

Our study area is located on the northern slope of the Himalayas and serves as the headwaters of the YZR. The study area is defined based on the geographical extent of topographical maps (see Section 3.2), spanning $\sim {30}^{\circ }{0}^{\prime }-{30}^{\circ }{25}^{\prime }$N and $\sim {82}^{\circ }{0}^{\prime }-{82}^{\circ }{40}^{\prime }$E (Figure 1). The average altitude of this area exceeds 5000 meters a.s.l., and the climate is semi-arid with dry and cold characteristics. This aridity is due to the Himalayas blocking moisture from the Indian monsoon and the Karakoram range blocking moisture from the westerlies. The average annual air temperature in this area is approximately $3.5°\mathrm{C}$, and the annual total precipitation is about 157 mm, concentrated in January-April and June-August [39].

There are about 25 glaciers in this region with areas larger than 1 ${\mathrm{km}}^2$. These glaciers are classified as subcontinental type and flow northwards onto the southwestern Tibetan Plateau. According to the second Chinese glacier inventory (CGI2), almost all glaciers in this region are debris-free and primarily land-terminating. However, a few glaciers terminate in proglacial lakes, making them lake-terminating (Figure 1). This region, therefore, provides an ideal setting to compare the spatio-temporal changes between land-terminating and lake-terminating glaciers.

For detailed analysis of changes in mass, velocity, and flow regime, we select four glaciers. Two are land-terminating glaciers (GLIMS IDs: G082213E30171N and G082128E30241N), which we refer to as N171 Glacier and N241 Glacier, respectively. The other two are lake-terminating glaciers (GLIMS IDs: G082163E30214N and G082215E30085N), known as Jiemayangzong Glacier and Asejiaguo Glacier, which we refer to as JMYZ Glacier and ASJG Glacier. The sizes of these glaciers are relatively similar.

3. Data and Methods

3.1. Glacier Outlines and Types

We used the glacier outlines from the CGI2, which was derived from the Landsat images [40]. The glaciers in our study region were classified as either land-terminating or lake-terminating types. We identified the lake-terminating glaciers by combining data from CGI2 and the HMA glacial lake inventory [41], the latter of which provides glacial lake boundaries for 1990 and 2018. Glaciers were classified as lake-terminating if contact between the terminus and a proglacial lake was evident in both 1990 and 2018. The identified lake-terminating glaciers were further verified through manual interpretation of Landsat images obtained in 1990 and 2018 to ensure that the glaciers terminated into the lakes. Based on the CGI2 glacier outlines, we manually digitalized the glacier terminus between 1990 and 2019 with a time interval of approximately 5 years (see Figure 7).

3.2. Digital Elevation Models

In this study, we used multiple DEMs to assess glacier surface elevation changes, including TOPO DEM, SRTM DEM, and TanDEM-X DEM. Six TOPO DEMs were digitized from $1:\mathrm{50,000}$ scale topographic maps produced by the Chinese Military Geodetic Service in October 1974, based on aerial stereo pairs. Contours and elevation points from these maps were used to construct a triangulated irregular network, which was then used to generate the TOPO DEMs. The spatial reference systems of the maps were transformed from the Beijing Geodetic Coordinate System 1954 and Yellow Sea 1956 datum to the World Geodetic System 1984 (WGS84) and Earth Gravity Model 1996 (EGM96) using a seven-parameter transformation method.

The SRTM DEM, derived from C-band radar images with an original resolution of 3 arc seconds (approximately 90 m), was collected by NASA in February 2000. Due to its continuous spatial coverage, the C-band DEM was selected for this study. The SRTM DEM has been widely used in geodetic glacier mass balance estimation [42,43,44]. Its horizontal and vertical datums are WGS84 and EGM96, respectively.

We also used the TanDEM-X DEM, a global product with a 90 m spatial resolution. The TanDEM-X mission aimed to generate accurate 3D land surface models. It used SAR instruments aboard the TanDEM-X and TerraSAR-X satellites. DEM data for this region was obtained from multiple radar images taken around 2017. The absolute horizontal and vertical errors of the TanDEM-X DEM were both less than 10 m. The WGS84 ellipsoid heights of the TanDEM-X DEM were converted to EGM96 geoid heights using the MATLAB geoidheight function.

3.3. Glacier Mass Changes

In this study, we calculated glacier mass changes for the period 1974–2000 by differencing the TOPO and SRTM DEMs. To extend these results to the past two decades, we incorporated an open-access global glacier elevation change product derived from time series of DEMs with a spatial resolution of 100 m [45]. We refer to this as the Hugonnet dataset. For our analysis, we utilized the elevation change maps from the Hugonnet dataset for the periods 2000–2010, 2010–2020, and 2000–2020. The methods used for calculating glacier mass changes and the associated uncertainty assessment are described below.

Co-registration of DEMs is a prerequisite before performing the differencing. First, both DEMs were resampled to a 90 m spatial resolution and reprojected to the same projection system (WGS84, UTM 44N). We co-registered the TOPO DEM to the SRTM DEM over stable terrain following the Nuth and Kääb method [46]. Stable terrain pixels were selected outside the glaciers with slopes between $5^{\circ }$ and ${45}^{\circ }$. This method iteratively computes the shift vector by relating elevation differences to terrain slope and aspect.

The SRTM DEM, derived from C-band radar images, can be biased due to radar wave penetration into ice and dry snow [47]. To correct this bias, we applied a linear model that estimates radar penetration depth ($d_p$) based on altitude (z) [48]. This model assumes a linear relationship between altitude and radar penetration depth across a $1^{\circ }\times 1^{\circ }$ grid in HMA. For our study area, the local relationship $d_p=0.00355z-18.2$ was used. The average penetration depth in this region was estimated to be around $2.0$ m.

After performing DEM differencing, we removed outliers where the elevation change ($\Delta h$) exceeded 100 m or the slope was steeper than ${45}^{\circ }$. Additionally, we filtered out pixels with elevation differences outside the $31.7\%$ to $68.3\%$ quantile range, following the approach recommended by [49]. Finally, we converted glacier elevation changes to mass changes ($\Delta M$) using the following equation:

$$\Delta M={\rho }_c{A}_g\Delta h$$

(1)

where ${\rho }_c=850$ kg ${\mathrm{m}}^{-3}$ is the density conversion factor from volume to mass [50] and $A_g$ is the glacier area.

We calculated the uncertainties in mass change by considering four main sources of error: uncertainty in elevation change (${\sigma }_{\Delta h}$), uncertainty in glacier area (${\sigma }_{A_g}$), uncertainty in the density conversion factor (${\sigma }_{{\rho }_c}$), and uncertainty in the C-band penetration depth (${\sigma }_{d_p}$). Assuming these four error components are independent, the total uncertainty in mass change, ${\sigma }_{\Delta M}$, is given by:

$${\sigma }_{\Delta M}=\left|\Delta M\right|\sqrt{{\left({\displaystyle \frac{{\sigma }_{\Delta h}}{\Delta h}}\right)}^2+{\left({\displaystyle \frac{{\sigma }_{A_g}}{A_g}}\right)}^2+{\left({\displaystyle \frac{{\sigma }_{{\rho }_c}}{{\rho }_c}}\right)}^2+{\left({\displaystyle \frac{{\sigma }_{d_p}}{\Delta h}}\right)}^2}.$$

(2)

We quantified ${\sigma }_{\Delta h}$ using the normalized median absolute deviation (NMAD) of the elevation change estimates over stable terrain [51], which is given by:

$${\sigma }_{\Delta h}={\displaystyle \frac{{\sigma }_p}{\sqrt{N_{\mathrm{eff}}}}},$$

(3)

where ${\sigma }_p$ is the NMAD of the elevation changes on stable terrain, and $N_{\mathrm{eff}}$ is the effective number of independent values on stable terrain, calculated as:

$$N_{\mathrm{eff}}={\displaystyle \frac{N_{\mathrm{tot}}·P_s}{2·l}},$$

(4)

where $N_{\mathrm{tot}}$ is the total number of values, $P_s$ is the pixel size of the DEM differencing image, and l is the spatial autocorrelation distance (500 m). The area uncertainty, ${\sigma }_{A_g}$, is assumed to be $8\%$ of the glacier area [52]. The uncertainty in the density conversion factor, ${\sigma }_{{\rho }_c}$, is $\pm 60$ kg ${\mathrm{m}}^{-3}$[50]. The penetration depth uncertainty, ${\sigma }_{d_p}$, is estimated for our glacierized region (elevation range $4700-6500$ m), accounting for the slope error ($\pm 2.7\times {10}^{-4}$) and intercept error ($\pm 1.301$ m) from the linear relationship.

3.4. Glacier Surface Velocities

We used ITS_LIVE (Inter-Mission Time Series of Land Ice Velocity and Elevation) product to examine the spatio-temporal changes of glacier surface velocity. The ITS_LIVE data product used here is a set of compilations of annual mean surface velocities for HMA. It is derived from Landsat 4, 5, 7, and 8 imagery using feature tracking algorithm through the auto-RIFT pipeline [53]. It spans the time period from 1985 to 2018 with a spatial resolution of 240 m. We extracted the average surface velocities along the glacier centerlines using bilinear interpolation method between 1985 and 2018.

Although ITS_LIVE glacier velocity product has been widely used in glaciological studies [54,55], its accuracy remains unclear, particularly in the Himalayas. Here we used in-situ surface horizontal velocity measurements on the JMYZ Glacier to assess the applicability of ITS_LIVE in our study region. Repeated stake measurements were carried out on October 1, 2010 and May 26, 2011 using differential global positioning system (Figure 1). Due to logistical difficulty and remoteness, only seven stakes were installed in the ablation zone, and the measurements did not cover a complete hydrological year with data gap in summer. Satellite-derived seasonal glacier velocity in the western Himalayas suggested that the peak flow speed occur in spring, which are much larger than the summer speed [56,57]. To derive the annual average glacier velocities from stake measurements, we thus multiplied the measured velocities by a constant factor $0.7$ following [56,57].

3.5. Glacier Thicknesses

We used the ice thickness data as an input for ice flow model to simulate the glacier flow regimes. The ice thickness data utilized in this study was derived from the composite product developed by [58], which provides a consensus estimate of the ice thickness distribution for global glaciers. This dataset was generated through a combination of up to five different models utilizing principles of ice flow dynamics to invert for ice thickness from surface characteristics, such as elevation and slope. The models were assessed for performance against field measured ice thicknesses through a cross-validation scheme, and the final composite solution is obtained by applying inverse variance and bias weighting to the individual model outputs. This ice thickness dataset enhances our understanding of glacier volume and provides critical inputs for numerical ice flow models.

3.6. Ice Flow Model

We used the numerical ice flow model PoLIM to inferred the flow regimes of four typical glaciers in our study region (see Section 2) for two different years, approximately 2000 and 2017, by assimilating the ITS_LIVE velocity dataset. By comparing these two flow regimes, we investigated the impacts of lake evolution and surface elevation changes on glacier dynamics. PoLIM is a Blatter-Pattyn type higher-order flowband model with inclusion of longitudinal stress [59], which is essential for accurately modeling mountain glacier dynamics. Here, we briefly describe this model and our simulation strategies.

3.6.1. Model Description

We first define a Cartesian coordinate system ($x,y,z$), where x denotes the direction along the flowline, y represents the direction transverse to the flowline, and z is the vertical coordinate perpendicular to both x and y. For PoLIM, the simplified form of Stokes equation is described by

$${\displaystyle \frac{\partial }{\partial x}}\left(2{\tau }_{xx}+{\tau }_{yy}\right)+{\displaystyle \frac{\partial {\tau }_{xy}}{\partial y}}+{\displaystyle \frac{\partial {\tau }_{xz}}{\partial z}}={\rho }_{i}g{\displaystyle \frac{\partial s}{\partial x}},$$

(5)

where ${\tau }_{ij}$ is the deviatoric stress tensor, ${\rho }_i$ is the ice density, g is the accelerated gravity, and s is the surface elevation. The lateral drag (${\tau }_{xy}$) is parameterized by the glacier’s half widths (W) as

$${\tau }_{xy}=-{\displaystyle \frac{\eta u}{W}},$$

(6)

where u is the horizontal velocity and $\eta$ is the ice viscosity. For modeling the creep behavior of glacier ice, a Glen-type flow law is used as the constitutive equation:

$${\tau }_{ij}=2\eta {\dot{ϵ}}_{ij},\eta ={\displaystyle \frac{1}{2}}{A}^{-1/n}{\dot{{ϵ}_e}}^{(1-n)/n},$$

(7)

where ${\dot{ϵ}}_{ij}$ is the strain rate tensor, ${\dot{ϵ}}_e$ is the effective strain rate, A is the flow rate factor, and n is the flow law exponent.

At the glacier surface, a stress-free surface boundary condition is assumed. For PoLIM, this boundary condition reads:

$$(2{\tau }_{xx}+{\tau }_{yy}){\displaystyle \frac{\partial s}{\partial x}}-{\tau }_{xz}=0.$$

(8)

At the glacier bedrock, a linear friction law is applied which relates the basal drag ${\tau }_b$ to the sliding velocity $u_b$:

$${\tau }_b=-{\beta }^2{u}_b,$$

(9)

where ${\beta }^2$ is the basal friction parameter and is positive.

3.6.2. Simulation Strategy

Longitudinal glacier topographies, including ice thickness, surface and bed elevations, are necessary inputs for the PoLIM model. We first extracted the surface elevation from the SRTM DEM (referred to as $s_{\mathrm{srtm}}$) and the ice thickness (referred to as $H_{2000}$) along the centerline. The bedrock elevation b was then derived by substracting the ice thickness $H_{2000}$ from the SRTM elevation $s_{\mathrm{srtm}}$. Assuming a fixed bedrock topography, we calculated the ice thickness in 2017 ($H_{2017}$) using the equation $H_{2017}=s_{\mathrm{tandem}}-b$, where $s_{\mathrm{tandem}}$ represents the TanDEM-X DEM elevation. It is worth noting that we consider the year used here to be the calendar year.

We utilized the Robin inversion method implemented in PoLIM to infer glacier flow regimes [60]. This method is a numerical approach developed to solve the inverse problem of determining optimal basal conditions from glacier topography and surface velocity observations [61,62,63]. It iteratively solves the higher-order approximated equations with two types of boundary conditions, i.e., a Neumann condition (see Equation 8) and a Dirichlet condition representing the free surface and measured velocities, respectively. The method uses a cost function to minimize the mismatch between modeled and observed velocities, and iteratively adjusts the basal friction parameter ${\beta }^2$ to achieve this minimization. The Dirichlet condition at the surface is set as

$$u\left(s\right)=u_{\mathrm{itslive}},$$

(10)

where $u\left(s\right)$ is the horizontal surface velocity and $u_{\mathrm{itslive}}$ is the ITS_LIVE velocitity. To reduce the uncertainty inherent in ITS_LIVE, we adopted the mean velocities from the periods 1999-2001 and 2016-2018 as representative values for the years 2000 and 2017, respectively. For simplicity, we negeleted the thermo-mechanical coupling effect and assumed a constant flow rate factor A. Through assimilating ITS_LIVE velocities and adjusting the basal friction parameters, we therefore modeled the horizontal velocity fields using the given glacier topographies.

4. Results

4.1. Geodetic Glacier Mass Balance Changes

Almost all glaciers experienced thinning during 1974–2000, with ice loss rates significantly intensifying in the first two decades of the 21st century, particularly at the glacier tongues (Figure 2a-b). The average regional geodetic mass balances for these periods were $-0.24\pm 0.11$ m w.e. ${\mathrm{a}}^{-1}$ and $-0.57\pm 0.05$ m w.e. ${\mathrm{a}}^{-1}$, respectively, indicating a twofold increase in the ice loss rate. Figure 2c illustrates the region-wide average mass balance profile versus elevation. Notably, the 2000–2020 geodetic mass balance profile was considerably steeper compared to that of 1974–2000. The differences in ice loss rates between the two intervals markedly increased with decreasing elevation, peaking at nearly $2.0\pm 0.05$ m w.e. ${\mathrm{a}}^{-1}$ below an elevation of 5000 m. We observe that the loss rates in the 1974–2000 profile were relatively steady across elevations of $5600–6000$ m, with small fluctuations around the mean value of $-0.16\pm 0.05$ m w.e. ${\mathrm{a}}^{-1}$. Below an elevation of 5600 m, the loss rates increased significantly, while in the 2000–2020 profile, the loss rates sharply decreased below an elevation of 6250 m.

Figure 3 shows the temporal variations of geodetic mass balances along the centerlines of four selected glaciers over the periods 1974–2000, 2000–2010, and 2010–2020. We found that the geodetic mass balances of the four glaciers at lower elevations became significantly more negative since the start of the 21st century. However, the signals of temporal variations in ice thickness changes at the upper areas of individual glaciers were relatively noisy, particularly for the profiles derived from the Hugonnet dataset. These profiles displayed jagged oscillations in geodetic mass balances, which changed rapidly over short spatial distances. For example, the 2000–2010 geodetic mass balance of JMYZ Glacier showed sharp fluctuations at distances of approximately $1.5$ to $4.0$ km from the glacier headwall, reaching a peak of $4.4\pm 0.05$ m w.e. ${\mathrm{a}}^{-1}$ at around $1.6$ km, before dropping to a minimum of $-7.8\pm 0.05$ m w.e. ${\mathrm{a}}^{-1}$ at approximately $3.2$ km. This phenomenon was likely caused by uncertainties in DEM production, resulting from low radiometric contrast in areas at high elevations covered by snow and firn. Therefore, we will focus our analysis on the mid- and low-elevation areas.

The mass balance profiles at the tongues of lake-terminating glaciers (Figure 3c-d) were notably less smooth than those of land-terminating glaciers (Figure 3a-b), indicating the impact of proglacial lakes on ice thickness changes. Generally, the lake-terminating glaciers showed higher mass loss level at the termini than the land-terminating glaciers. However, N241 Glacier also exhibited higher mass loss pattern similar to lake-terminating glaciers, indicating that land-terminating glaciers may also have comparable loss rates with lake-terminating glaciers. The highest mass loss of JMYZ Glacier reached during the period 2000-2010, while ASJG Glacier reached during the period 2010-2020, which indicated that the mass loss of lake-terminating glaciers was likely related to the evolutional stages of proglacial lakes. Although the land-termianting glacier N171 experienced intensified mass loss from the period 1974-2000 to the period 2010-2020, it generally exhibited stably increased mass loss rates towards the terminus.

The mass balance profiles at the tongues of lake-terminating glaciers (Figure 3c-d) were notably less smooth compared to those of land-terminating glaciers (Figure 3a-b), highlighting the influence of proglacial lakes on ice thickness changes. In general, lake-terminating glaciers showed greater ice loss at their termini than land-terminating glaciers. However, N241 Glacier exhibited a similar pattern of accelerated ice thinning as seen in the lake-terminating glaciers (Figure 3c), suggesting that land-terminating glaciers can also experience comparable rates of ice loss. The most significant ice reduction for JMYZ Glacier occurred during the 2000–2010 period, while ASJG Glacier’s greatest ice loss rate was recorded during 2010–2020. This suggests that the mass loss of lake-terminating glaciers may be linked to the evolutionary stages of proglacial lakes. Although the land-terminating N171 Glacier experienced intensified ice loss from 1974–2000 to 2010–2020, its geodetic mass balance decreased steadily toward the terminus (Figure 3d).

4.2. Glacier Velocity Changes

Before we quantify the glacier velocity changes using ITS_LIVE product in our study region, we first evaluate its performance in comparison to measured surface velocities on JMYZ Glacier. As the time span of stake measurements overlaps the acquisition dates of annual ITS_LIVE velocities in 2010, 2011, and 2012, we extract the ITS_LIVE velocities at the stakes and along the centerline in these three years using bilinear interpolation method. Figure 4 compares the mean annual ITS_LIVE velocities between 2010 and 2012 with our measurements. It shows that the ITS_LIVE velocities are in good agreement with the measurements except at the stakes A and B (Figure 4a). The mean absolute deviation between ITS_LIVE and measured velocities is 7.2 m ${\mathrm{a}}^{-1}$ corresponding to 25% of the largest measured velocities. However, the mean absolute deviation significantly decreases to 2.7 m ${\mathrm{a}}^{-1}$ (corresponding to 10% of the largest measured velocities) after excluding the data at the stakes A and B. The mismatch at these two stakes is possibly due to coarse resolution of ITS_LIVE or local failure of feature tracking. The stakes were roughly along the glacier centerline (Figure 1). Figure 4b shows that the ITS_LIVE velocities generally capture the spatial pattern of surface flow speed changes on JMYZ Glacier. Through our comparisons, we think that ITS_LIVE product is relatively reliable and can be applied to our study.

We investigated the spatio-temporal glacier speed changes in our study region using time series of ITS_LIVE velocities for the period of 1990-2018. Figure 5 shows the multi-year averaged velocity maps for the periods 1990-1999, 2000-2009, and 2010-2018. Here we used the 2010-2018 velocity map to demonstrate the regional flow speed characteristics (Figure 5c). We can see that large glaciers generally flowed faster than small glaciers. The highest flow speed was more than 40 m ${\mathrm{a}}^{-1}$ occurring at the JMYZ Glacier. The glacier-wide mean velocity across this region ranged from 0.4 m ${\mathrm{a}}^{-1}$ to 11.1 m ${\mathrm{a}}^{-1}$. It also shows that the surface velocities of lake-terminating glaciers were substantially higher than those of land-terminating glaciers. Lake-terminating glaciers occupy a relatively large surface area, and their surface velocities were substantially higher than those of land-terminating glaciers.

To characterize the regional velocity trends, we examined the temporal velocity changes for the four selected glaciers along the centerlines (Figure 6). We can see that the spatial patterns of lake-terminating glacier velocity profiles were more complex than those of land-terminating glaciers, which generally showed a single maximum velocity pattern. Although the selected glaciers were exposed to similar environment settings, their temporal velocity trends were significantly different. The surface velocities of JMYZ Glacier were at a relatively low level less than 23 m ${\mathrm{a}}^{-1}$ in 1990s (Figure 6a). In 2000s, however, the surface velocities substantially increased, particularly at the glacier terminus (6.0-11.0 km). The maximum velocity increased by 54% and reached 34.7 m ${\mathrm{a}}^{-1}$. During the period 2010-2018, the terminus velocity decreased slightly, but the velocity in the upper profile (2.2-7.7 km) showed a pronounced increase. For the ASJG Glacier, its velocity trends during the investigated period were more complicated than JMYZ Glacier (Figure 6b). It exhibited a high flow regime in 1990s with a mean value of 17.2 m ${\mathrm{a}}^{-1}$ and a maximum value of 38.4 m ${\mathrm{a}}^{-1}$. In 2000s, we observed a significant velocity slowdown along most sections of the profile. Velocity increase only occurred in the upper section (2.1-3.1 km) and at the glacier tongue (10.5-11.5 km). The surface velocities during the period 2010-2018 showed further slowdown, particularly at the glacier front. For the land-terminating glaciers, i.e., N241 Glacier and N171 Glacier, they exhibited relatively simple patterns of longitudinal velocities compared to lake-terminating glaciers (Figure 6c-d). N241 Glacier showed a significant speedup since the 2000s (Figure 6c). The maximum velocity increased by about 60% from 12.6 m ${\mathrm{a}}^{-1}$ in 1990s to 20.1 m ${\mathrm{a}}^{-1}$ for the period 2010-2018. However, the terminus velocities showed a gradual slowdown over the investigated period. For the N171 Glacier, its mean surface velocity was approximately two times higher than N241 Glacier. It also showed speedup since 2000s particularly at the central section. During the period 2010-2018, a slightly speed increase was observed in the upper section (0.5-3.7 km), but the peak velocity decreased. The terminus velocity of N171 Glacier also exhibited gradual reductions since 1990s.

4.3. Glacier Terminus Changes

Figure 7 shows the terminus positions of the four glaciers between 1990 and 2019 with a temporal resolution of approximately five years. It is evident that all glaciers have significantly retreated over the past three decades. The recession of lake-terminating glaciers (JMYZ and ASJG) was larger in both magnitude and rate compared to land-terminating glaciers (N241 and N171). Specifically, the lake-terminating glaciers retreated more than 1000 m during the study period 1990-2019, while the land-terminating glaciers retreated less than 750 m. Among them, the length change of land-terminating glacier N171 was the smallest, approximately $-516$ m over the study period. In contrast, the ASJG Glacier had the most significant length change, approximately $-1600$ m, which was three times the length change of the N171 Glacier. The ASJG Glacier also exhibited the fastest retreat rate compared to the other three glaciers in all epochs between 1990 and 2019. The retreat rates of the glacier termini showed no clear overall trend during the study period. However, the maximum retreat rates for all glaciers occurred between 2005 and 2010, with the ASJG Glacier receding at a rate of more than 70 m ${\mathrm{a}}^{-1}$. During the period 2015-2019, the lengths of all glaciers significantly reduced, and the retreat rate of the land-terminating glacier N241 ($40.5$ m ${\mathrm{a}}^{-1}$) was comparable to that of the lake-terminating JMYZ Glacier ($39.4$ m ${\mathrm{a}}^{-1}$).

Figure 7. Temporal changes in glacier termini from 1990 to 2019 for four glaciers. (a) JMYZ Glacier, (b) ASJG Glacier, (c) N241 Glacier, (d) N171 Glacier. The background in each panel is a Landsat 8 image.

Figure 7. Temporal changes in glacier termini from 1990 to 2019 for four glaciers. (a) JMYZ Glacier, (b) ASJG Glacier, (c) N241 Glacier, (d) N171 Glacier. The background in each panel is a Landsat 8 image.

4.4. Modelled Glacier Flow Regimes

Figure 8 and Figure 9 show the modelled horizontal velocity fields in 2000 and 2018 along the longitudinal profiles for the lake- and land-terminating glaciers, respectively. Unlike the surface velocity profiles shown in Figure 6, the velocity fields are capable of capturing changes in interior and basal flow speeds responding to variations in glacier geometry, bed condition, and proglacial lake. To help understanding the flow regime changes, the modelled surface and basal velocities through robin inversion were also presented (Figure 10).

We observed that the selected glaciers display distinct flow regimes regardless of their terminus types. For the lake-terminating glaciers, the changes in their flow patterns showed significant differences. The flow regime of JMYZ displayed a clear speedup between 2000 and 2018 (Figure 8a and c). The high flow speed zone migrated upstream from the terminus towards the headwall during the study period. The increased flow speeds occurred not only on the glacier surface as shown in Figure 10a but also within its interior (Figure 8c). For the other lake-terminating glacier ASJG, its flow velocity field generally decreased with substantial surface thinning (Figure 8b and d). In 2000, the highest velocity zone was located at the glacier terminus. However, by 2018, the terminus speeds significantly decreased, and the high flow zone shifted to an up-glacier section at approximately $\sim 2.0-4.0$ km along the longitudinal profile.

For the land-terminating glaciers, i.e., N241 Glacier and N171 Glacier, their flow regime changes can be characterized as accelerating and stable. The velocity fields of N241 showed a substantial transition from a slower to a faster flow pattern during 2000-2018 (Figure 9b and d). In 2000, the high flow zone was located in the central section $\sim 1.5-2.5$ km of the longitudinal profile. However, by 2018, the high velocity zone had significantly expanded covering the section $\sim 1.0-4.0$ km along the longitudinal profile, and the basal layer also exhibited fast flow speeds. The N171 Glacier showed a relatively stable flow structure even though the glacier surface largely thinned during this period (Figure 9a and c).

5. Discussion

5.1. Reasons for Heterogeneous Glacier Flow Regime Changes

Our numerical modeling revealed highly heterogeneous patterns of glacier flow regimes within our study region, indicating that glacier dynamic changes were controlled by different factors and complex feedbacks. To elucidate the mechanisms underlying the evolution of flow patterns, we performed a force balance analysis along the longitudinal profile (Figure 11 and Figure 12). For stress equilibrium, the driving stress (${\tau }_d$) is balanced by resisting forces, which include basal shear stress (${\tau }_b$), longitudinal drag (${\tau }_l$), and side drag (${\tau }_w$) [64]. It should be noted that the longitudinal stress can be either positive or negative.

The spatial patterns of driving stress (${\tau }_d$) on the JMYZ Glacier remained consistent from 2000 to 2018, with only a slight decrease in ${\tau }_d$ observed at the terminus. The slowdown in terminus velocity was due to decreased driving stress and increased basal shear stress, while the accelerated flow in the upper section of the profile can be attributed to a significant reduction in basal shear stress (Figure 8a, Figure 11a and Figure 12a). Longitudinal coupling may have also contributed to the upstream acceleration of ice flow by transmitting the high terminus velocity observed in 2000. Despite the distance between the terminus and the upper fast-flow region being ten times greater than the average ice thickness of the JMYZ Glacier (approximately 170 m), longitudinal coupling could still influence flow variations upstream of the terminus due to ineffective basal resistance and side drag (Figures S1a and S2a) [64,65].

The driving stress of the ASJG Glacier rapidly decreased between 2000 and 2018, particularly in the lower ($7.0-11.0$ km) and upper ($0.0-2.0$ km) sections (Figure 11b), primarily due to reduced ice thickness. This significant reduction in driving stress contributed to the deceleration of terminal velocity (Figure 8b and Figure 10b). The increased flow in the upper section ($2.0-4.0$ km) was driven by a decrease in basal resistance, which dropped by more than $50\%$ during this period (Figure 12b). While ice-marginal lakes are generally believed to exacerbate mass loss at the terminus, steepen the ice surface, and increase driving stress [28], our findings suggest that sustained mass loss leading to significant thinning of the ice can largely reduce driving stress and slow down flow speeds.

The fast-flow region of the land-terminating N241 Glacier expanded both upstream and downstream between 2000 and 2018, with significant intensification during this period (Figure 9a and c). Force balance analysis revealed that increased driving stress between approximately $\sim 1.65-1.95$ km and a substantial decrease in basal shear stress between $\sim 1.9-5.0$ km were responsible for this expansion and acceleration (Figure 11c and Figure 12c). Additionally, the surface slope of the N241 Glacier steepened around the $2.0$ km from 2000 to 2018 (Figure 9a and c), further increasing driving stress. Longitudinal coupling played a crucial role in transmitting these increased speeds both upstream and downstream. An extended flow zone became evident in the upper section of the N241 Glacier (Figure S1c). In contrast, the spatio-temporal patterns of the velocity field and force balance terms for the N171 Glacier remained consistent over the study period (Figure 9d, Figure 10d, Figure 11d and Figure 12d and Figures S1d-S2d), indicating relatively stable flow dynamics for this land-terminating glacier. However, there was a slight increase in surface velocities upstream and a decrease downstream around 4.0 km (Figure 10d), which can be attributed to increased basal sliding (Figure 12d).

Our observations and simulations of heterogeneous changes in ice flow regimes within a specific region appear to contrast with widespread reports of glacier slowdown in High Mountain Asia (HMA) [66]. The dynamic responses of individual glaciers to climate change are complex and influenced by various factors. However, large-scale satellite studies often emphasize glacier-wide velocity changes without adequately accounting for detailed spatio-temporal variations. This study suggests that the flow dynamics of individual glaciers can be significantly influenced by increased basal sliding due to surface meltwater input and steepened surface slopes. To improve understanding, future large-scale satellite monitoring of glacier velocity changes should place greater emphasis on capturing the variations among individual glaciers across both space and time.

5.2. Implications for glacier evolution models

Our study reveals that changes in glacier surface velocity are influenced by several factors: surface topography, ice thickness, basal lubrication, and longitudinal coupling. These factors affect variations in force balance, thereby impacting the ice flow regime. Current large-scale glacier evolution models can predict future changes in individual glaciers [37,67,68,69], but they often oversimplify ice dynamics. Some advanced models employ the shallow ice approximation to compute ice flux [67,68], yet they fail to account for longitudinal stress-gradient coupling, and none consider the impacts of subglacial hydrology on ice dynamics in a changing climate.

In HMA, glacial lakes are projected to increase in number, area, and volume by the end of the century [34]. This expansion will undoubtedly influence glacier mass balance and dynamics. However, current glacier evolution models generally overlook frontal ablation in water-terminating glaciers [70], instead treating glaciers as land-terminating [71]. To improve frontal ablation modeling, glacier evolution models should incorporate ice-lake interactions. Given limited data on marginal lake depth and temperature, developing a parameterization for ice-lake interaction is crucial to reduce computational costs.

Our study also emphasizes that changes in surface topography, particularly in surface slope, can significantly affect driving stress and ice velocity. Accurately modeling changes in surface topography requires careful consideration of both mass balance and ice dynamics. Many models currently compute elevation-banded surface mass balance [69,70], but this approach introduces uncertainties in spatial surface elevation change. Future models should integrate longitudinal coupling, frontal calving dynamics, and distributed surface mass balance to better capture interactions between ice and marginal lakes, thereby improving predictions of ice flow changes. It is essential for next-generation glacier evolution models to incorporate higher-order or full-Stokes dynamics to accurately simulate the complex physical processes responding to climate change [72,73].

5.3. Limitations

In this study, we used multi-sourced datasets, including field measurements, DEMs, satellite optical images, glacier velocity products, and glacier elevation change data, to analyze the changes in land- and lake-terminating glaciers at the headwaters of the YZR. These datasets, however, are not temporally consistent, which created challenges in comparing changes over time in mass balance, terminus position, and velocity. For example, DEMs generated from topographic maps can date back to the 1970s, while other datasets are more reliable only after the 1990s due to the increasing availability of satellite images. Additionally, some of the datasets are open-access global products typically available at relatively coarse spatial resolutions. While these datasets are useful for large-scale analysis, they may lack the spatial detail necessary to capture localized glacier dynamics. The ITS_LIVE velocity product, for instance, provides long-term, high-temporal-resolution data across broad areas, but its 240 m resolution limits the detection of fine-scale velocity variations. To address this limitation, we focused on four relatively large glaciers for detailed analysis.

To characterize and compare the mass and flow changes between land- and lake-terminating glaciers, we selected two of each type from our study region. Ideally, for a meaningful comparison, the selected glaciers would be similar in size (e.g., area, length). However, both lake-terminating glaciers in this region are larger than the land-terminating ones in both area and length. While this selection highlights the dynamic differences between glacier types, some of these differences, particularly in velocity changes, could be influenced by the disparity in glacier size. It is challenging to find glaciers of similar size within the same region that represent both types. Nonetheless, our comprehensive analysis was able to reveal significant differences between the two glacier types by examining the spatial patterns and temporal changes in surface velocities and flow regimes.

We inferred glacier flow regimes by assuming steady-state geometry and incorporating ITS_LIVE velocity data, using an inversion method to estimate the initial dynamic state of ice flow. However, glacier flow in any given year may not be in equilibrium with its present geometry, as transient effects from past dynamics can still influence current behavior. Our modeling treated the glaciers as isothermal, with a constant flow rate factor, A. In reality, the glaciers in this region are likely polythermal, as ground-penetrating radar detected temperate ice at the basal layer of JMYZ Glacier [74]. Accurately modeling glacier flow would require thermomechanical modeling that considers transient thermal conditions and surface elevation changes. Despite these limitations, the goal of this study was to capture long-term changes in glacier flow regimes, rather than replicating the precise dynamics of each glacier.

6. Conclusions

This study presents the changes in geometry and flow speed of both land- and lake-terminating glaciers at the YZR headwaters over the past four decades, using a combination of field measurements, remote-sensing observations, and numerical ice flow modeling. We found that glacier behavior in this region exhibited significant heterogeneity, not only between different glacier terminus types but also among glaciers of the same type, even those located close to each other.

From 1974 to 2020, we observed widespread, accelerated glacier thinning and mass loss in our study region. The region-wide average ice mass loss rate more than doubled in the first two decades of the 21st century ($-0.57\pm 0.05$ m w.e. ${\mathrm{a}}^{-1}$) compared to the period 1974–2000 ($-0.24\pm 0.11$ m w.e. ${\mathrm{a}}^{-1}$). Although lake-terminating glaciers lost more mass than land-terminating glaciers, we found that the land-terminating N241 Glacier exhibited a comparable mass loss rate to lake-terminating glaciers. Between 1990 and 2019, lake-terminating glaciers showed significantly greater and faster retreat, exceeding 1000 m, whereas land-terminating glaciers retreated by less than 750 m. The highest retreat rates occurred between 2005 and 2010, with further substantial length reductions observed between 2015 and 2019.

The ITS_LIVE velocity dataset proved reliable in capturing glacier surface velocities, aligning well with in situ measurements. Our analysis of glacier velocity changes from 1990 to 2018 revealed that larger, lake-terminating glaciers consistently exhibited higher flow speeds than smaller, land-terminating glaciers. Additionally, lake-terminating glaciers demonstrated more complex and variable velocity trends over time, while land-terminating glaciers generally showed more straightforward acceleration. These results emphasize the distinct dynamic behaviors of different glacier types.

Numerical modeling revealed distinct flow regime changes in both lake- and land-terminating glaciers from 2000 to 2018. Lake-terminating glaciers displayed contrasting behaviors: JMYZ Glacier experienced significant speedup, with the high flow speed zone expanding from the terminus to the headwall, while ASJG Glacier showed a decrease in flow velocity, with the high flow zone shifting up-glacier. Among the land-terminating glaciers, N241 Glacier showed substantial acceleration, expanding the high flow zone along its length, while N171 Glacier maintained a stable flow regime despite notable surface thinning.

Our study underscores the significant and likely increasing impact of lake-terminating glaciers on ice mass loss and downstream meltwater supply. Accurate modeling of these glaciers is crucial due to their role in providing meltwater to Himalayan rivers. Our findings emphasize the need for improved glacier models that account for changes in surface topography, including variations in slope, as well as longitudinal coupling and frontal calving dynamics. Incorporating higher-order or full-Stokes dynamics is essential to capture the complex processes driving glacier behavior in a changing climate. This study provides new insights into glacier dynamics associated with proglacial lake formation, advancing our understanding and highlighting the need for future research to refine predictions of glacier changes and their hydrological impacts.