Velocity variations and hydrological drainage at Baltoro Glacier, Pakistan

Glacial meltwater directly influences glacier dynamics. However, in the case of debris-covered glaciers, the drivers of glacier velocity and the influence of supraglacial lakes have not yet been sufficiently analysed and understood. We present a spatio-temporal analysis of key glacier characteristics for Baltoro Glacier in the Karakoram from October 2016 to September 2022 based on Earth observation data and climate parameters extracted from the High Asia Refined analysis (HAR) data set. For the glacier variables, we used surface velocity, supraglacial lake extent, melt of snow and ice, and proglacial run-off index. For climate variables, we focused on air temperature and precipitation. The surface velocity of Baltoro Glacier was characterized by a spring speed-up, summer peak, and fall speed-up with a relative increase in summer of 0.2–0.3 m d⁻¹ (75 %–100 %) in relation to winter velocities, triggered by the onset of or an increase in basal sliding. Snow and ice melt have the largest impact on the spring speed-up, summer velocity peak, and the transition from inefficient to efficient subglacial drainage. The melt covered up to 64 % (353 km²) of the entirety (debris-covered and debris-free) of Baltoro Glacier and reached up to 4700 m a.s.l. during the first melt peak and up to 5600 m a.s.l. during summer. The temporal delay between the initial peak of seasonal melt and the first relative velocity maximum decreases downglacier. Drainage from supraglacial lakes (3.6–5.9 km²) contributed to the fall speed-up, which showed a 0.1–0.2 m d⁻¹ (20 %–30 %) lower magnitude compared to the summer velocity peak. Most of the run-off can be attributed to the melt of snow and ice. However, from mid-June onward, the lakes play an increasing role, even though their contribution is estimated to be only about half of that of the melt. The observed increase in summer air temperatures leads to a greater extent of melt, as well as to a rise in the number and total area of supraglacial lakes. This tendency is expected to intensify in a future warming climate.

1 Introduction

Glacial meltwater is an important control over glacier dynamics. In particular, the seasonal evolution and variation of glacier flow are strongly influenced by the timing and amount of meltwater (Glasser, 2013; Iken and Bindschadler, 1986). With the melt onset, surface meltwater is formed and possibly drains to the glacier bed via crevasses or englacial pathways and is stored in cavities and distributed channels (Röthlisberger, 1972; Cuffey and Paterson, 2010). This influx of meltwater into the subglacial drainage system leads to an increase in basal water pressure. When subglacial water pressure approaches ice-overburden pressure, basal traction decreases and sliding increases as the ice decouples from the bed. The glacier – in the case of a hard glacier bed – lifts up and accelerates (Weertman, 1964; Lliboutry, 1968; Iken and Bindschadler, 1986; Nolan and Echelmeyer, 1999; Sugiyama et al., 2011; Hoffman et al., 2016; Benn et al., 2019). Additionally, the inflow of water through small and incipient subglacial channels generates frictional heat which melts the ice walls and thus expands the channels, leading ultimately to the formation of an efficient drainage system (Röthlisberger, 1972; Flowers, 2015). In an efficient subglacial channel system, water storage is reduced, water pressure in the channels lessens, and glacier flow velocity decreases (Benn et al., 2019). In the absence of meltwater, the effective pressure is larger and leads to a closure of the channels through regelation (Cuffey and Paterson, 2010) and creep (Benn et al., 2019; Flowers, 2015; Jiskoot, 2011).

Glacier movement is governed by gravitational driving stress and balanced by the resistive stresses. The driving stress is influenced by gravitational acceleration, ice density, thickness, and surface slope and varies slowly. The resistive stresses arise from the drag at the glacier and by means of dynamical flow resistance. Glacier flow processes can be divided into ice deformation, basal sliding, and sediment deformation. Sliding and bed deformation occur only in the case of temperate and polythermal glaciers, with a higher contribution by basal sliding than internal deformation velocities (Boulton and Hindmarsh, 1987; Jiskoot, 2011).

Glacier speed-ups in spring (Mair et al., 2001; Macgregor et al., 2005; Nanni et al., 2023), summer (Iken and Bindschadler, 1986; Copland et al., 2003; Bartholomaus et al., 2008; Quincey et al., 2009; Hewitt, 2013; Werder et al., 2013; Van Wychen et al., 2014; Armstrong et al., 2017; Nanni et al., 2023; Rada Giacaman and Schoof, 2023), and winter (Burgess et al., 2013; Hart et al., 2022) have been widely observed, and their changes have been linked to an increase in air temperature, surface melt, and subglacial hydrology. In the case of debris-free glaciers, meltwater creation can be directly associated with warmer air temperatures and surface melt. Debris-covered glaciers, however, have a more complex, non-linear ablation with enhanced melting in areas of thin debris cover (a few centimetres), thermal insulation in areas of thick debris coverage (Østrem, 1959; Nicholson and Benn, 2006), and melting hotspots at ice cliffs (Brun et al., 2018; Buri et al., 2021) and in supraglacial lakes (Miles et al., 2020). In an inefficient subglacial drainage system, supraglacial lake discharge can support basal sliding and hence cause higher glacier velocities (Sakai and Fujita, 2006; Sakai, 2012; Watson et al., 2016; Benn et al., 2017; Miles et al., 2020). A few studies have observed that supraglacial lake drainage could lead to a transition from an inefficient to an efficient subglacial drainage system and hence lead to a slowdown of the glacier velocity (Vincent and Moreau, 2016; Stevens et al., 2022).

