Surge dynamics and lake outbursts of Kyagar Glacier , Karakoram

The recent surge cycle of Kyagar Glacier, in the Chinese Karakoram, caused formation of an ice-dammed lake and subsequent glacial lake outburst floods (GLOFs) exceeding 50 and 40 million m3 in 2015 and 2016, respectively. GLOFs 5 from Kyagar Glacier reached double this size in 2002 and earlier, but the role of glacier surging in GLOF formation was previously unrecognised. We present an integrative analysis of the glacier surge dynamics from 2011 to 2016, assessing surge mechanisms and evaluating the surge cycle im10 pact on GLOFs. Over 80 glacier surface velocity fields were created from TanDEM-X, Sentinel-1A and Landsat satellite data. Changes in ice thickness distribution were revealed by a time series of TanDEM-X elevation models. The analysis shows that during a quiescence phase lasting at least 14 years, 15 ice mass built up in a reservoir area at the top of the glacier tongue and the terminus thinned by up to 100 m, but in the two years preceding the surge onset this pattern reversed. The surge initiated with the onset of the 2014 melt season, and in the following 15 months velocity evolved in a manner con20 sistent with a hydrologically-controlled surge mechanism. Dramatic accelerations coincided with melt seasons, winter deceleration was accompanied by subglacial drainage, and rapid surge termination occurred following the 2015 GLOF. Rapid basal motion during the surge is seemingly controlled 25 by high water pressure, caused by input of surface water into either an inefficient subglacial drainage system or unstable subglacial till. The potential lake volume increased to more than 70 million m3 by late 2016, as a result of over 60 m of thickening at the terminus. Lake formation and the evolution 30 of the ice dam height should be carefully monitored through remote sensing to anticipate large GLOFs in the near future.

2) We have made a new observation about a supraglacial lake on the eastern edge of the glacier, near the confluence, which is well visible with TanDEM-X backscatter images. This lake becomes visible about a month before the GLOF in 2015, was present during the GLOF but had disappeared by the day after. In 2016 this lake again appeared, a bit over a month before the first 2016 GLOF and disappeared after the outburst. We see the presence of this small lake as an expression of high englacial water pressure before the GLOFs, at which time glacier velocity was heightened. The disappearance of the small lake after the GLOFs indicates the establishment of more efficient en/subglacial drainage and corresponds with the deceleration of the glacier following the GLOF. We have added a paragraph about this observation in the lake drainage results section 4.4, p.11, a mention in the discussion section 5.1.2 p.13, and Figure 1 below is to be included in the supplementary figures. We hope that you agree this is a valuable addition to the paper, even if it comes rather late in the submission process. Figure 1: This image is to go in the supplementary figures with the following caption: A series of TanDEM-X backscatter intensity images showing a small lake on the eastern edge of the glacier (indicated by the inset box on the full image) which was present before and absent after the GLOF events in 2015 and 2016. The presence of the supraglacial lake is an expression of high englacial water pressure during the surge and its disappearance after the GLOFs strongly indicates the establishment of more efficient en/subglacial drainage.
Many thanks again, on behalf of the authors, Vanessa Round from Kyagar Glacier reached double this size in 2002 and earlier, but the role of glacier surging in GLOF formation was previously unrecognised. We present an integrative analysis of the glacier surge dynamics from 2011 to 2016, assessing surge mechanisms and evaluating the surge cycle im-10 pact on GLOFs. Over 80 glacier surface velocity fields were created from TanDEM-X, Sentinel-1A and Landsat satellite data. Changes in ice thickness distribution were revealed by a time series of TanDEM-X elevation models. The analysis shows that during a quiescence phase lasting at least 14 years, 15 ice mass built up in a reservoir area at the top of the glacier tongue and the terminus thinned by up to 100 m, but in the two years preceding the surge onset this pattern reversed. The surge initiated with the onset of the 2014 melt season, and in the following 15 months velocity evolved in a manner con-20 sistent with a hydrologically-controlled surge mechanism. Dramatic accelerations coincided with melt seasons, winter deceleration was accompanied by subglacial drainage, and rapid surge termination occurred following the 2015 GLOF. Rapid basal motion during the surge is seemingly controlled 25 by high water pressure, caused by input of surface water into either an inefficient subglacial drainage system or unstable subglacial till. The potential lake volume increased to more than 70 million m 3 by late 2016, as a result of over 60 m of thickening at the terminus. Lake formation and the evolution 30 of the ice dam height should be carefully monitored through remote sensing to anticipate large GLOFs in the near future.

