Topography and accumulation rate as controls of asynchronous surging behaviour in the eastern and western branches of the Western Kunlun Glacier, Northwestern Tibetan Plateau

ABSTRACT

The western Kunlun main peak region is among the areas in High Mountain Asia where surge-type glaciers are highly concentrated. Here, we analyse the surging characteristics of the eastern and western branches of the Western Kunlun Glacier and the factors controlling the asynchronous behaviour of their surges. The eastern branch entered an unstable state in 1999 and remained so until the culmination of its surge in the summer of 2019, spanning 21 years. Conversely, the surge of the western branch commenced in the summer of 2020. The surge duration for this glacier was four years, characterized by a rapid acceleration and deceleration process. Based on the glacier surge characteristics, we posit that western branch of Western Kunlun Glacier was influenced by hydrological mechanisms, while eastern branch was affected by subglacial thermal processes. These process intensifies crevasse formation on the glacial surface, providing conduits for surface meltwater to reach the glacier bed, thus elevating subglacial water pressure. The difference of subglacial hydrology and thermal processes caused by different subglacial topography and mass accumulation rates was the main factor of the asynchronous behaviour of the west and east branches of the West Kunlun glacier.

KEYWORDS:

• Mechanism of glacier surge
• glacier surface elevation
• glacier velocity
• Western Kunlun
• remote sensing

1. Introduction

Glacier surges, characterized by periodic and rapid movement of glaciers over relatively short periods, are a phenomenon observed in certain glaciers referred to as ‘surge-type glaciers (STGs)’ (Meier and Post 1969). Recent statistical analyses of literature have revealed a global total of approximately 1850 STGs, which exhibit a distinct clustered distribution pattern (Guo et al. 2022). High Mountain Asia (HMA) is among the regions worldwide where STGs are concentrated (Guillet et al. 2022). In recent years, the frequency of glacier surge-related hazards, including collapses, and glacial lake outburst floods, has been on the rise in this region, making it a focal point in glaciological research (Gao et al. 2021). Over the past decade, numerous scholars have utilized multisource remote sensing data to investigate the spatial distribution of STGs at regional or overall scales in the HMA based on various criteria, such as glacier terminus advance, abnormal surface flow velocity, and significant changes in glacier surface elevation, and corresponding inventories of STGs have been published (Bhambri et al. 2017; Goerlich, Bolch, and Paul 2020; Guillet et al. 2022; Guo et al. 2023; Yao et al. 2023). These efforts have provided invaluable datasets that contribute to our deeper understanding of the spatiotemporal distribution patterns of STGs in the HMA and enhance our disaster mitigation capabilities.

Nonetheless, our understanding of glacier surges remains limited, primarily due to the abrupt and complex nature of surge events, as well as the remoteness of the regions where STGs are typically found. For example, why do neighbouring glaciers, subject to similar climatic conditions, surge asynchronously? Why do some glaciers exhibit prolonged periods of acceleration during surging events, while others experience shorter durations? Addressing these questions is challenging, primarily because obtaining field observations of crucial glacier internal and subglacial parameters is extremely difficult (Guo et al. 2022). Despite these challenges, many glaciologists have proposed various mechanisms to elucidate the diversity in glacier surge behaviour. Currently, the most widely applied mechanisms are hydraulic and thermal theories, originally formulated based on field observations during the surges of Variegated and Trapridge glaciers in the 1980s (Clarke 1976; Kamb et al. 1985; Nolan et al. 2021). However, since the beginning of the twenty-first century, studies have indicated that glaciers exhibiting these different surging mechanisms could coexist in the same regions, such as Greenland and the Karakoram (Quincey et al. 2015). In recent years, some researchers have introduced conceptual glacier surge models based on limited observational data and assumptions, with the aim to explain the relationships among glacier surges, glacier geometry, and climatic conditions and reveal the evolution of glacier surges (Sevestre and Benn 2015). Examples include the structural glaciology-based surge glacier evolution model proposed by Clarke and Hambrey (2019) and the conceptual models grounded in enthalpy theory advanced by Benn et al. (2019). However, the complexity and numerous parameters involved in these models have made their widespread application challenging.

Currently, advances in remote sensing technology and the accumulation of remote sensing data have made it possible to comprehensively study the details of glacier surges (Guillet et al. 2022). This, in turn, allows for a deeper understanding of the control mechanisms governing surging glacier behaviour in remote regions (Guo et al. 2022; Lovell, Carr, and Stokes 2018). The West Kunlun Glacier (WKG), situated in the northwest of the main peak of the West Kunlun Mountains, is a typical cold-type glacier. However, we have observed strikingly different glacier surges between its eastern and western branches, indicating significant disparities. In this manuscript, we utilize Digital Elevation Models (DEMs) and glacier surface velocity data derived from multiple remote sensing sources to elucidate the mechanisms and predominant factors governing the contrasting surges of the two glacier branches. Our findings are poised to significantly enhance our understanding of surging glacier characteristics and their underlying mechanisms in the Western Kunlun Main Peak region.

2. Study sites

The Western Kunlun Main Peak (WKMP, 80°E∼82°E, 35°N∼36°N) represents one of the largest contemporary glacier regions on the Tibetan Plateau (Figure 1a) (Yasuda and Furuya 2015). Its southern slope features gently sloping terrain and is predominantly characterized by wide-tailed valley glaciers, with no glaciers below 5200 m in elevation. On the northern slope, the glaciers are relatively longer and steeper, comprising mainly compound valley glaciers and dendritic valley glaciers, with some glaciers extending below 5000 m in elevation (Figure 1b) (Guan et al. 2021). According to recent research, in 2020, the WKMP harboured a total of 440 glaciers covering an area of 2964.59 ± 54.87 km² (Zhang et al. 2023). Guan et al. (2021) conducted a comprehensive reassessment of the distribution of STGs in the region from 1972 to 2020 using multisource remote sensing data. Their findings indicate that among the 440 glaciers, ten are confirmed to be STGs, three are likely STGs, and five are potential STGs (Guo et al. 2023). Located in the northwestern part of the peak area, the WKG consists of two main branches (Figure 1c). According to the Second Chinese Glacier Inventory, the eastern branch covers an area of 51.48 km² with an average slope of 9.2° (Guo et al. 2015). In contrast, the western branch glacier has an area of 76.92 km² and exhibits a relatively gentle slope (Guo et al. 2015). The surge in EWKG began before 2010, while detailed research on the surge in WWKG is currently lacking. Additionally, this region is primarily influenced by the mid-latitude westerlies, characterized by a cold and semiarid climate. The mean annual temperature around the snow line (5930 m) is −13.9°C, with an annual precipitation of 300 mm. Precipitation is predominantly in the form of solid precipitation and is concentrated mainly from May to August, representing a summer-recharge type glacier system (Zhang et al. 2023).