Previously, we presented a time series of annual and seasonal glacier surface velocities derived from multi-mission synthetic aperture radar (SAR) data for Baltoro Glacier, located in the Karakoram, Pakistan, from 1992 to 2017 (Quincey et al., 2009; Wendleder et al., 2018). We could show that, in some years, the acceleration lasted longer and affected a larger glacier area than in others. In years with higher velocities, the supraglacial lakes mapped from Landsat and ASTER imagery were characterized by a larger number and a larger total area as well. However, only one image for each summer was available for mapping, and this did not provide sufficient insight into the seasonal evolution of the supraglacial lakes or their link to surface velocity. Therefore, we developed an approach using multi-temporal and multi-sensor Earth observation data to provide a dense, almost daily summer time series of supraglacial lakes (Wendleder et al., 2021a). Nevertheless, there is still a lack of detailed process understanding of whether and how the development of supraglacial lakes triggers basal sliding.

In this study, we combined relevant glacier variables derived from Earth observation data and climate records to assess the extent to which Baltoro Glacier responds, spatially or temporally, to a given climatic forcing. The processes and relationships of the variables were examined by statistical analysis. For the glacier variables, we used (1) surface velocity derived by intensity offset tracking from Sentinel-1 time series; (2) supraglacial lakes mapped by a random forest classifier applied to Sentinel-2, PlanetScope, Sentinel-1, and TerraSAR-X data; (3) melt of snow and ice detected using a change detection algorithm based on Sentinel-1; and (4) run-off index estimated as surface areal coverage of the proglacial stream (given in km²) and used as a proxy of the quantitative run-off from Sentinel-2 and PlanetScope imagery. For the climate variables, we focussed on air temperature and precipitation extracted from the High Asia Refined analysis (HAR) data set. The analysis focused on the period from October 2016 to September 2022, providing a dense and continuous time series.

2 Study site

Baltoro Glacier is located in the eastern Karakoram in the northern part of Pakistan (Fig. 1). The glacier has a length of about 63 km and, together with its tributary glaciers, covers an area of approximately 554 km². Above Concordia (4600 m a.s.l.), the two major tributaries, Godwin Austen Glacier and Baltoro South Glacier, converge to the main Baltoro Glacier and change their flow direction westward. Several major tributary glaciers merge with the main branch along its northern and southern margins (Mayer et al., 2006; Quincey et al., 2009). Surface velocities range between 0.5 to 0.65 m d⁻¹ in summer and 0.3 to 0.4 m d⁻¹ in winter between Concordia and Urdukas (3900 m a.s.l.) and decrease to 0.03 to 0.1 m d⁻¹ near the terminus at 3400 m a.s.l. (Quincey et al., 2009; Wendleder et al., 2018). Approximately 38 % of Baltoro Glacier is debris covered, with a thin layer of 5–15 cm at Concordia, 30–40 cm at Urdukas, and on the order 1 m thickness near the glacier terminus (Mayer et al., 2006; Quincey et al., 2009). Furthermore, the debris thickness varies spatially due to advection, debris and meltwater movement, and slow cycles of topographic inversion (Nicholson et al., 2018; Huo et al., 2021), which impedes the determination of the area proportion enhancing surface melt. Debris-covered glaciers are characterized by the presence of ice cliffs (Mayer et al., 2006; Evatt et al., 2017) and supraglacial lakes (Wendleder et al., 2021a). On Baltoro Glacier, ice cliffs are found between the terminus and Gore. Supraglacial lakes are located on the main glacier from the terminus up to Concordia. They usually fill between mid-April to mid-June and drain between mid-June to mid-September.

The Karakoram has a mid-latitude, high-mountain climate with cold winters and mild summers. The regional climate is influenced by winter westerly disturbances with dominant winter and spring snowfall (Mayer et al., 2014; Dobreva et al., 2017); Indian summer monsoons with higher liquid precipitation, temperatures, and cloud coverage (Thayyen and Gergan, 2010; Bookhagen and Burbank, 2006); and the predominantly stable Tibetan anticyclone. In the case of an irregular weakening of the Tibetan anticyclone, which thus causes an incursion of the Indian summer monsoon, large amounts of summer precipitation can be observed (Dobreva et al., 2017). The climate variables are strongly determined by the topography. Hence, precipitation increases with altitude, reaching a mean annual rate of approximately 1600 mm yr⁻¹ at 5300 m a.s.l. (Godwin Austen region) and at 5500 m a.s.l. (Baltoro South region). The average daytime air temperature during summer is close to the freezing point at 5400 m a.s.l.; consequently, most of the precipitation deposits as snow above this elevation (Mayer et al., 2006).

Figure 1 Overview of Baltoro Glacier, Pakistan. Important place names are indicated. Red crosses mark the locations at Urdukas, Yermanendu, Gore, and Concordia, where the glacier surface velocity was extracted. The red point indicates the location of Urdukas, which differs by 3 km from the extracted location for the glacier velocity. The area of the mapped run-off is displayed with a white polygon. The western, central, and eastern sectors for the mapping of the snow and ice melt are displayed in yellow, blue, and orange, respectively. The air temperature was extracted at Urdukas, and the precipitation was extracted at Baltoro South Glacier, which is marked with a yellow point (centre of the HAR pixel). The background image is a Sentinel-2 composite (shortwave infrared, near-infrared, red) acquired on 22 July 2019.

3 Material and methods

This section describes the used Earth observation data, known methods to process them into level-2 data, and the used reanalysis data. Furthermore, the analyses of the spatio-temporal relationships of the variables, the Pearson correlation, and the linear regression to quantify the direction and strength of the relationship and the dependency are explained.

3.1 Glacier surface velocity