Introduction
Glacier surges are dynamic instabilities affecting about 1% of glaciers worldwide . They consist 35 of periodically alternating long quiescent phases, characterised by years to decades of slow flow, and short active surge phases, characterised by months to years of acceleration and mass transport down the glacier (Meier and Post, 1969). During the active surge phase, the glacier typically 40 experiences dramatic lengthening or thickening at the terminus with potentially hazardous consequences, in particular ice-dammed lake formation (Harrison et al., 2014). While surging glaciers in North America and Svalbard have been investigated in considerable detail, the large concentration of 45 surge-type glaciers existing in the central Asian mountains, including the Karakoram (Copland et al., 2011) are less extensively studied. Improved understanding of surge glacier dynamics in this region can assist anticipation of glacier behaviour and hazard development in the future. The recent un-50 precedented collapse of two surging glaciers in Tibet (GAP-HAZ, 2016) highlights the potentially unexpected nature of glacier instabilities in the region.
Kyagar (Keyajir) Glacier, situated on the northern slopes of the Karakoram Mountains, occasionally causes glacial 55 lake outburst floods (GLOFs) with devastating impacts on downstream communities along the Yarkant River in northwestern China (Zhang, 1992;Hewitt and Liu, 2010;Haemmig et al., 2014). The lake forms when ice at the glacier terminus impounds the river in the Upper Shaksgam Val-60 ley. Owing to the remote location of Kyagar Glacier, about 450 km upstream of the Yarkant floodplain ( Fig. 1), the origin of these floods was poorly understood in the past and they arrived without warning (Zhang, 1992; Hewitt and Liu, ). An automated monitoring station was placed at Kyagar Glacier in 2012 to assess lake formation (Haemmig et al., 2014), at which time there was no lake as the river flowed through subglacial channels at the terminus. From mid-2014, camera images from the station showed dramatic vertical 5 thickening of the glacier terminus followed by lake formation. This rapid thickening indicated a possible glacier surge and increased the potential ice-dammed lake volume.
Although it was already recognised in the 1990s that Kyagar Glacier sometimes dammed the river in the Upper Shaksgam Valley and that there had been periods of advance or thickening in the late 1920s and 1970s (Zhang, 1992) and 1990s (Hewitt and Liu, 2010), the possibility of Kyagar being a surge-type glacier wasn't realised until recently (Gardelle et al., 2013;Haemmig et al., 2014) and no surge of the 15 glacier has ever been documented. An investigation of historic GLOF occurrences from Kyagar Glacier shows that there have been three main periods of flooding in the last 60 years, with peak volumes larger than 130 million m 3 in 1961, 1978and 1999. At least the latter two of these 20 periods coincide with periods of suspected advance or thickening. Recurring GLOFs linked to periods of glacier surging have also been observed for other surging glaciers (e.g. Hoinkes, 1969).
Surging affects both temperate and polythermal glaciers 25 with a variety of geometries and settings (Clarke et al., 1986;Clarke, 1991;Jiskoot et al., 2000) and on vastly different timescales. The general mechanism is as follows: an unstable profile develops during quiescence, as slow velocities fail to redistribute accumulated mass from the upper to the lower 30 part of the glacier, causing steepening of the surface and increasing basal shear stresses, to a point at which surging occurs through dramatically accelerated basal sliding (Raymond, 1987). The proposed mechanisms by which the accelerated basal sliding occurs are various and not completely 35 understood, particularly because the subglacial environment is so difficult to observe. A switch in basal thermal conditions has been identified as a surge mechanism for some polythermal glaciers, with surging occurring when cold basal conditions switch to temperate (Clarke et al., 1984;Murray et al., 40 2000;Fowler et al., 2001). On the other hand, for temperate glaciers and many polythermal glaciers that are already temperate at the base (Sevestre et al., 2015), surging has been explained by a hydrological switch mechanism, by which a surge occurs when the subglacial drainage system becomes 45 inefficient, raising subglacial water pressure and facilitating rapid sliding (Kamb et al., 1985;Björnsson, 1998). Rapid deformation within subglacial till, in response to disturbance of the hydrological system within the till and increased effective water pressure, has also been proposed as an important 50 possible surge mechanism, and is the largest uncertainty in surge understanding (Boulton and Jones, 1979;Truffer et al., 2000;Harrison and Post, 2003). In all cases, a number of positive feedback mechanisms may enhance basal motion during a surge, for instance feedbacks between deformation, 55 frictional heating and subglacial water pressure (Weertman, 1969;Clarke et al., 1984;Sevestre et al., 2015).  Zhang (1992) and Chen et al. (2010). Volumes from 2006 -2016 are estimated from lake extent on satellite images.
Glacier surging in the Karakoram region has mainly been studied by satellite remote sensing, due to the difficulty of field access. Glaciers have been classified using visible morphological features (Barrand and Murray, 2006;Copland et al., 2011) and surge dynamics have been observed through 5 surface velocities (Quincey et al., 2011(Quincey et al., , 2015Mayer et al., 2011), producing some contradicting conclusions on surge mechanisms. Quincey et al. (2011) interpreted a lack of seasonal control on surge initiation as an indication of thermally controlled surges, whereas Mayer et al. (2011) pro-10 posed a hydrological switch mechanism for North Gasherbrum Glacier. Quincey et al. (2015) concluded that Karakoram glacier surging must be quite heterogeneous with a spectrum of surge mechanisms at play, having observed surges exhibiting a surge-front like down-glacier acceleration as 15 well as surges showing simultaneous glacier-wide acceleration. In the nearby West Kunlun Shan, two glacier surges showed a clear seasonal modulation of velocities during the active phase with winter velocities up to 200% higher (Yasuda and Furuya, 2015). The main limitations of these stud-20 ies were data gaps meaning that various stages of the surge development, such as surge initiation, weren't observed, and changes to ice mass distribution during surging also weren't investigated. For such an investigation, digital elevation models (DEMs) from before and after the surge would be re-25 quired.
In this study the combination of optical and synthetic aperture radar (SAR) satellite data reveals the lead up, the onset and termination of the surge, as well as velocity modulations in relation to the seasonal cycle during the surge phase. A 30 DEM time series exposes the ice mass distribution changes caused by the surge and allows us to examine the mass buildup which ultimately drives the surge. Our analysis of the most up-to-date available satellite tools provides a synthesis of the dynamics of a Karakoram glacier in unprecedented de-35 tail, showing the relationships between surging and external factors such as seasonal melt cycles and lake drainage events.
In addition, we assess the impact of surging on the GLOF hazard posed by Kyagar Glacier in the recent past and into the future. GLOF hazard is largely determined by the lake 40 volume and its drainage rate (Björnsson, 2010). The presented time series of glacier DEMs allows for the estimation of the potential lake volume and projection of potential GLOF volumes in the near future, and high-resolution satellite images reveal the drainage mechanism. 45 2 Study site and data

Study site
Kyagar Glacier is a polythermal glacier spanning from 4800 to over 7000 m a.s.l., consisting of three upper glacier tributaries 6-10 km in length which converge to form an 8 km 50 long glacier tongue, approximately 1.5 km wide (Fig. 3). The total glacier area is 94 km 2 (Randolph Glacier Inventory Version 5.0, 2015) and the average surface slope is approximately 2 • over the tongue and 4.5 • -20 • over the branches above the confluence. The surface of the glacier tongue is 55 characterised by ice pinnacles (Fig. 4) up to 40 m high and as narrow as 10 m, indicating cold ice and low shear deformation (Haemmig et al., 2014).
The tongue of Kyagar Glacier is most likely carved into the brown/black shales and cherty limestones of the 3 km 60 thick Perminan-Jurassic Shaksgam sedimentary formation, while the mountain range forming the southern margin of the glacier catchment consists of the Aghil formation limestone and perhaps dolomite (Desio et al., 1991). The Shaksgam Valley follows the Shaksgam fault which passes under 65 the terminus of Kyagar Glacier (Searle and Phillips, 2007).
Fieldwork at Kyagar Glacier is limited because of its remoteness and political restriction of access. Since a Sino-Swiss expedition in 2012 (Haemmig et al., 2014), in situ observations became available from an automated monitoring 70 station about 500 m upstream of the Kyagar Glacier terminus (Fig. 3), which operated from 7 September 2012 until being drowned by the growing lake on 29 June 2015. The northern Karakoram is largely influenced by westerly weather patterns and snow accumulation mainly in winter, while rainfall 75 (at lower altitudes) peaks between May and September (Kapnick et al., 2014). Balanced or slightly positive mass balances for Karakorum glaciers between 1999-2011 (Gardelle et al., 2013) contradict global trends of decreasing glacier mass balance in line with global warming, but may be explained 80 by regional increases in winter precipitation (Kapnick et al., 2014).

Data
In situ data from the automated observation station included daily camera images of the glacier terminus, showing the up-85 stream face of the ice dam (Fig. 4). Meteorological variables included air temperature and precipitation amount and type, among others, recorded at hourly intervals until the station became submerged on 29 June 2015. Further meteorological data and river water level measurements were available from monitoring stations on the Yarkant River located at Cha Hekou and Kuluklangon, 190 and 500 km downstream from 5 Kyagar Glacier, respectively (Fig. 1).
Three different satellite systems were used to determine surface velocities from the end of 2011 to mid-2016 ::: the ::: end :: of :::: 2016: the SAR systems Sentinel-1A and TanDEM-X and the Landsat-8 optical system. TanDEM-X is a formation of two tandem satellites, TanDEM-X and TerraSAR-X, data from both of which are used for all velocity and elevation analyses. In addition, Sentinel-2 optical images were used for visual assessment of lake formation but not for velocity analysis. Acquisition details of the three main satellite systems are presented in Table 1 (for a complete list of acquisitions see supplementary material).

Image co-registration
All satellite scenes were co-registered to a common mas-20 ter scene to allow accurate comparison of images from the same orbit. The master scene for Sentinel-1A and Landsat-8 images was the first available image from each orbit. For TanDEM-X, another co-registration algorithm was used where the master was updated progressively as the average 25 of all previously co-registered scenes in order to temporally smooth out moving features or strongly changing patterns such as snowmelt. The scenes used for image co-registration covered an area of approximately 30 × 50 km 2 extending north from Kyagar Glacier, over which local offsets were 30 calculated for patches of 512 × 512 pixel 2 . To remove offsets resulting from patches covering moving glaciers, a planar function was fitted to the offset-fields and large outliers  were removed, before again fitting a planar function to the filtered offset fields. The scenes were then resampled according to the final fitted function, resulting in a stack of images with sub-pixel co-registration accuracy.

Glacier surface velocity 5
Glacier surface velocities were determined using offset tracking, through which the ground offset between corresponding patches of co-registered repeat-pass satellite image pairs is computed (Strozzi et al., 2002;Luckman et al., 2007). Intensity cross-correlation was applied to paired patches from 10 the SAR images, while for Landsat optical data phase crosscorrelation was used to better deal with variable illumination conditions (Zitova and Flusser, 2003). The resulting offset field covering the glacier and its surroundings was then converted to surface velocity by dividing by the elapsed time 15 between the paired images and scaling by ground range resolution. Longitudinal velocity profiles were determined along a manually determined central glacier flowline (as shown in Fig. 3) in the velocity offset patch coordinates. The patch size and spacing are presented in Table 2. Patch 20 sizes were selected to optimise the superior correlation ability of larger patches with the superior spatial resolution of smaller patches. Larger patches were required for the SAR systems, despite their finer resolution, to compensate for radar speckle. Velocity fields were filtered to remove offsets 25 calculated with low correlation quality, as determined by the height of the correlation-function peak over the noise. Offsets with high divergence from neighbouring values and outliers with velocities 50% larger than the maximum offset over the glacier were also removed.

30
The accuracy of the offset tracking procedure was assessed by calculating the patch offsets over a 1 × 2 km 2 area of stable ground next to the glacier terminus. Since no offsets are expected over stable ground, offsets represent local inaccuracies in the co-registration of images caused by slight changes 35 in imaging geometry and, hence, scene projection, as well as the inherent inaccuracy in the sub-pixel determination of the correlation-function peak. The root-mean-square error for the offsets over stable ground was 0.08 pixels or less for almost all the used image pairs from the three satellite sys-40 tems, similar to the 0.05 pixel error estimated by Strozzi et al. (2002). The velocity difference from assuming a horizontal surface for velocity calculations was only about 0.06% over the 2 • sloping glacier tongue and 0.4% over the 5 • slope just above the confluence. Due to the side-looking radar imag-45 ing geometry, steep slopes in the range direction result in a different pixel spacing hence biased velocities. However as the glacier is not very steep and flows predominately in the azimuth direction, this is not a problem. 50 Digital elevation models were derived using data from the TanDEM-X satellite formation (Krieger et al., 2007(Krieger et al., , 2013 using single-pass SAR interferometry. SAR interferometry allows accurate DEM generation if the absolute interferomet-ric phase can be successfully determined from the wrapped interferometric phase measured between 0 and 2π. Determination of the absolute phase requires phase unwrapping algorithms to be applied to the interferogram (e.g. Goldstein et al., 1988;Zebker and Yanping, 1998). The phase unwrapping can be simplified by first subtracting a synthetic interferogram, based on a reference DEM, and adding it back after unwrapping is completed (e.g. Dehecq et al., 2015). This can help minimise phase-wrapping errors which can easily be recognised in the interferograms when an accurate reference 10 DEM is used. If the phase difference to the measured data does not exceed 2π, phase unwrapping can even be avoided entirely.

Digital elevation models
The phase gradient and hence the DEM accuracy depends on the perpendicular interferometric baseline B ⊥ , the com-15 ponent of the distance between the two SAR satellites which is perpendicular to both the line-of-sight and the flight direction. Large baselines provide a better height accuracy with phase cycles of 2π corresponding to smaller height of ambiguity (HoA, see p. 3320 and Eq. 37 in Krieger et al., 2007), 20 but on the other hand large baselines are more prone to phase unwrapping errors and signal decorrelation due to scattering volumes (Zebker and Villasenor, 1992) and noise contained in the reference DEM.
The reference DEM used for this study was composed giving large HAs :::: HoAs : of 250-400 m. For noise reduction, an adaptive filter was applied to the interferograms (Goldstein and Werner, 1998). After phase unwrapping and conversion to height, the height corrections were averaged and added to the SRTM DEM to form the reference DEM which 35 was downsampled to a resolution of 8 × 8 m 2 . DEMs for each acquisition from orbit 75D were created by converting the phase difference against the reference DEM into a height change ∆h, which was then added to the reference to obtain an absolute DEM for each acquisition date. The ex-40 tremely rough glacier surface topography, with ice pinnacles up to 40 m high and 20-40 m apart (estimated from shadow lengths and the observations from Haemmig et al., 2014), caused strong decorrelation and phase wraps within the coherence window of 15 × 15 m 2 for large baselines, meaning 45 that DEMs could not be created over the glacier tongue with baselines B ⊥ >200 m (HoAs below 20 m).
The generated DEMs contain errors from processing uncertainties as well as from microwave penetration into snow. Processing uncertainties include phase noise due to low cor-50 relation in the interferograms, global offsets due to geometric errors, and errors of the SRTM DEM. Errors due to phase noise were estimated from the differences between the eight DEMs used for the reference DEM. The standard deviation was below 4 m. The SRTM DEM is specified with an abso- 55 lute vertical accuracy of about 10 m (Farr et al., 2007), but for comparison of DEMs systematic vertical shifts or tilts were corrected for by referencing DEMs to a common reference height of flat terrain near the tongue of Kyagar Glacier. The remaining relative error between different DEMs was esti-60 mated from four flat valley planes and resulted in a maximum height error of ± 1 m (standard deviation 0.65 m).
The error due to microwave penetration into dry snow can reach up to 6 m (Dehecq et al., 2015) for a microwave frequency of 9.65 GHz (TanDEM-X), but penetration is negli-65 gible over wet snow (more than 1% volumetric water content, Leinss et al., 2015, Fig. 5). Microwave penetration leads to potential underestimation of the actual surface height over dry snow and ice surfaces. For Kyagar Glacier, penetration depths of up to two meters have been estimated by distin-70 guishing between dry and wet snow conditions based on backscatter intensity (Nagler and Rott, 1998;Small, 2012;Nagler et al., 2016), and determining the apparent height difference between DEMs from wet versus dry conditions. Over the tongue of Kyagar Glacier, the backscatter intensity 75 changed little between seasons (<5 dB), because infrequent snowfall means that the bare ice surface roughness dominates the backscatter signal from the tongue. Penetration is therefore expected to be negligible over the glacier tongue. In contrast, large seasonal changes in backscatter intensity indicate 80 changing water content and thus varying penetration depths over the accumulation basin. Backscatter decreased by more than 10 dB at the onset of snowmelt in 2015 over the accumulation areas, and an apparent surface height increase of less than 2 m was calculated between two large baseline interfer-85 ograms from before snowmelt (2015-06-02) and at the onset of snowmelt (2015-06-13). This indicates a TanDEM-X penetration depth of 2 m or less in dry snow conditions over the upper glacier. The relatively small penetration depths in the accumulation area might be a result of strongly-scattering 90 high-density firn with ice inclusions, formed by refreezing after strong melt events extending to over 6000 m a.s.l. in August ( Supplementary Fig. 1), a phenomenon also observed by Dehecq et al. (2015). The backscatter intensity changes are shown in Supplementary Figures 1 and 2. 95 The penetration error in the SRTM DEM should be slightly larger than for the TanDEM-X DEMs, as the SRTM DEM was acquired with a C-band radar with 5.3 GHz (Farr et al., 2007) during winter (Feb. 2000). For C-band radars, expected penetration depths are 1-2 m into exposed ice (Rig-100 not et al., 2001) and 5-10 m into dry snow (Rignot et al., 2001;Fischer et al., 2016). However, because in the accumulation area the penetration for X-band is <2 m, we estimate a penetration of <4 m in C-band (cf. Fig. 9 in Fischer et al., 2016), and over the glacier tongue 1-2 m. 105 In summary, systematic shifts are removed when comparing DEMs but differences in penetration must be considered in particular when comparing the SRTM to the TanDEM-X DEMs or when comparing DEMs from different seasons. Over the glacier tongue, penetration errors are <2 meters and 110 over the accumulation area they are estimated to be <4 meters.

Calculation of positive degree days
Positive degree days (PDD) at the glacier terminus were calculated as a proxy for potential melting. Positive air tem-5 perature measurements were summed with each measurement weighted by the fraction of a day which it represented (Vaughan, 2006), such that one hourly measurement of 6 • would contribute 0.25 PDD. The hourly air temperature data from the station at Kyagar Glacier were used for comput-

Lake volume estimation
Lake volumes were calculated using the DEM of the empty lake basin from TanDEM-X data acquired on 18 Aug. 2016, 20 together with the lake extent and thus lake surface altitude from optical (Landsat or Sentinel-2) or SAR backscatter images (Sentinel-1A and TanDEM-X). In addition, the initial lake formation during the winter of 2014/15 was observed by the in situ camera, as the small initial volumes were not 25 seen on the satellite images but were important for assessing possible subglacial drainage. Potential lake volumes were estimated by calculating the lake volume which would result if the lake basin would be filled to the 90% of the ice dam height, as determined from the DEM of the glacier terminus. In the 2.5 years before surge onset, a gradual but clear acceleration occurred, greatest over the middle of the glacier 45 tongue (between km 3 and km 6) with an increase in velocity from 0.1 m d -1 in winter 2011/12 to over 0.4 m d -1 in winter 2013/14 (Fig. 5). The location of the maximum veloc- Figure 5. Longitudinal surface velocity profiles showing (from bottom to top) the pre-surge acceleration that occurred in the 2.5 years before the main surge onset. The profiles, derived from TanDEM-X data, follow the longitudinal path from Fig. 3. Gaps above km 9 indicate failed velocity calculation owing to the poor surface contrast providing no clear correlation. The labels state the time period over which each velocity calculation was averaged, in this case ranging from 2-13 months. ity moved from above the confluence at km 10 at the end of 2011 to over the glacier tongue at km 5 in 2013/2014. Apart 50 from this early shift, the spatial pattern of acceleration over the glacier tongue was quite uniform with no evidence of a surge front moving down the glacier, as observed for some other Karakoram glaciers (Mayer et al., 2011;Quincey et al., 2015). The presence of seasonal modulation could not be as-55 sessed due to the coarse temporal resolution of the six presurge acquisitions, but it can be seen that acceleration continued over the winter immediately before surge initiation (Fig. 5, fastest velocity profile, from Oct. 2013-Feb. 2014). This gradual pre-surge acceleration may indeed have already been 60 under way prior to 2011, with acceleration between annual velocities from 2004/2005 to 2010/2011 based on Landsat velocity analysis by Heid and Kääb (2012).
The pre-surge acceleration appears insignificant in comparison to the main surge phase, which started at the end of 65 April 2014. Rapid acceleration first became evident between 28 Apr. and 30 May 2014 ( Fig. 6a-b), with a doubling of maximum velocity from 0.5 m d -1 to over 1 m d -1 within 32 days. Velocities continued increasing steadily (Fig. 6c) to a peak of almost 2.5 m d -1 (Fig. 6d) between the 19 Sep. and 70 5 Oct. 2014. The maximum instantaneous velocity is likely to have been higher than the calculated values which are averages over 16-day periods. The surge caused a six-fold acceleration in the five months following May 2014, and more than a 20-fold acceleration since 2011/12. 75 After the surge peak in Sep. 2014, there was a slight deceleration which continued during winter (Fig. 6e) until maximum velocities had dropped to about 1.2 m d -1 in April 2015. This was followed by a new phase of acceleration through  May-July 2015 to almost 2 m d -1 in late July, slightly slower than the peak velocity in summer 2014. This acceleration came to an abrupt halt between 27 July and 7 Aug. 2015, causing the most rapid change observed with a halving of velocities over the tongue within 22 days (Fig. 7a-c). This 5 abrupt slow-down was aligned with the lake drainage on 27 July, as indicated by the arrow in Fig. 8.
Deceleration continued over autumn 2015 and winter of 2015/16, and velocities almost returned to pre-surge levels with a maximum less than 0.5 m d -1 in March 2016. There was a slight acceleration after April 2016 but velocities were still significantly below the previous two summers, remaining below 1 m d -1 . ::: and ::::::::: decreasing ::: to ::::: below ::::::::: 0.5 m d -1 files along the glacier, shows that the surge mainly affected the tongue of the glacier, between km 1 and km 8, while above the confluence (> km 8) the effect of the surge was small.

Glacier surface elevation 20
Four DEMs based on TanDEM-X data acquired before the surge (2012)(2013)(2014) and eight DEMs from after the main part of the surge (Oct.-Dec. 2015) were compared with each other and with the SRTM DEM from 2000, to reveal the dramatic changes in glacier surface elevation and, hence, ice mass dis-25 tribution over Kyagar Glacier caused by the surge. Maps of elevation change over the glacier during the quiescence and surge periods are shown in Fig. 9 and 10, respectively. Longitudinal profiles of surface elevation are shown in Fig. 11a and the elevation change rates in Fig. 11b. Before the surge, between Feb. 2000 and Nov. 2012, the surface elevation decreased gradually over most of the glacier tongue at a rate of 5 m a -1 ( Fig. 9 and Fig. 11b), resulting in an elevation loss of over 60 meters at the glacier terminus. At the same time, elevation increased over the 5 western branch by up to 30 meters just above the confluence, while over the eastern branches the surface elevation increased more moderately with a maximum gain of 10 meters (Fig. 9). This observed pattern is typical of a surging glacier in the quiescence phase, with downwasting over the glacier tongue and ice build up in a reservoir area, which for Kyagar Glacier forms just above the confluence on the western branch.
Between 2012 and 2014, in the two years preceding the surge, there was already a slight reversal of the quiescence 15 pattern seen in the previous 12 years, with minor elevation loss just above the confluence and mass gain over the tongue (Fig. 11), indicating mass transport down the glacier from the reservoir. During the surge in 2014/15, this mass transport from the reservoir area intensified dramatically with ice sur-20 face elevation increasing at a rate of almost 40 m a -1 over the lowest parts of the glacier tongue, causing thickening in excess of 60 meters at the terminus since Feb. 2014. At the same    Fig. 9 and 10, and the black arrow indicates the location of the surface hollow remaining after the surge.
There were significant adjustments to surface slope throughout the course of the surge, particularly over the glacier tongue. In 2000, the average slope over the first eight kilometres of the glacier was about 1.4 • and by 2012/2014 it had increased to about 2.3 • . By the end of 2015, after the surge, slope had decreased again to about 1.6 • .
The post-surge glacier surface in Oct. 2015 showed the presence of a surface hollow, approximately 12 m deep and up to 1 km wide, at the very beginning of the western branch 15 above the confluence just before the slope significantly steepens (Fig. 11a, indicated by arrow). Such a depression is an unusual feature but could have formed as a consequence of the surge transporting ice from the reservoir area faster than the rate of replacement from above, owing to the observed 20 flow disparity between the tongue and the glacier branches.

Mass balance and equilibrium line altitude
Although not directly related to the surge characterisation, we provide a geodetic mass balance estimate for Kyagar Glacier between 2000 and 2015. The average volume 25 difference between the SRTM DEM and eight TanDEM-X DEMs from between Oct.-Dec. 2015 was calculated and converted to mass change assuming an ice density of 850 ± 60 kg m -3 (Huss, 2013). The mass balance was found to be -0.24 ±0.22 m w.e a -1 . For the uncertainty, the radar 30 penetration difference between the SRTM and the TanDEM-X DEMs dominates and was estimated to be a conservative 3 m systematic error over the whole glacier. As the penetration for the SRTM C-band microwaves is deeper than the TanDEM-X X-band, our calculation may slightly un-35 derestimate mass loss. On the other hand, the area used for calculation (61 km 2 ) missed some of the steepest portions of the accumulation area due to lack of interferometric coherence affecting DEM creation, possibly leading to an under-representation of the accumulation area and exag-40 gerated mass loss. For comparison, Gardelle et al. (2013) reported an average mass balance of +0.11 ±0.14 m w.e a -1 for glaciers in the east Karakoram region between 2000 and 2008.
The equilibrium line altitude (ELA) estimated from the lo-45 cation of the snow line at the end of the ablation period observed from Landsat and TanDEM-X images, was 5350±80, 5400±80 and 5510±80 m a.s.l. over the western, middle and eastern branches, respectively.

Meteorological observations 50
Temperatures remained below 0 • C between mid-October and late April according to data from the meteorological station at the glacier terminus (at 4800 m a.s.l.). The warmest months, July and August, experienced average daily maximum temperatures of 4-7 • C and monthly PDDs exceeding 55 150 at the glacier terminus. By taking into account the glacier surface elevation and an lapse rate of about -0.006 • C m -1 , it can be inferred that over the whole glacier tongue, PDDs are positive between May and October, whilst over the bulk of the accumulation area (about 900 m above the terminus) 60 melt potential was only significant from June to August. Evidence of high-altitude melt is also seen in the TanDEM-X backscatter images from August 2015 ( Supplementary Fig.  1). Annual PDDs at the glacier terminus were 647 • C, 481 • C, 552 • C and 528 • C in 2013, 2014, 2015 and 2016, respec-65 tively. The melt rate at the terminus is estimated to be around 5 m a -1 , according to the terminus surface elevation decrease during quiescence (Fig. 11) and the melt rate of icebergs left in the empty lake basin after lake drainage in 2009 (Haemmig et al., 2014). Combining this melt rate and the annual 70 average of 552 PDDs gives a realistic degree day factor of about 9 mm w.e. • C -1 d -1 .

Lake formation and drainage
Images from the monitoring station at Kyagar Glacier showed that a lake initially began forming in the river basin 75 upstream of the glacier terminus in early December 2014. During January and February 2015 the lake appeared to fill Figure 12. Radar backscatter images of the glacier terminus showing (a) the lake 11 days before drainage, (b) just after the start of drainage and (c) after the lake drainage. Lake drainage clearly occurred through subglacial channels, rather than through dam collapse or over topping. Images from TanDEM-X data provided by DLR.
faster, before remaining at a constant size (still less than 1 million m 3 ) during March. In April 2015 the lake size increased again in line with the onset of spring melting, continuing more rapidly during the summer until reaching an estimated volume of 53 million m 3 before draining through 5 subglacial channels on the 27 July 2015, as observed by TanDEM-X acquisitions (Fig. 12).
Following the drainage in July, a new lake started forming in September 2015 and remained at a volume of approximately 1.5 million m 3 between Oct. and Dec. 2015.

Discussion
Based on the results, we discuss possible surge mechanisms for the observed behaviour of the glacier before and during 35 the main surge phase, and rule out mechanisms which contradict the observed behaviour. The effect of the surge cycle on the GLOF hazard posed by Kyagar Glacier in the past and future is assessed to provide an outlook for its hazard potential. The observed pre-surge acceleration could have arisen through increased internal ice deformation and/or increased basal sliding, both of which may be expected following 45 the steepening of the glacier tongue between 2000 and 2012 (Fig. 11a). The contribution of internal ice deformation u d to surface flow can be estimated with the parallelsided slab assumption with the plain strain approximation (Greve and Blatter, 2009), as where the strain rate factor A = 2.4 × 10 −24 s −1 Pa −3 (for temperate ice, a conservative estimate), ice density ρ = 900 kg m −3 , Glen's exponent n = 3, and gravitational acceleration g = 9.8 m s −2 , leaving the key variables surface 55 slope, α, and ice thickness H (Cuffey and Paterson, 2010). Assuming a constant glacier thickness of 250 m, an estimation on the high side according to the glacier bed profile presented by Haemmig et al. (2014), the 1.4 • surface slope over the glacier tongue in 2000 would result in a surface velocity 60 of 4 mm d -1 . The increased slope of 2.3 • in 2012 would give deformation velocities of around 18 mm d -1 , approximately one order of magnitude lower than the observed 0.1 m d -1 between 2011-2012 (Fig. 5). Hence, it seems that basal motion significantly contributed to flow of the glacier tongue al-65 ready prior to the surge, indicating that the base of the glacier tongue was already temperate and contradicting the thermal mechanism in which a switch from cold to temperate base causes surge onset. Conditions are different above the confluence where the surface slope of around 4.5 • in 2012 could 70 cause surface velocities on the order of 0.1 m d -1 through internal deformation alone, in the same order of magnitude as observed velocities. Pre-surge velocities above the confluence could therefore feasibly occur in a cold-based situation through internal deformation without the contribution of 75 basal motion. However, basal motion upstream of the confluence is not ruled out with this simple calculation. The effect of increased surface slope on basal shear stress, τ b , during the quiescence can be estimated from if the glacier base is assumed to mirror the surface slope (Cuffey and Paterson, 2010). If all variables except the slope 5 are considered constant, the increase from 1.4 • to 2.3 • between 2000 and 2012 over the glacier tongue would have caused a 64% basal shear stress increase, from about 54 to 88 kPa. The thickness increase over the reservoir area would further increase τ b at the upper part of the glacier tongue. 10 Given the potential sensitivity of the subglacial hydrological drainage system efficiency to basal stress (Eisen et al., 2005), the conditions at the end of the quiescence could favour a switch to an inefficient drainage system. However, the slope began decreasing between 2012 and 15 2014 (Fig. 11) while velocity continued to increase. This contradicts the idea that increasing slope alone could have driven the acceleration. Positive feedback mechanisms triggered by increasing basal stress must therefore play a role in the continued acceleration during the late quiescence, and 20 ultimately in bringing the glacier into a critical state before surge initiation. These could include increased frictional heating enhancing melt water production and water pressure at the glacier base (Dunse et al., 2015;Weertman, 1969), or increased basal deformation closing subglacial drainage 25 channels, thus trapping water and increasing water pressure (Clarke et al., 1984;Kamb et al., 1985). Processes within the subglacial till such as a positive feedback between till deformation, consolidation and water pressure (Boulton and Zatsepin, 2006) could also play a role.

30
The continuation of acceleration during winter 2013/2014 (Fig. 13), rather than another slow-down as observed by Haemmig et al. (2014) during the previous winter, may indicate the presence of an inefficient subglacial drainage system. Such winter acceleration was observed prior to the 1982/83 35 surge of Variegated Glacier and was attributed to the establishment of an inefficient linked cavity drainage system with higher water pressure, in part due to low water flux allowing drainage channels to close (Kamb et al., 1985). The rapid spring acceleration observed in May 2014 can be explained by the input of surface meltwater increasing water pressure in an inefficient subglacial drainage system, which was suspected to be present during the preceding winter. Continued acceleration through the summer could re-45 flect increasing subglacial water pressure as meltwater input continued. The deceleration after reaching maximum velocities in Oct. 2014 indicates decreasing subglacial water pressure, perhaps through the gradual evolution of a subglacial drainage system towards the end of summer followed by de-50 clining melt water input and possible subglacial drainage. Evidence for the drainage of en-or sub-glacially stored water during wintertime comes from the observed lake formation starting in December 2014, at a time when temperatures consistently well below zero exclude surface water sources. The 55 lake growth and, hence, the drainage of subglacial water, appeared to end in January 2015. This indicates that most subglacial water was already drained or that subglacial drainage channels closed towards the end of the winter. Closing of the subglacial channels would again put the subglacial system 60 into a state very sensitive to surface water input and allow summer-onset acceleration again in 2015.

Summary of surge mechanism
The various phases of the surge were facilitated by a basal motion mechanism very sensitive to subglacial water pressure, controlled by meltwater input in summer, reduced input and perhaps drainage of most of the subglacial water in early 45 winter, and rapid subglacial drainage during the GLOF in summer 2015. It seems that the surge is well explained by the presence of an inefficient basal drainage system facilitating high subglacial water pressure, corresponding to the mechanism proposed by Kamb et al. (1985). However, the sea-50 sonality observed at Kyagar Glacier is different to the often cited winter initiation associated with closure of subglacial channels in the hydrological switch mechanism (Eisen et al., 2005;Kamb et al., 1985). In the case of Kyagar Glacier, development of an inefficient drainage system in winter does 55 not necessarily facilitate increased subglacial water pressure until the beginning of the melt season, due a lack of liquid water in winter. Surge initiation in winter should not be considered a precondition of hydrologically controlled surges (see eg. Jiskoot and Low (2011)). 60 We note that the idea of surge initiation through formation of an inefficient drainage system as discussed in Section 5.1.2, could be replaced by that of a layer of subglacial till in which increased water pressure reduces till strength to a deformation threshold (Boulton and Jones, 1979;Cuffey and 65 Paterson, 2010). It is likely that the tongue of Kyagar Glacier is underlain by a permeable till, owing to fine-grained sedimentary rock on which the glacier tongue lies (Desio et al., 1991;Searle and Phillips, 2007). We can not speculate further on the exact nature of the subglacial drainage system as 70 there is no field evidence, but conclude that Kyagar Glacier is a system very sensitive to water in-and outputs during the surge, rather than being purely internally regulated.
The spatial pattern of acceleration and elevation change over Kyagar Glacier provides further information about the 75 nature of the surge, in particular that it was the tongue of the glacier which primarily underwent surging, evidenced by the velocity increase (Fig. 8) and the steepening of the profile over the glacier tongue during quiescence (Fig. 11). The build-up of an ice reservoir at the confluence represents the 80 intersection between the surging tongue and the tributaries, which maintain more steady flow and support the recharge of the ice reservoir during quiescence. We note also that looped moraines do not form at Kyagar Glacier because there is no surging of upper tributaries into a non-surging part of the 85 glacier (see supplementary video). Surging confined mainly to the flatter, lower part of the glacier has been observed for a number other Karakoram surges (Mayer et al., 2011;Quincey et al., 2015).

Future outlook for Kyagar Glacier and lake formation
The potential volume of the glacier-dammed lake at Kyagar Glacier depends both on the height of the ice dam at 25 the glacier terminus and whether the subglacial channels through which the lake drains are open or closed. Thickening of over 60 m at the glacier terminus caused potential GLOF volume to increase more than 40-fold since early 2014, to over 70 million m 3 according to the August 2016 DEM of the 30 glacier. GLOF hazard potential is expected to remain high for a number of years as the still slightly elevated tongue velocity continues to transport mass to the terminus area, potentially increasing the height of the ice dam until mass transport to the terminus area falls below the ablation rate. The height of 35 the ice dam is expected to decrease at an estimated rate of 5 m a -1 once mass transport to the terminus ceases, according to the ablation rate of the latest quiescence phase and the estimated melting rate of icebergs left in the empty lake basin after lake drainage in 2009/2010 (Haemmig et al., 2014). Un-40 less the mass balance significantly changes, it would be expected that the next quiescence phase would last until around 2030 based on an estimated 15-20 year return period. The size of future GLOFs depends largely only on the potential lake volume, the volume reached if the lake filled to 45 90% of the ice height. However the actual volume reached may be less it the lake drains before the potential volume is reached. Despite the ice dam being about 5 meters higher in 2016 than 2015, the lake volume was less in 2016 as the outburst occurred at about 85% of the ice dam 50 height. This can be explained by subglacial channels from the 2015 lake outburst providing a weaker, preferential pathway for lake drainage and thus earlier outburst in 2016 fol-lowed by a smaller second outburst. We note that GLOFs > 80 million m 3 have never been followed by a significant 55 drainage event in the next year (Fig. 2), which perhaps indicates that large floods cause formation of subglacial channels large enough to remain open until the following year. Meteorological factors, such as temperature during the GLOF, may also influence the peak flood discharge (Ng et al., 2007). 60 It is important that the height of the ice dam and the lake evolution is monitored through satellite imagery each summer in the years following the surge, to assess imminent GLOF threat and allow affected communities to be better prepared for flood impacts.

Conclusion
Our integrative picture of the recent surge of Kyagar Glacier, built from satellite surface velocity maps, terrestrial station images and DEMs, provides an extraordinary insight into glacial surging in connection with surface hydrology and 70 glacier-dammed lake formation and outburst. After gradual surface velocity increase through the last few years of the quiescence, the glacier entered a state highly sensitive to surface water input. Two dramatic acceleration phases occurred in concurrence with the onset of the surface meltwater pro-75 duction in the seasons of 2014 and 2015, indicating a surge mechanism related to the evolution of the basal hydrological system and associated changes in subglacial water pressure, rather than to an internally controlled switch to temperate basal temperatures. Between the acceleration phases, 80 deceleration was accompanied by drainage of subglacial water, evidenced by the filling of the glacier-dammed lake during the winters of 2014/15 and 2015/16. Lake drainage in July 2015 caused instantaneous deceleration over the whole glacier tongue, indicating that a sudden drainage of the sub-85 glacial water under a large part of the glacier tongue occurred with the lake outburst event.
Surging of Kyagar Glacier is the main driver of icedammed lake formation and GLOFs. The thickening of over 60 m at the glacier terminus during the surge caused poten-90 tial GLOF size to increase almost 40-fold since early 2014, to over 70 million m 3 at the end of summer 2016. The hazard potential of large GLOFs remains high in the next few years, potentially larger than the 2015 and 2016 GLOFs, but the actual magnitude depends on the timing of lake drainage. 95 Remotely sensed data, in particular from TanDEM-X, is essential to the observation of the surge phenomenon and the assessment of hazard formation. The remote sensing of the glacier should be continued to monitor lake formation and the evolution of the ice dam height.