Figure 1. (a) Location of the WKG. Background: Shuttle Radar Topography Mission (SRTM DEM and its hillshade) and glacier outlines derived from the Randolph Glacier Inventory V6.0. (b) Glacier distribution on the WKMP (c) The background of the main image is a false colour composite of Landsat 9 Operational Land Imager (OLI) (Band 6, 5 and 3, LC09_L1TP_145035_20230915_20230915_02_T1). The curved scale bar up the WKG indicates the longitudinal profile used for surface velocity and elevation analysis.

3. Data and methods

3.1. Data

We obtained one Landsat 2 Multi-Spectral Scanner (MSS) image on July 12, 1977, 20 Landsat 5 Thematic Mapper (TM) images from 1989 to 2011, three Landsat 7 Enhanced Thematic Mapper Plus (ETM+) images from 2000 to 2012, and 118 Landsat 8/9 Operational Land Imager (OLI) images from 2013 to 2023, with minimal cloud and snow cover. Furthermore, to characterize the evolution of crevasses during the glacier active phase, we acquired eight scenes of Sentinel-2 imagery. These images were sourced from the United States Geological Survey (USGS, http://glovis.usgs.gov/) and the European Space Agency (ESA, https://sentinel.esa.int). It was worth noting that this study exclusively utilized Level-1 products of Landsat series and Sentinel-2 optical imagery. These datasets have undergone geometric and radiometric correction by USGS and ESA, and no further processing was conducted in this study. On the other hand, we procured NASADEM, AW3D30 DEM, and COP30 DEM from multiple sources, including the Land Processes Distributed Active Archive Center of the National Aeronautics and Space Administration (NASA, https://cmr.earthdata.nasa.gov), the Earth Observation Research Center of the Japan Aerospace Exploration Agency (JAXA, https://www.eorc.jaxa.jp), and the Copernicus programme of the ESA (https://spacedata.copernicus.eu). These datasets were utilized to assess the changes in glacial surface elevations. The ASTER-L1A data are sourced from NASA Earth Observation Data and were employed in the production of the latest generation of DEMs. Additional datasets, such as the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) and albedo, were obtained from the University of California, Santa Barbara (https://chc.ucsb.edu/data/chirps) and NASA, respectively. This study uses the Inter-Mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) data (Gardner, Fahnstock, and Scambos 2019) obtained from the U.S. National Snow and Ice Data Center (NSIDC, https://nsidc.org) to evaluate the interannual flow velocity changes in the WKG. We utilized the glacier velocity dataset derived from Sentinel-1 radar imagery as outlined by Friedl, Seehaus, and Braun (2021). This dataset offers monthly mean glacier velocity information from October 2014 to July 2021 Table 1.

Table 1. Overview of satellite images and data sources.

Acquisition Date	Sources	Pixel size (m)	Amount	Application
2000	NASADEM	30	1	Estimation of glacier elevation change and mass balance
2006–2011	AW3D30 DEM	30	1
2010–2015	COP30 DEM	30	1
2023	ASTER-L1A	15	2
1990–2018	ITS_LIVE	240	–	Annual glacier surface flow velocity
2014–2021	Sentinel-1 glacier velocity data	200	82	Monthly glacier velocity
1977–2023	Landsat MSS/TM/ETM+/OLI	60/30/15	158	Identification of glacier outlines Glacier surface flow velocity
2016–2023	Sentinel-2	10	8	Assessment of glacier surface crevasse
1981–2022	CHIRPS	5566	–	Estimation the air precipitation in the study area
2000–2022	MCD43A3	500	–	Estimation the albedo in the study area

3.2. Methods

3.2.1. Glacier delineation

The glaciers in the WKMP lack debris cover, which is advantageous for glacier outline identification (Zhang et al. 2023). In this study, we employed ArcGIS 10.8 to manually digitize the boundaries of the WKG. Research has indicated that manual digitization of glacier outlines may lead to inconsistent and nonreproducible results (Hall et al. 2003). Therefore, we adopted the approach recommended by Guo et al. (2015) to assess glacier area errors. The formula is as follows: (1) $\epsilon =N\cdot A$(1) where ϵ was the area error (km²); N was the perimeter of the glacier boundary; and A was the length of half a pixel.

3.2.2. Determination of glacial motion

To obtain high temporal resolution velocity data for the WKG's east and west branches from 2013 to 2023, this study employed the COSI-Corr (Coregistration of Optical Sensing Images and Correlation) software package for coregistration and correlation of Landsat OLI Band 8 images, a tool extensively described by Leprince et al. (2007). Simultaneously, we supplemented ITS_LIVE data with Landsat OLI Band 8 data to extend the record up to 2023. To ensure consistency with the ITS_LIVE data, we adopted the mosaic approach used in the ITS_LIVE data to derive annual surface velocity. This method computes the optimal estimate of glacier surface velocity by performing a weighted average of all data available for the same year (Gardner, Fahnstock, and Scambos 2019). The specific steps are outlined as follows. Initially, we performed coregistration of the acquired Landsat OLI images to ensure consistent spatial referencing with a deviation of less than one pixel. Subsequently, the frequency-domain cross-correlation module within the tool was applied to enhance result accuracy and mitigate image contrast issues. To obtain daily scale velocity data, we employed a methodology involving the utilization of an initial 128 × 128-pixel window, a final 64 × 64-pixel window, and a 2-pixel step size. Regarding annual-scale velocity data, we opted for the same cross-correlation parameters employed in the production of ITS_LIVE data. In addition, to eliminate noise introduced by errors during the coregistration process, all pixel values with a signal-to-noise ratio below 0.9 were systematically discarded. In the end, the outcomes obtained through COSI-Corr represent displacements in the north‒south (NSD) and east‒west (EWD) directions for the image pairs. The total displacement (D) can be calculated using Equation (2(2) $D(x,y)=\sqrt{\mathrm{NSD}{(x,y)}^2+\mathrm{EWD}{(x,y)}^2}$(2) ). (2) $D(x,y)=\sqrt{\mathrm{NSD}{(x,y)}^2+\mathrm{EWD}{(x,y)}^2}$(2) Hence, in combination with the acquisition time interval (Day_i) of the image pairs, the glacier surface daily scale velocity (GVD) and annual scale velocity (GVY) can be obtained, as expressed by the following formula: (3) $\mathrm{GVD}(x,y)=D(x,y)/\mathrm{Da}{y}_i$(3) (4) $\mathrm{GVY}(x,y)=D(x,y)\times 365/\mathrm{Da}{y}_i$(4) Here, GVD (x, y) and GVY (x, y) represent the velocity values in metres per day (m/d) and metres per year (m/yr), respectively, for each pixel. D (x, y) denotes the displacement magnitude for each pixel. A year in this study is considered 365 days. Furthermore, following established practices in previous studies, a 3 × 3 Gaussian low-pass filter kernel is applied to the final velocity results obtained from the glacier region speed images. This filtering step aims to reduce the impact of decorrelation on the results. Notably, in this study, all the processes, including the conversion of displacements into velocities and subsequent postprocessing steps, were implemented using Python.

3.2.3. DEMs and changes in glacier surface elevation, volume and mass balance

This study utilizes DEMs to assess surface elevation changes, volume variations, and mass balance in the east and west branches of the WKG. The DEMs employed include the KH9 DEM, NASADEM, AW3D30 DEM, COP30 DEM, and ASTER DEM. Notably, the differential results between the KH9 DEM and NASADEM were provided by Guo et al. (2023). The ASTER DEM was generated based on ASTER-L1A data from March 30, 2023, using the MicMac ASTER (MMASTER) algorithm proposed by Girod et al. (2017).

Based on the above information, we first resampled all DEMs to a uniform 30 m resolution. Subsequently, we employed the method proposed by Nuth and Kääb (2011) to scrutinize geometric displacements, systematic biases, and elevation-dependent deviations among the DEMs. Based on statistical principles, this method corrects elevation differences in glacier areas by relating the elevation difference to the slope and aspect in stable terrain regions. The equation was as follows: (5) $\mathrm{dh}=a\cdot \cos (b-\phi )\cdot \tan (\alpha )+\overline{\mathrm{dh}}$(5) where dh was the individual elevation difference, a is the magnitude of the horizontal shift, b is the direction of the shift vector, α is the terrain slope, φ was the terrain aspect and $\overline{\mathrm{dh}}$ was the overall elevation bias between the two elevation datasets. After eliminating the errors caused by translation, normalizing the vertical deviation dh with the tangent of the slope at that pixel was achieved, establishing a cosine relationship with the slope and aspect: (6) $\mathrm{dh}/\tan (\alpha )=a\cdot \cos (b-\phi )+\overline{\mathrm{dh}}/\tan (\alpha )$(6) The coefficients a, b, and c were obtained through the least squares fitting, while the translation offsets of various DEM data in the east–west direction (x), north–south direction (y), and vertical direction (z) can be expressed as follows: (7) $x=a\cdot \sin (b),y=a\cdot \cos (b),z=c\cdot \tan (\beta )$(7) where β was the slope. We implemented all the recommended steps of this method in Python.

Following the correction of DEM pairs, we performed differencing to derive changes in glacier surface elevation over the respective time intervals. Consequently, in conjunction with the spatial resolution of the DEMs, we computed the total mass transport during the glacier surge period, expressed by the following formula: (8) $\mathrm{\Delta }V=S_p{\sum }_{i=1}^N\mathrm{\Delta }{h}_i$(8) In the equation, ΔV represents the volume change, S_p denotes the area of a single pixel, N stands for the total number of pixels within the glacier region, and Δh_i signifies the surface elevation change for the i^th pixel.

Converting glacier volume changes into their water equivalent variations yields the glacier mass balance, representing the surplus or deficit between accumulation and ablation over the given period: (9) $B_N={\displaystyle \frac{\mathrm{\Delta }V}{\overline{S}}\cdot {\displaystyle \frac{\rho }{{\rho }_{\mathrm{water}}}={\displaystyle \frac{\mathrm{\Delta }V}{1/2(S_{t0}+S_{t1})}\cdot {\displaystyle \frac{\rho }{{\rho }_{\mathrm{water}}}}}}}$(9) In the equation, B_N represents glacier mass balance in metres water equivalent (m w.e.), ρ denotes the volume-mass conversion factor, which signifies the average density of ice/snow, ρ_water is the density of water, and $\overline{S}$ represents the average glacier area during the period from t₀ to t₁ (Zemp et al. 2013). In this study, we adopt the recommended value of 850 ± 60 kg m⁻³ for the volume-mass conversion factor, with an associated error of 60 kg m⁻³ in ice/snow density, as suggested by Huss (2013). In addition to considering elevation errors, the uncertainty of glacier mass balance must also account for uncertainties in glacier area and the conversion of glacier volume to glacier mass, which are influenced by uncertainties in ice and snow density. The representative formula is given as follows: (10) ${\delta }_B=\sqrt{{\left({\displaystyle \frac{{\rho }_{\mathrm{ice}}}{{\rho }_w\times A}\times {\delta }_V}\right)}^2+{\left({\displaystyle \frac{V}{{\rho }_w\times A}\times {\delta }_{\mathrm{ice}}}\right)}^2}$(10)

In the equation, δ_B represents the uncertainty in mass change, ρ_ice stands for the mean ice density, δ_ice represents the uncertainty in ice density, ρ_w is the density of water, and V and δ_V represent glacier volume change and its uncertainty, respectively.

4. Results

4.1. Variations in glacier surface flow velocities at different time intervals

Figure 2a₁ and b₁ depict the interannual glacier surface velocity variations for WWKG and EWKG, respectively. One of the most notable features is the asynchronous surging behaviour between WWKG and EWKG under similar climatic conditions. In contrast, WWKG experienced a delayed and more abrupt surge onset, while EWKG exhibited a considerably longer acceleration phase. EWKG remained active in its upper and middle parts from 1999 to 2010 (Figure 2b₁). Subsequently, the glacier entered an active phase, with peak velocities reaching 420 ± 1.34 m/yr in 2012 (Figure 2b₁). As the surging front propagated downstream, the peak velocity reached 452 ± 2.31 m/yr in 2012, representing a 14-fold increase from the glacier's corresponding area in 2010 (Figure 2b₁). Thereafter, the glacier gradually entered a deceleration phase. Examination of the velocity variations at different sections of the EWKG during the establishment phase of the surge revealed average speeds of less than 60 m/yr (Figure 2b₃). During the active phase, some sections of the glacier exceeded an average velocity of 150 m/yr. The maximum magnitude of glacier velocity surpassed 52 times that of the quiescent stage (Figure 2b₃). Due to data limitations prior to 2013, we were unable to characterize the intrayear processes of the EWKG surge (Figure 2b₂). However, our data illustrate the details of the glacier deceleration phase, which were prolonged, spanning 5–6 years and gradually stabilizing in 2018–2019, with surface velocities mostly decreasing to below 0.3 m/d in most regions.

Figure 2. The interannual and daily scale velocity variations in both WWKG (a₁ and a₂) and EWKG (b₁ and b₂) and the average velocities and velocity magnitudes of different cross-sections within WWKG (a₃) and EWKG (b₃) during various stages of glacier surging. Reference Figure 1 for flowline and cross-section locations.

Compared to EWKG, WWKG exhibits higher velocities. We observed that WWKG initiated its surge during the summer of 2020, with an initial peak velocity of 1.32 ± 0.22 m/d (Figure 2a₂). In the summer of 2021, the glacier surface velocity reached its maximum at 2.61 ± 0.18 m/d (Figure 2a₂). Subsequently, during the autumn of 2022, glacier surface velocities decreased to below 1 m/d. As of the summer of 2023, the glacier surface velocities did not exceed 0.2 m/d (Figure 2a₂). Regarding its interannual velocity variations, WWKG's surge persisted for 4 years, with maximum speeds reaching 611 ± 4.33 m/yr, which is 1.45 times that of EWKG (Figure 2a₁). Figure 2a₃ illustrates the average velocity changes at different sections of WWKG, revealing that the glacier's maximum velocity magnitude is smaller than that of EWKG, even during the active phase. Additionally, it is noteworthy that the impact zone of the WWKG surge is primarily concentrated within 0–9.6 km from the glacier terminus, while EWKG exhibits a more extensive active area (Figure 2a₂ and b₂). Moreover, we observed that both glaciers, at approximately 10 and 14 km from their termini, show nearly zero glacier surface velocities during both the quiescent and active phases (Figure 2a₃ and b₃).

In contrast to optical imagery, synthetic aperture radar (SAR) data is less affected by weather conditions, enabling the capture of more detailed glacier surges. Leveraging the Sentinel-1 glacier velocity data provided by Friedl, Seehaus, and Braun (2021), this study extracted monthly and seasonal variations in surface velocity for the EWKG and WWKG. Our analysis confirms previous findings based on optical remote sensing imagery: EWKG exhibited a gradual deceleration of motion transitioning to a quiescent phase, while WWKG experienced a surge onset during the summer of 2020 (Figure 3a). From the seasonal variation of the surface velocity of the two glaciers, EWKG exhibited higher velocities during autumn, the deceleration phase, compared to other seasons (Figure 3c), whereas WWKG showcased elevated velocities during winter, the surge initiation phase, relative to other seasons (Figure 3b).

Figure 3. The monthly (a) and seasonal (b amd c) variations in surface velocity of the eastern (c) and western (b) branches of the West Kunlun glaciers from 2014 to 2021. Reference Figure 1 for flowline and cross-section locations. The data originate from Sentinel-1 satellite observations, as provided by Friedl, Seehaus, and Braun (2021).

4.2. Glacier surface elevation changes and mass transport

Apart from increased velocities, glacier surges often result in changes in surface elevation and the migration of ice volume. Prior to 2000, both the upper regions of WWKG and EWKG exhibited a thickening trend, while the glacier terminus regions experienced thinning, resulting in a lowering of glacier surface elevation by 27.58 ± 7.50 m and 72.91 ± 7.50 m, with corresponding mass losses of 0.086 ± 0.062 km³ and 0.112 ± 0.064 km³, respectively (Table 2 and Figure 4a). Entering the twenty-first century, there was a notable thickening of the EWKG terminus region (Figure 4b), with an average increase in glacier surface elevation of 23.49 ± 3.91 m and a corresponding mass gain of 0.155 ± 0.097 km³ (Table 2). This is a well-known indicator of glacier surging. In contrast, WWKG continued to thin, albeit at a significantly lower rate than before 2000 (Table 2). However, the differential results from the COP30 DEM and AW3D30 DEM showed a distinct change in the continuous thinning trend of WWKG since then (Figure 4c). Moreover, the reservoir area (Res. A) thickened by an average of 8.411 ± 0.43 m during this period, accompanied by a further reduction in the thinning rate in the receiving area (Rec. A) (Table 2). This finding provided the necessary material conditions for the initiation of WWKG surge. Following the triggering of the EWKG surge, glacier material continued to be transported from the Res.A. to the Rec.A., resulting in continued thickening at the glacier terminus (Figure 4c and Table 2). The WWKG glacier experienced a surge event in 2020. During this surge period, the glacier transported a total ice volume of 0.634 ± 0.135 km³ downstream, resulting in an average thickening of over 70 m at its terminus (Figure 4c and Table 2). In comparison, the EWKG glacier transported an accumulated ice volume exceeding 0.6 km³ (Table 2). Examining the elevation changes along the centrelines of both glaciers during their recent surge events, it is noteworthy that the maximum elevation changes at the glacier terminus exceeded 200 m (Figure 4e).

Figure 4. Changes in glacier surface elevations during the periods (a) KH9 DEM – NASADEM, (b) NASADEM – AW3D30 DEM, (c) AW3D30 DEM – COP30 DEM and (d) COP30 DEM – ASTER DEM, and vertical profiles of elevation change along the glacier centreline for the WWKG (e₁) and EWKG (e₂). Reference Figure 1c for flowline location.

Table 2. Mean surface elevation and ice volume changes in the Res.A. and Rec.A. of the WWKG and EWKG during the KH9 DEM – NASADEM, NASADEM – AW3D30 DEM, AW3D30 DEM – COP30 DEM, and COP30 DEM – ASTER DEM.

Glacier	WWKG				EWKG
DEMs	Res.A		Rec.A		Res.A		Rec.A
	Mean dH (m)	Volume (km^3)	Mean dH (m)	Volume (km^3)	Mean dH (m)	Volume (km^3)	Mean dH (m)	Volume (km^3)
KH9 DEM (1970s) – NASADEM (2000)	−0.717±7.5	−0.002±0.05	−27.58±7.5	−0.086±0.062	−2.901±7.5	−0.008±0.058	−42.91±7.5	−0.112±0.064
NASADEM (2000) – AW3D30 DEM (2006/2011)	−0.079±3.91	−0.0006±0.075	−3.154±3.91	−0.028±0.094	−0.65±3.91	−0.006±0.088	23.49±3.91	0.155±0.097
AW3D30 DEM (2006/2011) – COP30 DEM (2010/2015)	8.411±0.43	0.067±0.008	−0.649±0.43	−0.005±0.01	−8.843±0.43	−0.088±0.009	31.57±0.43	0.195±0.01
COP30 DEM (2010/2015) – ASTER DEM (2023)	−3.968±5.56	−0.027±0.109	70.25±5.56	0.634±0.135	6.99±5.56	0.06±0.126	31.25±5.56	0.252±0.138

4.3. The evolution of glacier termini and surface crevasses

The manually digitized glacier outlines reveal that both WWKG and EWKG experienced a significant advance at their terminus between 1977 and 2023. In 1977, WKG covered an area of 133.81 ± 0.84 km² (Figure 5b). However, due to rising regional temperatures, this glacier has been in a continuous state of retreat (Figure 5b). By 2005, WKG had experienced terminus separation (Figure 5a), resulting in a reduction in glacier area by over 5 km² compared to 1977 (Figure 5b). Subsequently, the EWKG underwent a surge, advancing its terminus by over 1000 m (Figure 5d). Furthermore, following the triggering of the WWKG surge in 2020, it advanced by 804 m (Figure 5c). By 2022, the termini of WWKG and EWKG reconnected, forming a single glacier once again. At this point, WKG covered an area of 132.42 ± 0.84 km² (Figure 5b), reflecting an increase by 3.65 ± 0.85 km² compared to 2005. However, despite this increase, it remains in a state of overall retreat when compared to the 1970s.

Figure 5. Glacier terminus (a) and its area (b) changes in WKG. Advancement of glacier terminus in WWKG (c) and EWKG (d). Figure (a), (c), and (d) base maps are derived from Band 8 of Landsat OLI and ETM + . Figure (b) Results obtained through manual digitization of Landsat data, with blue shadows indicating glacier area uncertainty.

Examining the evolution of glacier surface crevasses during glacier surges, it becomes evident that both WWKG and EWKG experienced pronounced crevasse expansion. Particularly noteworthy is the case of WWKG, where in August 2018, significant lateral crevasses had developed in the central part of the glacier (Figure 6a₁ and a₂). This finding indicated a sharp increase in instability at the glacier's base during this period. As the glacier flowed downstream under the influence of gravity, these crevasses gradually widened (Figure 6a₃). Furthermore, after WWKG entered its active phase in 2020, these crevasses expanded rapidly and extended to the lower portion of the glacier. Similarly, crevasses on the surface of EWKG became noticeably more abundant during its surge, albeit at a slower rate compared to WWKG. Moreover, these crevasses tended to be concentrated at the glacier's edges, leaving the central part of the glacier relatively smooth in appearance (Figure 6b).

Figure 6. Changes in glacier surface crevasses of WWKG and EWKG during the active phase. Except for the base maps in (b₁) and (b₂), which are sourced from Landsat OLI Band 8 (with a spatial resolution of 15 m), the base maps for the other subfigures are all derived from Sentinel-2 Band 4 (with a spatial resolution of 10 m).

5. Discussion

5.1. Uncertainties

The uncertainties in glacier surface elevation and velocity changes are critical factors influencing the reliability of our study's findings. In glaciology research, errors in results obtained from remote sensing are typically categorized as relative and absolute errors (Lv et al. 2019; Nolan et al. 2021). However, due to the remote and harsh environments of glaciers, most scholars evaluate the accuracy of their findings using relative or cross-validation methods. Therefore, following precedent, this study assumes that the terrain outside of the glacier areas remains stable (Shean et al. 2020). Then, we computed the mean elevation residual (MD), standard deviation (STDV), and the number of pixels after spatial autocorrelation removal for elevation changes (N) in stable terrain surrounding glaciers to calculate the standard error of the mean (SE). The formula is provided below: (11) $\mathrm{SE}=\mathrm{STDV}/\sqrt{N}$(11) We computed the overall error of glacier surface elevation changes based on the error propagation method. Detailed error information is presented in Table 2. (12) ${\delta }_h=\sqrt{M{D}^2+S{E}^2+{\delta }_P^2}$(12) In the equation, δ_P denotes the penetration depth of NASADEM into ice and snow. NASADEM was derived from data obtained through the SRTM-C band (Nolan et al. 2021). In this study, we performed a penetration depth correction using X-band data acquired concurrently with the C-band data. The detailed errors of glacier surface elevation changes at different time periods are presented in Table 3 of this study.

Table 3. Statistics of vertical errors between multisource DEMs.

Item	MD (m)	STDV (m)	N	SE (m)	δP (m)	δh (m)
NASADEM-KH9 DEM	3.95	11.62	575	0.48	3.52	7.50
AW3D30 DEM-NASADEM	−0.12	7.77	411	0.38	3.52	3.91
COP30 DEM-AW3D30 DEM	0.34	5.44	411	0.27	–	0.43
ASTER DEM-COP30 DEM	5.43	23.94	411	1.18	–	5.56

In current research, image cross-correlation and offset tracking are automated methods employed to identify and track glacier surface displacements before and after their occurrence based on optical or SAR data. These methods extract similar features between two image scenes to derive glacier surface velocity, with theoretical observation precision reaching 1/10 of a pixel size. In this study, we employed the same approach as glacier elevation change error assessment to evaluate the uncertainty of the obtained velocity data. Our findings indicate that the errors in daily scale velocity range from 0.05 m/d to 0.59 m/d, while errors at the interannual scale range from 1.01 m/yr to 6.3 m/yr (Figure 7). In addition, to assess the accuracy of surface velocity extracted from the WKG using Landsat OLI imagery, we initially compared it with velocities derived from Sentinel-2 imagery obtained during the same period. We observed a correlation of 0.88 between the two datasets, with a root mean square error of 0.015 (Figure 8). Consequently, we thought that Landsat OLI demonstrates high reliability in capturing glacier surface velocities.

Figure 7. Displacement on stable terrain observed in (a) daily scale velocity and (b) interannual scale velocity. Box plots illustrate the overall distribution of horizontal displacements in stable terrain areas during the respective periods. The dashed line represents the overall mean.

Figure 8. Comparison of glacier surface velocity between Landsat and Sentinel-2.

5.2. WWKG and EWKG surge characteristics

Recently, many scholars have summarized the characteristics of glacier surge events controlled by different glacier surge mechanisms. Based on these characteristics, glacier surges have been classified into two main categories: Alaskan-type STGs and Svalbard-type STGs (Quincey et al. 2015; Yasuda and Furuya 2015). Some scholars argue that Alaskan-type STGs are influenced by glacier internal hydrological conditions. The acceleration and deceleration of such surges occur rapidly, typically initiating in the winter when the glacier internal drainage system is obstructed and concluding in the summer when the drainage system regains its efficiency (Kamb et al. 1985). On the other hand, Svalbard-type STGs are believed to be influenced by changes in glacier thermodynamic conditions. These surges exhibit a delayed response, with the acceleration phase peaking several years after initiation and the deceleration phase extending over several years, with the surge's onset and termination being less constrained by seasonal conditions (Frappé and Clarke 2007; Murray et al. 2003).

In this study, we obtained comprehensive temporal characteristics of WWKG and EWKG surges using high temporal resolution data at both annual and daily scales. We found that the WWKG surge commenced in the summer of 2020, with a noticeable deceleration in glacier surface velocity during the winter of 2021. Upon merging with EWKG, the glacier surface velocity returned to presurge levels. In terms of annual velocity changes, the WWKG glacier displayed a 1-year acceleration phase followed by a 3-year deceleration phase, mirroring the surge characteristics of Alaskan-type STGs. In contrast, the surge characteristics of EWKG closely resemble those of Svalbard-type STGs. This glacier experienced a surge lasting over 10 years, with an 8-year deceleration phase. Examination of the glacier daily velocity revealed that the surge in EWKG terminated in the summer. In comparison with typical STGs in other regions of HMA, the surge behaviour of EWKG was akin to that of cold-type glaciers, such as Monomah Glacier in the Eastern Kunlun Mountains, as described by Guo et al. (2020), which are controlled by thermal mechanisms. Conversely, WWKG's surge characteristics align with glaciers, such as South Rimo Glacier (Jiang et al., 2021), Shisper Glacier (Beaud et al. 2022), Kyager Glacier (Bazai et al. 2021) in the Karakoram, Karayaylak Glacier (Zhang et al. 2022) in the Eastern Pamirs, and Sabche Glacier (Lovell, Carr, and Stokes 2018) in the Himalayas. Notably, research indicates that the surge of these glaciers was governed by subglacial hydrological conditions. In summary, based on the surge characteristics of WWKG and EWKG, we propose that the EWKG was controlled by thermal mechanisms, while the western branch was influenced by hydraulic mechanisms. However, it is essential to acknowledge that the seasonality of glacier surges in HMA differs from that of Alaskan-type and Svalbard-type STGs. Therefore, we propose that, in addition to hydrothermal conditions, other factors, such as climate and topography, play significant roles in governing the asynchronous surge behaviour of these two glaciers Table 4.

Table 4. Summary of the time characteristics of glacier surges in different STGs.

Glacier	Duration of glacier surge (yr)	Acceleration duration of glacier surge (yr)	Deceleration duration of glacier surge (yr)	Initial season of glacier surge	End season of glacier surge	References
Alaskan-type STGs	1–3	Rapidly	Rapidly	winter	summer	Quincey et al. (2015)
Svalbard-type STGs	>3	Slow	Slow	any season	any season	Quincey et al. (2015)
Sabche Glacier	3–5	3	2	winter	–	Lovell, Carr, and Stokes (2018)
South Rimo Glacier	2	1	1	summer	summer	Jiang et al. (2021)
Shisper Glacier	3	1	2	Spring	Spring	Beaud et al. (2022)
Monomah Glacier	9	5	4	summer	Spring	Guo et al. (2020)
Kyager Glacier	3	1	2	Spring	summer	Bazai et al. (2021)
Karayaylak Glacier	1.5	1	0.5	Spring	Autumn	Zhang et al. (2022)
WWKG	4	1	3	summer	winter	In this study
EWKG	20	12	8	–	summer	In this study

5.3. Influence of climate change on surge behaviours

Previous research has indicated that while climate change is not the fundamental trigger for glacier surges, it is the outcome of disturbances in climate/mass balance that are subsequently adjusted through glacier dynamics (Farinotti et al. 2020). Recent studies have further demonstrated the significant influence of climate change on glacier surge periodicity and characteristics (Guan et al. 2021). Consequently, changes in glacier thermal and hydrological states are inevitably influenced by meteorological factors. In this study, we observed a clear decreasing trend in precipitation in this region before 2000 (Figure 9a). This finding resulted in negative mass balances of −0.10 ± 0.03 m w.e.a⁻¹ and −0.15 ± 0.05 m w.e.a⁻¹ for WWKG and EWKG, respectively (Figure 9a). After 2000, the mass balance of EWKG became 0.56 ± 0.23 m w.e.a⁻¹ (Figure 9a). Consequently, the exceptional acceleration of EWKG motion during this period can be understood (Figure 2b₁), as the increase in ice thickness led to elevated basal stresses and a basal state near the pressure melting point, a major driver of velocity increase (Cuffey and Paterson 2010). In contrast, WWKG exhibited a near-zero mass balance during this time (Figure 7a). However, with increasing regional precipitation, the glacier's mass accumulation rate rapidly rose, resulting in a mass balance twice that of EWKG. Nevertheless, WWKG did not surge in approximately 2010, whereas EWKG experienced a rapid increase in velocity in 2011 (Figure 2b₁). Undoubtedly, this was a response to the rapid increase in glacier mass. After 2013, there was a noticeable increase in cumulative precipitation in the region, with mass balances reaching 1.45 ± 0.13 m w.e.a⁻¹ and 0.99 ± 0.18 m w.e.a⁻¹ for WWKG and EWKG, respectively (Figure 9a). This changing trend directly triggered the surge of WWKG. In conclusion, we contend that under the backdrop of climate change, the rate and intensity of glacier mass accumulation in the region are driving factors in the asynchronous surge behaviour of neighbouring glaciers.

Figure 9. (a) Anomalous precipitation and mass balance variations in the WWKG and EWKG regions from 1981 to 2022. The relationship between glacier surface albedo and velocity for WWKG (b) and EWKG (c) from 2013 to 2023. In the figure, the solid blue line represents the velocity, and the solid red line represents the albedo.

On the other hand, studies have underscored the pivotal role of water in glacier motion, as basal sliding is a fundamental component of glacier motion (Cuffey and Paterson 2010). The velocity of basal sliding depends on basal shear stress and bed roughness as well as on the magnitude of water pressure at the glacier base (Benn et al. 2019). This water primarily originates from various sources, including surface meltwater, rainfall, internal meltwater, and basal meltwater (Guo et al. 2022). Among these, surface meltwater constitutes the predominant water source for glacier sliding (Yasuda and Furuya 2015). Nevertheless, quantifying its impact on glacier dynamic processes remains exceedingly challenging without direct field observations. Therefore, this study employs high temporal-resolution relationships between glacier surface velocity and albedo to qualitatively unveil the influence of surface meltwater on WWKG and EWKG in the context of climate change. We discovered a strong positive correlation between albedo and glacier velocity in WWKG before the glacier surged, signifying that as albedo increased, glacier velocity likewise rose. However, during the active phase, a distinct negative correlation emerged between the two (Figure 9b). Although we could not obtain complete high-temporal-resolution velocity data for EWKG during its active phase, its characteristics appeared similar to those of WWKG, as evidenced by the relationship between albedo and glacier velocity during its deceleration phase (Figure 9c). Glacier surface albedo serves as an indicator of surface accumulation and ablation (Di Mauro and Fugazza 2022). An increase in albedo implies a larger snow-covered area on the glacier surface, while a decrease indicates an elevated shortwave radiation received by the glacier surface (Cuffey and Paterson 2010). Consequently, based on the Glen flow law (Glen 1955), we infer that during the quiescent phase, glacier dynamics for both WWKG and EWKG are predominantly influenced by glacier ice deformation. In contrast, during the active phase, considering the increased fragmentation and the presence of crevasses on both glacier surfaces (Figure 6), glacier dynamics are primarily driven by basal sliding induced by surface meltwater.

5.4. Outlook on surge behaviour

Based on the dynamic evolution characteristics of WWKG and EWKG (Figure 2), we observe that the surging behaviour of WWKG aligns closely with hydrologically controlled Alaskan-type STGs (Kamb et al. 1985), whereas the surging behaviour of EWKG exhibits similarities with thermally conditioned Svalbard-type STGs (Murray et al. 2003). Additionally, according to prior research, till deformation plays a crucial role in glacier surges (Harrison and Post 2003). However, due to the challenges in on-site observations of subglacial/subice processes, we are unable to quantitatively expound on their impact on glacier surging. Nevertheless, according to Coulomb-plastic rheology for till deformation, elevated basal water pressure is a necessary condition for substantial till deformation (Harrison and Post 2003). Therefore, the surging mechanisms of WWKG and EWKG should involve a process capable of providing high basal water pressure. Fortunately, the relationship between glacier velocity and surface albedo provided in this study (Figure 9b and c) demonstrates a strong correlation between the initiation and termination of surges in WWKG and EWKG with surface ablation. Hence, undoubtedly, surface meltwater, facilitated through crevasses in the ice surface (Figure 6), provides an ample water source for the elevation of subglacial/subice water pressure. This aligns with the observations of Yasuda and Furuya (2015). Further analysis reveals that the velocities of WWKG and EWKG are directly proportional to surface albedo during the quiescent phase. According to prior research, this suggests a logical sequence of increased solid precipitation – higher albedo – glacier thickening – enhanced deformational speeds – accelerated flow velocities (Guo et al. 2020). Thus, we infer that during the quiescent phase, the glacier motion of WWKG and EWKG is dominated by glacier creep influenced by glacier thickness.

However, this still does not explain why, under similar climatic conditions, the surges of the two glaciers are asynchronous and exhibit such significant differences in duration. This study investigates the evolution of the surface slope and subglacial topography of WWKG and EWKG using open-source DEM, and glacier thickness simulation data provided by Farinotti et al. (2019). We observe that in 2000, the upper slopes of the WWKG and EWKG surfaces are both less than those at the glacier terminus (Figure 10a₁). As glacier mass increases (Figure 9a), the slope at the glacier terminus gradually increases and exhibits a discontinuous spatial pattern (Figure 10a₂). Based on COP30DEM, under the context of a rapid increase in glacier surface mass, the surface slope transitions to a discontinuous state from top to bottom. However, in comparison to EWKG, the upper slopes of WWKG are consistently less than 3° and exhibit a more continuous spatial distribution. Therefore, considering the start-stop times of surges for WWKG and EWKG and the shallow ice approximation theory, we conclude that the asynchronous behaviour of the two glaciers during surges is controlled by subglacial topography. At the glacier scale, the upper and lower terrains of WWKG are notably flatter in comparison (Figure 10b₁), and due to their relatively larger size, they possess a larger accumulation area than EWKG. Consequently, under similar climatic conditions, this glacier's accumulation area can accommodate more glacier mass, maintaining a more stable state. As mass rapidly accumulates in the accumulation area and the glacier transports mass downstream under the influence of gravity, variations in glacier thickness result in changes to basal stress, leading to a fracture in the glacier approximately 7 km from the terminus (Figure 6a). Subsequently, surface meltwater enters the glacier base, increasing basal water pressure. In contrast, the bottom terrain of the EWKG exhibits a distinct serrated pattern (Figure 10b₂). For such subglacial topography, basal meltwater in the glacier must flow around obstacles to the downslope side with lower pressure and higher-pressure melting points, requiring a longer time to raise basal water pressure to the threshold triggering glacier surges. Consequently, although EWKG entered an unstable state in 1999, actual acceleration occurred in 2010 (Figure 2b₁). Therefore, subglacial topography determines the spatial distribution of internal glacier stress and the rate at which basal water pressure increases, leading to significant differences in surge behaviour and characteristics, despite similar glacier terrain and climatic backgrounds.

Figure 10. Spatial evolution of glacier surface slopes during glacier surges (a) and slope and subglacial topography of centreline profiles in WWKG (b₁) and EWKG (b₂). The vertical red dashed line in the figure denotes the location of the fracture in WWKG. The glacier bed was derived by subtracting NASADEM from Farinotti et al. (2019), providing support for certain aspects of our study.

6. Conclusion

In investigating the asynchronous controls on the surging behaviours of the eastern and western branches of the WKG glacier, we derived temporal characteristics, velocity variations, elevation changes, and surface morphological evolution during glacier surges using high-temporal-resolution multisource DEM and glacier surface velocity data. We observed that under similar climatic conditions, the surge in EWKG preceded that of WWKG, persisting for more than 20 years. Considering its acceleration and deceleration phases, we posit similarities between the surge characteristics of EWKG and those controlled by thermodynamic mechanisms, akin to Svalbard-type STGs. In contrast, WWKG's surge commenced in the summer of 2020, exhibiting a relatively rapid acceleration and deceleration process reminiscent of surges governed by hydraulic mechanisms, akin to Alaskan-type STGs. Analysing the relationships between glacier mass accumulation rates and precipitation, as well as albedo and flow velocity, revealed a positive correlation between the initiation time of glacier surges and the rate of glacier mass accumulation. Furthermore, during quiescent phases, glacier albedo was proportional to surface velocity, while during active phases, it exhibited an inverse relationship. We propose that the gradual transport of surface meltwater into glacier beds, accompanying the development of surface crevasses, promotes an increase in subglacial water pressure. Additionally, we hypothesize that subglacial topography controls the initiation time and duration of surges in WWKG and EWKG glaciers by influencing spatial variations in subglacial shear stress due to changes in glacier thickness and the rate of increase in subglacial water pressure.

Acknowledgements

We thank the United States Geological Survey (USGS) for the Landsat data (https://www.usgs.gov/), the National Aeronautics and Space Administration (NASA) for the NASADEM, ITS_LIVE, MCD43A3 and ASTER-L1A (https://cmr.earthdata.nasa.gov/), Japan Aerospace Exploration Agency (JAXA) for the AW3D30 DEM (https://www.eorc.jaxa.jp/), the European Space Agency (ESA) for the COP30 DEM (https://spacedata.copernicus.eu/), and University of California, Santa Barbara for the CHIRPS (https://chc.ucsb.edu/data/chirps).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study can be obtained from the corresponding author upon reasonable request.

Additional information

Funding

This study was supported by the Top-notch Talent of the Qinghai Province ‘Kunlun Talent. High-end Innovation and Entrepreneurship Talent’ program [grant no 2023-QLGKLYCZX-001], the Postdoctoral research project, Faculty of Geography, Yunnan Normal University [grant number 01300205020516079], the Project Supported by Key Laboratory of Resource Environment and Sustainable Development of Oasis, Gansu Province [grant no GORS202302], the Enhancement Program for Research Capabilities in the Faculty of Geography [grant no 01300205020516083/002], the Postdoctoral research start-up project of Yunnan Normal University [grant number 01300205020503329], and the Researcher Development program of Qinghai University of Science and Technology [grant number 2023011wys007].