The glacier surface velocity was calculated from the Sentinel-1 Interferometric Wide (IW) swath single-look complex (SLC) data with a spatial resolution of 10 m (Table A1 in the Appendix). We only used data from the ascending orbit as only these were continuously available for the complete observation period. The intensity offset tracking algorithm was applied to consecutive pairs of co-registered SAR intensity images (Strozzi et al., 2002; Friedl et al., 2018). Feature tracking and subsequent processing steps were implemented in GAMMA (release 20211201, Wegmüller et al., 2016). Tracking patch sizes and step sizes were adapted to sensor specifications and expected displacement, lying at 250×50 pixels and 50×10 pixels, respectively. The Copernicus digital elevation model (release DEM, 30 m, version 2022_1, ESA, 2019) provided the topographic reference for geocoding and orthorectification of the surface velocity maps. The mean accuracy of the velocity maps is 0.06 m d⁻¹, resulting in a standard deviation of 0.042 m d⁻¹, which is calculated according to Strozzi et al. (2002) and Friedl et al. (2021). For the analysis, we extracted the surface velocity values along the glacier centreline at four different locations (Fig. 1), namely at Urdukas (11 km), at the confluence with Yermanendu Glacier (16 km), at Gore (23 km), and at Concordia (33 km distance from terminus). The four different points were selected to best reflect the spatial variation along the glacier.

3.2 Supraglacial lake mapping

Supraglacial lakes were mapped based on a multi-sensor and multi-temporal summer time series from 2016 to 2022, acquired by the optical sensors Sentinel-2 and PlanetScope, as well as by the SAR sensors Sentinel-1 and TerraSAR-X. The Sentinel-2 Multi Spectral Instrument (MSI) orthorectified level-1C top-of-atmosphere products (spatial resolution of 10 m) were atmospherically corrected to L2A products using MAJA (MACCS-ATCOR Joint Algorithm, release 4.2, Hagolle et al., 2017). The PlanetScope Analytic Ortho Scene products (level 3B, spatial resolution of 3 m) were downloaded as orthorectified and atmospherically corrected surface reflectance (SR) data. In our approach, a coregistration of the PlanetScope data is not needed since the assignment of the classified lakes over the season is performed using their centre coordinate. Sentinel-1 IW swath SLC (spatial resolution of 10 m) C-band and TerraSAR-X ScanSAR Multi Look Ground Range Detected (MGD, spatial resolution of 40 m) X-band data were processed to Analysis Ready Data (ARD) using the Multi-SAR system (Schmitt et al., 2015, 2020). The mapping of the seasonal lake evolution used a semi-automatic approach which is based on a random forest classifier applied separately to each sensor. To produce a consistent and internally robust time series, a combination of linear regression and the Hausdorff distance (Hausdorff, 1914) was used to harmonize SAR- and optical-derived lake areas. The advantage of the multi-sensor approach is that it highlights the strengths of each sensor and compensates for their weaknesses: Sentinel-2 provides a continuous, radiometrically stable time series with a temporal resolution of 12 d, which is filled by the high temporal sampling of PlanetScope data. During periods of cloud cover, SAR data provide important information. The disadvantages of the SAR data, i.e. lake area underestimation and missing data from side-looking radar geometry and undulating glacier surface topographies, can be compensated for by using optical data acquired on the same day. The time series has a temporal sampling of 2–4 d with a spatial sampling of 10 m. The mean relative root-squared error (RSE) is at 1.0 % (total area of 9.151 km² with an absolute RSE of 0.0945 km²). Detailed processing steps, as well as the distribution and lifespan of the supraglacial lakes for the years 2016 to 2020, were published by Wendleder et al. (2021b). For this study, the time series was extended to 2021 and 2022. The filling and draining of the supraglacial lakes are presented by their aggregated area.

3.3 Glacier run-off mapping

Since we only had Earth observation data available, we were limited to mapping the width of the glacier discharge stream and using the estimated surface areal coverage of the proglacial stream as a proxy of the quantitative run-off (given in km²). To avoid misunderstanding conventional run-off measurements given in volume, we are using the term run-off index. The mapping was only applied to the area of the upper proglacial stream (2.7 km² area and 5.6 km length) bordered to the south and north by alluvial fans and slopes (see white polygon in Fig. 1). From the low-flow to peak-discharge period, the water-filled channel width at the glacier's terminus increased 10-fold from 70 to around 700 m. We used the Sentinel-2 Multi Spectral Instrument (MSI) orthorectified and MAJA atmospherically corrected L2A products (spatial resolution of 10 m) and the PlanetScope Analytic Ortho Scene products (level 3B) surface reflectance (SR) data (spatial resolution of 3 m). Due to the high turbidity of the glacier discharge, water and sandy soil are hard to differentiate using the normalized different water index. Therefore, the optical data of each sensor were first fused to Kennaugh elements using the red, blue, and near-infrared bands. In radar, the Kennaugh elements describe the polarimetric information and enable the interpretation of physical scattering mechanisms. In the case of optical data, the Kennaugh elements are a fusion technique to combine multiple multi-spectral bands into one image while enhancing their spectral characteristics (Schmitt et al., 2020). Afterward, a K-means clustering grouped each image into the two classes “run-off” and “background”. The run-off class was assigned by a fixed, manually selected point on the western border of the mapping area that was always covered with water. To evaluate the accuracy, four reference data sets with different acquisition dates for the run-off area were digitized manually on the basis of the near-infrared band. The classification achieved a mean relative RSE of 0.12 % (total area of 3.1 km² with an absolute RSE of 0.14 km²). The time series has a temporal sampling of 5–15 d in the months of May and June, which are characterized by a higher cloud coverage, and 1–5 d in the months of July to September with less cloud coverage.

3.4 Mapping of snow and ice melt

We mapped the snow and ice melt for the entirety of Baltoro Glacier (see outline of the three polygons in Fig. 1), which includes the debris-covered and debris-free parts of the glacier from the Sentinel-1 IW swath C-band data with a spatial resolution of 10 m. The processing of the SLC to ARD data was performed with the Multi-SAR system (Schmitt et al., 2020). The additional use of image enhancement (specifically, multi-scale multi-looking) resulted in a smoothing of noise and hence a more homogeneous environment and a better classification (Schmitt et al., 2015; Wendleder et al., 2022). As cross-polarization (VH) has a greater absorption over wet surfaces than co-polarization (VV), leading to a lower backscatter and hence a better discrimination of wet surfaces (Rott and Mätzler, 1987), we used only this polarization. For every scene, the image difference in relation to a reference scene was calculated. To ensure cold air temperatures prevailed during the acquisition, the first scene in January of each year, the coldest month, was chosen as the reference scene. As the threshold, we used one-half of the signal power (−3 dB) (Shi and Dozier, 1995; Nagler and Rott, 2000; Scher et al., 2021). This approach achieves an overall accuracy of 81.4 %–85.4 % (Nagler and Rott, 2000). By intersecting the mapping of the melt with the Copernicus DEM (1 arcsec, version 2022_1, ESA, 2019), the 10th and 90th percentiles of the elevation of the aggregated melt area were used, changing over time. We selected only SAR acquisitions from the ascending orbit. Firstly, the ascending images were acquired during the afternoon (13:00 UTC, 18:00 local time) to better represent the maximum melt area compared to the descending images acquired during the morning (01:00 h UTC, 06:00 h local time). Secondly, the accumulation areas of the tributary glaciers were imaged with only minor influence of layover or radar shadow as the slopes faced away from the sensor (Wendleder et al., 2021a). The mapping has a temporal resolution of 12 d. For the mapping, we divided the glacier into the western (from glacier terminus to Urdukas), central (from Urdukas to Gore), and eastern areas (upward from Gore). As the area and elevation of the glacier increase from west to east, the temporal delay of the melt in the three sectors reflects the vertical gradient.

3.5 Air temperature and precipitation

The near-surface air temperature at 12:00 (local time) and the total precipitation data (daily sum) were obtained from the daily interpolated HAR data set (version 2) provided by the Chair of Climatology, TU Berlin. The air temperature data (2 m height) from an automatic weather station located at Urdukas (Mihalcea et al., 2008) were needed for the local downscaling of the HAR data. The weather station is equipped with an air temperature sensor (thermo-hygrometer with radiation shield) and took hourly measurements in 2011. The global European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-5, spatial resolution of 80 km) was dynamically downscaled by regional climate models for High Mountain Asia to produce the regional refined HAR data set with a spatial resolution of 10 km. HAR is available for the period 1980 to 2022 (Wang et al., 2021). Compared with in situ data, the modelled air temperatures show a mean bias of −0.58 K (Wang et al., 2021). The elevation model used for the downscaling, with its spatial resolution of 10 km, is too coarse to represent the complex terrain of the Karakoram, which leads to an underestimation of air temperatures in the lower glacier and an overestimation of temperatures on the accumulation area. To downscale the air temperatures, a mean lapse rate for each month was determined based on air temperature data from the automatic weather station using a fitting curve (third order). The derived lapse rates varied between 7 °C km⁻¹ (December) and 11 °C km⁻¹ (June). The precipitation was extracted at Baltoro South Glacier, one of the two major tributaries within the larger glacier area.

3.6 Temporal relationship

The first step is to analyse if there is a temporal relationship between the supraglacial lakes, melt of snow and ice, glacier velocity, run-off index, air temperature, and precipitation. For the presentation, we plotted all variables in one figure using the hydrological year (1 October–30 September) as a time reference.

3.7 Spatial relationship

If the temporal representation shows a relationship in time between the variables of supraglacial lakes, melt of snow and ice, and glacier velocity, the spatial representation demonstrates whether the variables are spatially related as well, i.e. whether the locations of drainage and velocity are nearby. Therefore, glacier velocity was displayed on a hexagonal grid known as H3 (Sahr, 2011; Brodsky, 2018; Fichtner et al., 2023). The advantage of a hexagonal grid compared to a rectangular grid is the simpler and more symmetric nearest neighbourhood and the better clarity in visualizations. The conversion from the rectangular presentation to the hexagonal presentation is based on resampling using the median value. Only the supraglacial lakes that arose or drained in the observed period of higher glacier velocity were displayed in order to analyse if the draining water could be a control on ice dynamics.

3.8 Pearson correlation

The glacier surface velocities were correlated with the supraglacial lake area, melt extent, run-off index, and meteorological parameters. As each variable had a different temporal resolution, the time series were resampled to daily values. The Pearson correlation coefficient (r) was applied each year for two different periods: (1) the first period was defined as being from the day with the first positive air temperatures at Urdukas until the time of maximum glacier velocity, and (2) the second period was defined as being from the time of maximum glacier velocity until the end of the melt season on 30 September (Berthier and Brun, 2019). The significance of the correlation coefficients with a confidence level of 95 % was estimated using the Student distribution (Obilor and Amadi, 2018).

3.9 Linear regression

To analyse separately the relationship between the data sets and the dependency of all variables for both periods, we applied a linear regression based on a least-square robust adjustment using R² and the regression coefficient (slope) to show the dependency of the variables. The significance of any observed trends was estimated using the Mann–Kendall test with a confidence level of 95 %.

4 Results

First, we present the results of the spatio-temporal relationship, followed by the results of the Pearson correlation, the linear regression, and the estimated temporal delay between melt and surface velocity.

Figure 2 Temporal relationship of the supraglacial lakes, run-off index, surface velocity, air temperature, precipitation, and melt for 1 October 2016 to 30 September 2022 (a–f). All images have the same legend shown in (a). The label of the left y axis is displayed on the outer-left side, and that of the right y axis is displayed on the outer-right side. The abbreviation U stands for Urdukas, Y stands for Yermanendu, G stands for Gore, and C stands for Concordia (surface velocity) whereas W stands for western sector, C stands for central sector, and E stands for eastern sector (melt).

Table 1 Glacier and climate variables with their relevant parameters for October 2016 to September 2022. The maximum area refers to the area at a specific point in time, and the total summer area refers to the summarized area during summer. The abbreviation W stands for western sector, C stands for central sector, and E stands for eastern sector.

4.1 Temporal relationship

The temporal relationship of the supraglacial lakes, run-off index, surface velocity, air temperature, precipitation, and melt for October 2016 to September 2022 is shown in Fig. 2. Table 1 lists the glacier and climate variables and their relevant parameters. The following subsections focus on selected seasonal effects and only consider the influencing parameters.

4.1.1 Warm spring season

During March (earliest: 5 March 2022, latest: 27 March 2019), daily noon air temperatures consistently rose into the positive range for the first time. In mid-April, there was usually a period of higher air temperatures, coinciding with a preceding precipitation event. At this time, the melt line and the zero-degree level (estimated from the air temperatures at Urdukas) were at 4700 m a.s.l.; hence, the precipitation above that altitude fell as snow. Afterward, melt area and air temperatures exhibited a linear relationship, which indicated that higher air temperatures were present in higher altitudes.

At the same time, the snow melted, and supraglacial lakes formed. The melt started between 27 March (earliest: 2019) and 21 April (latest: 2021) and terminated at the end of October. First, the snow melted in the lower altitudes (3600–4300 m a.s.l.) on the debris-covered main branch and debris-free tributary glaciers, and later, it melted only on the debris-free glaciers in the higher altitudes (up to 4700 m a.s.l.).

The debris cover leads to heterogeneous ablation (Nicholson et al., 2018) on the glacier surface with enhanced melting during warm periods on the ice cliffs (Brun et al., 2018; Buri et al., 2021) and in the supraglacial lakes (Miles et al., 2020). The formation of the lakes started in early April to early May and reached the peak of the maximum total area from end May or early June (2020, 2022) through to middle or end July (2017, 2018, 2019) or middle August (2021). It is possible that the surface meltwater could initiate lake formation, though this could not be proven due to the low spatial resolution of the remote sensing data. Spring 2022 was significantly affected by heatwaves (Otto et al., 2023). In summer 2022, the degree-day sum at Urdukas was at 2084 compared to the value range of 1822 (2019) and 1989 (2017), which led to an early supraglacial lake peak (28 May).

4.1.2 Spring and summer speed-up

Between May and September, the melt of snow and ice fluctuated gradually through time between minimum and maximum extent, reaching spatial maxima of 90 km² for the western sector (2018), 84 km² for the central sector (2017), and 190 km² for the eastern sector (2019). In July 2019, the total melt area had a maximum of 353 km² and also covered the debris-free part of the glacier. This corresponds to 64 % of the total area of Baltoro Glacier. A temporal delay from the lower (western) sector to the upper (eastern) sector was hardly discernible, which could probably be explained by an abrupt transition from winter to spring with higher air temperatures at higher altitudes. Differences only existed in the melt area as the eastern sector covers a larger glacier area than the western and central sectors. The first and second melt events occurred in the same period as the spring speed-up and lake evolution. In the area upward from Urdukas, which is characterized by less debris cover, the meltwater accumulates in the supraglacial lakes and drains into the subglacial system or flows laterally into larger streams. In the area downstream of Urdukas, there are no larger supraglacial channels, and the water can drain through moulins and crevasses to the ground.

The winter surface velocities (1 October until the first date with positive air temperatures, usually at the beginning or end of March) ranged, over the years, between 0.05 and 0.37 m d⁻¹ and were followed by a spring speed-up. This event usually lasted from the end of April to the beginning of May and led to a relative increase of 0.1–0.16 m d⁻¹ (30 %–35 %) above Yermanendu. In 2017, 2019, 2020, and 2021, the velocities at Gore exceeded even those at Concordia at this time. The spring speed-up was particularly visible in 2018, affecting all four locations. However, in 2022, the velocities at Urdukas were rarely impacted.

Table 2 Linear regression given as R² and the regression coefficient (in parentheses) of all variables for 2017–2019. Only significant trends according to a 95 % confidence interval based on the Mann–Kendall test are shown.

After the spring speed-up, the flow at Concordia accelerated continuously until the summer peak of 0.46–0.59 m d⁻¹ was reached in June to July. Only in 2022 was a short slowdown between the speed-up and the summer peak recognizable. The high lasted for 11–12 d in 2017, 2019, 2021, and 2022 and for 22 d in 2018 and 2020. However, not all four locations experienced the maximum velocity at the same time: in 2022, the summer peak affected all four locations; in 2021, only the area upglacier of Yermanendu was influenced; in 2020 and 2019, the velocities downglacier of Gore reached their maximum during the spring speed-up; in 2018, the velocities at Concordia and Gore had their high during summer peak, while Yermanendu and Urdukas experienced this during spring speed-up; in 2017, the velocities upward from Yermanendu reached their maximum during summer peak, and the velocities upward from Urdukas experienced this afterward. Additionally, the year 2021 was characterized by a continuous transition between spring speed-up and summer peak. Between June and July, the surface velocities between Concordia and Urdukas spread with a difference of up to 0.46 m d⁻¹, with an increase at Concordia and a decrease at Urdukas. In 2017 and 2021, this phenomenon happened after the summer peak, and in 2018–2020, it happened between speed-up and summer peak.

4.1.3 Autumn speed-up

Supraglacial lakes store large quantities of water and can hence influence the velocity pattern in the event of rapid drainage. The maximum area of the supraglacial lakes ranged between 1.9 km² (2019) and 3.1 km² (2022). The drainage started between early June and mid-July. In 2018 to 2020, it occurred in the period around the summer velocity peak and coincided with a velocity decrease. However, in 2017 and 2021, lake drainage started shortly before the fall speed-up, and in 2022, it started between spring speed-up and summer peak. The supraglacial lake drainage lead to a velocity acceleration from the winter velocities to 0.35–0.5 m d⁻¹ at Concordia in mid-September.

4.1.4 Run-off

The run-off intensity is an indicator of when the water began to flow continuously and thus marks the transition from inefficient to efficient subglacial drainage. It rose from July onward after the maximum of summer speed-up and the maximum of supraglacial lake area and during the melt. The pattern and peak of the run-off index changed every year: 2017 was characterized by a slow increase during May and June, with two distinctive peaks in mid-July (21 July) and early August (5 August); 2018 was characterized by a steady and slow rise, with lows and highs (12 June, 22 July, 10 August); 2019 was characterized by a significant increase in the beginning of July, a constant run-off until the end of July, and a subsequent peak in early August; 2020 was characterized by a slow gain from May until August and two peaks (3 and 25 August), with a continuous run-off index between both dates; 2021 was characterized by a continuous increase until the end of July, succeeded by a high (28 July); and 2022 was characterized by a continuous low run-off index until the end of June and then a sudden rise and, thus, a distinct transition to an efficient subglacial drainage.

4.2 Spatial relationship

The spatial relationship of supraglacial lakes, melt, and surface velocity for the two periods of spring speed-up (13 April–7 May) and summer peak (7 May–18 July) is presented as an example for 2018 in Fig. 3. This year was chosen because the events were clearly recognizable. In the first period (Fig. 3a), the surface velocity between Urdukas and Gore ranged from 0.2 to 0.4 m d⁻¹, and between Gore and Concordia, it ranged from 0.4 to 0.6 m d⁻¹. During this time, the main branch between Gore and Concordia and the neighbouring tributary glaciers up to 4700 m a.s.l. exhibited snow and ice melt covering an area of 237 km² (Fig. 3a, displayed in orange). The supraglacial lakes first formed between the glacier tongue and Yermanendu (increase of 1.1 km²). In the second period (Fig. 3b), the higher velocities of 0.4–0.6 m d⁻¹ expanded up to 5 km below Yermanendu, and the melt expanded upward to 5600 m a.s.l. and covered all neighbouring tributary glaciers, with an area of 327 km² (Fig. 3b, displayed in orange). However, the supraglacial lake formation extended up to Gore. During this period, the supraglacial lakes are characterized by a filling from 7 May to 2 July of 1.3 km², followed by filling (0.33 km²) and simultaneous draining (0.46 km²) from 2 to 18 July.

Figure 3 Spatial relationship of the maximum supraglacial lake area, aggregated melt, and median surface velocity from 13 April to 7 May 2018 (a) and from 7 May to 18 July 2018 (b).

4.3 Pearson correlation

Figure 4 lists the Pearson correlation of all variables for the first period, defined as being from the first day with positive air temperatures to the maximum glacier velocity, and the second period, defined as being from the maximum glacier velocity until the end of the ablation season. Correlation results can be found in the Appendix (Tables A2 and A3). In the two periods, the correlation of the variables of temperature and melt (0.7–0.95), as well as of temperature and supraglacial lake (0.58–0.89), showed a similar pattern with a correlation above 0.58, though with negligible minor changes in the Pearson correlation during both periods (0.1–0.2), except for the decrease of 0.4 for temperatures and supraglacial lakes in the second period of 2020. A high degree of correlation existed between melt and run-off index (0.65–0.93) and melt and supraglacial lakes (0.63–0.94), albeit with outliers in both relationships in the second period of 2020 (decrease of 0.4 and 0.9). The relationship between supraglacial lakes and velocity, as well as between melt and velocity, showed the same pattern, with higher correlations at Concordia and Yermanendu and lower correlations at Urdukas in both periods. However, the correlations of the second period vary significantly (0.5 to −0.8) over the years and generally decrease by 0.4–1.3 in the second period. Particularly noticeable was the variation between positive and negative correlations within the observation period for Urdukas in both parameters. The correlation between velocity and run-off decreased by 1.3 in the second period. However, compared to the two pairs of variables mentioned above, the second period is characterized by a heterogeneous pattern over the years.

Figure 4 Pearson correlation of glacier and climate variables for the first period, defined as being from the first day with positive air temperatures to the maximum glacier velocity (a), and for the second, period defined as being from the maximum glacier velocity until the end of the ablation season (b). The null line is displayed in grey. The abbreviation U stands for Urdukas, Y stands for Yermanendu, G stands for Gore, and C stands for Concordia.

4.4 Linear regression

Tables 2 and 3 list the results of the linear regression with R² and the regression coefficient. In 2018, 2020, and 2022, air temperature and melt, as well as air temperature and supraglacial lakes, had a high R² (0.65–0.83), but the dependency was low (0–7 °C km⁻²). Melt and surface velocity had a high R² for the first period in 2017 in the area upglacier of Urdukas (0.85–0.88), in 2019 in Concordia (0.84), in 2020 (0.78–0.87) and 2022 (0.74–0.9) upglacier of Gore, and in 2021 (0.7–0.8) upglacier of Yermanendu. The regression coefficient lies at 832–1897 km² d m⁻¹. The second period shows a low R². Supraglacial lakes and surface velocity show a high correlation upglacier of Urdukas in the first period of 2017 (0.73–0.87); upglacier of Gore in the first period of 2021 (0.71–0.94); and at Concordia in the first period of 2018 (0.86), 2019 (0.92), 2020 (0.91), and 2022 (0.75) and in the second period of 2019 (0.75). However, the dependency of supraglacial lakes in relation to velocity is minimal (4–11 km² d m⁻¹). R² and the dependency of precipitation in relation to surface velocity are not significant. Melt and run-off index, as well as supraglacial lakes and run-off index, show the pattern, with similar R² (0.77–0.88); however, the dependency of melt (171–370 km² km⁻²) is larger than that of supraglacial lakes (1–4 km² km⁻²). The dependency of the surface velocity in relation to run-off index is not significant (0 m d km⁻²).

Table 3 Linear regression given as R² and the regression coefficient (in parentheses) of all variables for 2020–2022. Only significant trends according to a 95 % confidence interval based on the Mann–Kendall test are shown.

4.5 Temporal delay between melt and velocity

Following the high correlation of melt area and surface velocity during the spring period almost every year and the process-oriented relation between melt onset, inefficient subglacial drainage, and basal sliding, we investigated the temporal delay between the first seasonal melt peak and the first relative velocity maximum at Urdukas, Yermanendu, Gore, and Concordia, which is shown in Fig. 5. The temporal lag between the two parameters is generally longest at Concordia (34±25 d), followed downglacier by Gore (12±19 d) and Yermanendu (12±14 d), and is shortest at Urdukas (10±16 d). During the period from 2017 to 2022, there is an indication that the temporal delay between the first seasonal melt peak and the first relative velocity maximum was decreasing. Thus the amount of meltwater causing glacier velocity variations is reached in a shorter time, which will be more supported by warmer climates in the future. As the measurement covers only 6 years, this observation requires a careful interpretation and needs further investigation.

Figure 5 Temporal delay between first seasonal melt and glacier velocity peak for Urdukas, Yermanendu, Gore, and Concordia. The linear trend – to indicate the tendency only – is displayed as dotted lines.

5 Discussion

In this section, we discuss the relationship between glacier velocity, basal sliding, and hydrological drainage, as well as the impact of melt, supraglacial lakes, and debris cover on this process.

5.1 Transition between efficient and inefficient subglacial drainage

Although each year exhibits a different behaviour, we present in the following a general conceptual model that tries to explain the observed seasonal patterns. The transition of winter to spring was typically observed between the end of March and the end of April when positive air temperatures begin to persist above 0 °C for periods of multiple days up to 4700 m a.s.l. (Fig. 2). The snow on the main branch has usually already melted at Gore at this time. This warm period in spring causes supraglacial lake formation up to 4300 m a.s.l. due to meltwater from snow and ice and the first snowmelt on the tributary glaciers and main branch up to 4700 m a.s.l. (Figs. 3 and 4a). The surface meltwater flows into the ice via crevasses and reaches the cavities at the glacier bed via vertical conduits. As the subglacial drainage system is inefficient, the channels cannot cope with the influx of the first melt, causing a high subglacial water pressure and separation of ice from the bed with a subsequent basal motion (Mair et al., 2001; Macgregor et al., 2005; Burgess et al., 2013; Nanni et al., 2023). The high correlation between melt and velocity confirms this phenomenon (Tables 2 and 3). The spring speed-up with a relative increase of 0.1–0.16 m d⁻¹ (30 %–35 %) affects the complete glacier area between Urdukas and Concordia. In the years 2017, 2019, 2020, and 2021, the velocity at Gore exceeds that at Concordia, suggesting that basal sliding is greatest at this location (Fig. 2). This phenomenon was already observable in 2008–2010 and 2014 (Wendleder et al., 2018).

In summer, persistently warm air temperatures support the formation and expansion of supraglacial lakes by producing large amounts of meltwater from snow and ice (Fig. 2). A detailed overview of the distribution and duration of the supraglacial lakes from 2017 to 2020 is shown in Wendleder et al. (2021a). With the rise in air temperatures, the melt expands up to 5600 m a.s.l. and provides a continuous water influx into the channels and a continuous velocity increase up to 0.59 m d⁻¹ from June to July. A notable indication of basal sliding during the summer months is the substantial increase, ranging from a factor of 2 to 6, in velocity compared to the winter (Mayer et al., 2006). These high summer values were also observed in summer 2004 and 2005 (Quincey et al., 2009) and summer 2008, 2009, and 2010 (Wendleder et al., 2018). The water pressure and the resulting frictional heat leading to a melting of the ice walls cause an expansion of the channels and thus lead to a transition from an inefficient to an efficient subglacial drainage (Schoof, 2010; Sundal et al., 2011; Benn et al., 2017). The increase in the capacity of the drainage system reduces the water pressure and slows the glacier motion at the end of the summer between Urdukas and Concordia down to the winter values of 0.25–0.39 m d⁻¹ (Fig. 2, Nienow et al., 1998; Hewitt, 2013; Vincent and Moreau, 2016).

Afterward, an increase in the time-delayed run-off at the glacier terminus is observable. Only in early July 2022 is an abrupt change from low to high run-off recognizable, indicating a faster transition to an efficient subglacial drainage (Fig. 2). In the other years, however, the run-off linearly increases or fluctuates gradually through time between minimum and maximum extent. Presumably, the higher summer velocities provoke the creation of crevasses and cracks and thus enable the drainage of the supraglacial lakes. The water of the lakes draining after the slowdown lead, in addition to the melt, to a water surplus and hence to an increase in water pressure and basal sliding during autumn (Fig. 2, Tables 2 and 3), but with 0.1–0.2 m d⁻¹ (20 %–30 %) lower magnitude than the summer velocity peak (Hewitt, 2013; Hart et al., 2022; Nanni et al., 2023).

5.2 Impact of melt and supraglacial lakes

In Wendleder et al. (2018), we inferred that warmer spring seasons (April–May) with higher precipitation or melt rates would lead to increased formation of supraglacial lakes and that the discharge of the supraglacial lakes would cause an increased basal sliding. Increased availability of Earth observation data since 2016 has now enabled a continuous mapping of the seasonal evolution of supraglacial lakes (Wendleder et al., 2021a) and melt and thus provides insights into the detailed process.

This study shows that snow and ice melt cover the largest area of 319–353 km² in mid-July to mid-August (Fig. 3). Assuming a snow depth of 1–1.5 m and a density for settled snow of 350 kg m⁻³ at the end of the ablation season (Mayer et al., 2014), this could be roughly estimated to be a snow water equivalent of 122–188 million m³. However, the supraglacial lakes cover a total annual area of 3.6–5.9 km² and a maximal volume of 60 million m³ assuming an average water depth of 10 m (Liu et al., 2015). These estimates show that the water volume of the melt is twice as large as that of the lakes and therefore has a significantly greater influence on the velocity and on the transition to an efficient subglacial drainage.

5.3 Impact of debris cover

Debris-covered glaciers show a different response to global warming because ice melt is reduced by the insulating effect of the debris cover over extensive parts of the ablation area (in the case where debris reaches a critical thickness). Another characteristic is their ability to store large amounts of water in supraglacial lakes, which are suddenly released when crevasses open, often resulting in basal sliding.

The important findings of this study are as follows: (1) the period of several days with positive air temperatures in spring leads to an increase in lake area and to a melt up to 4700 m a.s.l.; (2) the large melting area of snow and ice causes the spring speed-up and the high summer glacier velocities; and (3) Baltoro Glacier also shows an acceleration in autumn, caused by the drained supraglacial lakes. The area affected by the melt and its effects on glacier dynamics are therefore surprisingly large. Despite the insulation due to the debris cover on the main branch of Baltoro Glacier, the ice cliffs and supraglacial lakes act as hotspots that react sensitively to the rise in temperature and lead to an additional ice ablation.

6 Conclusions

To gain a more comprehensive understanding of whether, and how, meltwater and lake drainage events trigger or contribute to basal sliding, we conducted a spatio-temporal analysis. This involved integrating glaciological variables such as surface velocity, supraglacial lake extent, snow and ice melt, and run-off index derived from Earth observation data, along with climate variables including air temperature and precipitation obtained from the HAR data set. The multi-parameter time series was analysed by the Pearson correlation and linear regression. Our study showed that the period of several days with positive air temperatures from the end of March to the end of April led to an abrupt transition from winter to spring and that the ablation period was influenced by continued high air temperatures (Fig. 2). The higher air temperatures affected both snow and ice melt, as well as supraglacial lake formation (Fig. 4a). The snow and ice melt, covering up to 64 % of the total area of Baltoro Glacier (Fig. 3b) and lasting from April to November or December (Fig. 2), had the greatest impact on the basal sliding that led to the spring speed-up of 0.1–0.16 m d⁻¹ (30 %–35 %) in relation to the winter velocities and the summer velocity increase of 0.1–0.2 m d⁻¹ (20 %–50 %) in relation to the spring velocities (Tables A2 and A3). The subsequent transition from inefficient to efficient subglacial drainage lead to a glacier slowdown during August and September. Afterward, the supraglacial lakes with a total area of 3.6–5.9 km² (Fig. 3b) contributed to the fall speed-up, which has a 0.1–0.2 m d⁻¹ (20 %–30 %) lower magnitude than the summer velocity peak. The discharge from snow and ice melt accounts for the largest amount of run-off (Tables A2 and A3).

The year 2022 exhibits notable climatic impacts in the form of a sustained warm period (Table 1), giving rise to a series of consequential phenomena. This includes an increased number of supraglacial lakes, an advancement in the timing of maximal lake area, and a substantial amplification of snow and ice melt (Table 1). These factors collectively provoke two separate peaks in surface velocities during spring and summer (Fig. 2f), driven predominantly by the mechanism of basal sliding. Additionally, this altered dynamic triggers a more rapid transition toward an efficient subglacial drainage system (Fig. 2f).

There exists an indication that the temporal delay between the initial peak of seasonal snow and ice melt and the first relative velocity maximum is decreasing. This observation, covering only 6 years, warrants a careful interpretation, suggesting that the glacier's response is intensifying in relation to the processes of snow and ice melt and velocity. Consequently, this increases the glacier's potential for stronger responsiveness to climatic variations, in particular elevated air temperatures, thereby implying an augmented susceptibility to the effects of climate change.

Appendix A

Table A1 Observation period of glacier surface velocities derived from Sentinel-1 (orbit 27). Indicated are the dates of the first and second acquisitions used to derive the glacier velocity.

Table A2 Pearson correlation (r) of glacier and climate variables for 2017–2019. Only significant correlations according to a 95 % confidence interval are shown.

Table A3 Pearson correlation (r) of glacier and climate variables for 2020–2022. Only significant correlations according to a 95 % confidence interval are shown.