A Glacier Surge of Bivachny Glacier , Pamir Mountains , Observed by a Time Series of High-Resolution Digital Elevation Models and Glacier Velocities

Surge-type glaciers are characterised by relatively short phases of enhanced ice transport and mass redistribution after a comparatively long quiescent phase when the glacier is virtually inactive. This unstable behaviour makes it difficult to assess the influence of climate change on those glaciers. We describe the evolution of the most recent surge of Bivachny Glacier in the Pamir Mountains, Tajikistan between 2011 and 2015 with respect to changes in its topography and dynamics. For the relevant time span, nine digital elevation models were derived from TanDEM-X data; optical satellite data (Landsat 5, 7 and 8, EO-1) as well as synthetic aperture radar data (TerraSAR-X and TanDEM-X) were used to analyse ice flow velocities. The comparison of the topography at the beginning of the surge with the one observed by the Shuttle Radar Topography Mission in 2000 revealed a thickening in the upper part of the ablation area of the glacier and a thinning further down the glacier as is typically observed during the quiescent phase. During the active phase, a surge bulge measuring up to around 80 m developed and travelled downstream for a distance of 13 km with a mean velocity of 4400 m year−1. Ice flow velocities increased from below 90 m year−1 duringthe quiescent phase in 2000 to up to 3400 m year−1 in spring 2014. After reaching the confluence with Fedchenko Glacier, the surge slowed down until it completely terminated in 2015. The observed seasonality of the glacier velocities with a regular speed-up during the onset of the melt period suggests a hydrological control of the surge related to the effectiveness of the subglacial drainage system.

With the availability of temporally high-resolution remote sensing data, the chance to observe short-term surface changes of glaciers increased tremendously.Most detailed studies concentrated on the evolution of ice flow velocity during the surge (e.g., [15][16][17][18], with only few examples of time series of elevation changes during a surge (e.g., [19,20]).
In general, changes in glacier elevation, and thus volume, are interpreted as glacier mass changes, indicative of variations in the environmental conditions, and thus ultimately related to a change in climate (e.g., [21]).Surge-type glaciers, however, show distinct mass redistribution, independent of direct climate signals.In fact, surge mechanisms are attributed to changes either in the hydrological or thermal regime of the glacier [22].In both cases, the direct observation of the driving forces is almost impossible, because the mechanisms act predominantly at the base of the glaciers.In the case of thermally driven glacier surges, a switch from cold to temperate conditions in polythermal glaciers leads to a change in the basal force balance (e.g., [23,24]).Ice initially frozen to the bedrock is mobilised by the switch to temperate conditions and the subsequent existence of meltwater.For hydrologically-driven glacier surges, a sudden increase of englacial and basal meltwater meets an ineffective drainage basin, which results in high basal water pressures and consequently high sliding velocities [25,26].
The Pamir Mountain Range is known to host many surge-type glaciers [12].Based on optical space images, Kotlyakov et al. [5] identified 55 glaciers that surged between 1973 and 2006, and 215 glaciers that showed signs of recent instability.The best studied surging glacier in the region so far is the Medvezhiy Glacier, which was intensively studied regarding both velocity and elevation changes during a complete surge cycle from 1973 to 1989 [9] with possible hydrological forcing of the surges discussed in [27].In another study, surge-type glaciers were detected by comparing digital elevation models from 2000 and 2011 along the Pamir-Karakoram-Himalaya range [13].The identified glaciers feature large characteristic elevation changes with either thickening in the upper part and simultaneous thinning on the lower tongue in their quiescent phase, or vice-versa during the active phase.They found six actively surging glaciers in the surroundings of Fedchenko Glacier and at least seven quiescent ones, including Bivachny Glacier, whose most recent surge we will discuss in this paper.Bivachny Glacier is a tributary glacier of Fedchenko Glacier in the central Pamir.It is ~28 km long and covers an area of 96 km 2 .The glacier extends from 5100 m elevation in the south-eastern mountain range to 3400 m at its confluence with Fedchenko Glacier.The region is characterised by continental climate with cold winters and warm summers, where typical mean summer temperatures reach 2.4 • C and −17 • C during winter (based on historical observations at Gorbunov Station, 4200 m altitude, 17 km upstream of the Bivachny-Fedchenko confluence at the left bank of Fedchenko Glacier).The precipitation mainly occurs in winter and spring, while summers are very dry.The lower part of Bivachny Glacier is covered by supraglacial debris up to an elevation of about 3900 m.The late summer snow line is usually observed in an altitude range of 4600 to 4800 m, as observed in Landsat and on neighbouring Fedchenko in the field.
Bivachny Glacier is known to have surged thrice in the 20th century, in 1958,[1976][1977][1978] and "soon after" 1991 [5].The most evident sign of a surge at Bivachny is the displacement of the looped moraine of the MGU Glacier tributary (Figure 1; see Figure 2 for the most recent surge).The occurrence of a new surge in 2012-2015 was reported by Kotlyakov et al. [28] based mainly on the morphological interpretation of optical images taken aboard the International Space Station.
For this quantitative investigation of the recent surge, we combine high-resolution digital elevation models and multi-sensor feature tracking to describe and analyse the evolution of the active phase of the Bivachny Glacier surge between 2011 and 2015.In particular, the joint analysis of ice surface elevation and velocity changes for different periods during the surge gives valuable insights into surge dynamics and describes for the first time the surge-related ice volume transport.The looped moraines downstream of MGU Glacier are marked in white and black, respectively.

Materials and Methods
The study of the recent glacier surge at Bivachny Glacier is based on the interpretation of optical and synthetic aperture radar (SAR) remote sensing data, covering the entire surge phase from 2011 until 2015.Digital elevation models (DEMs) are derived from TanDEM-X data by interferometric SAR processing (InSAR) for the surge period.They are compared to the pre-surge shuttle radar  The looped moraines downstream of MGU Glacier are marked in white and black, respectively.

Materials and Methods
The study of the recent glacier surge at Bivachny Glacier is based on the interpretation of optical and synthetic aperture radar (SAR) remote sensing data, covering the entire surge phase from 2011 until 2015.Digital elevation models (DEMs) are derived from TanDEM-X data by interferometric SAR processing (InSAR) for the surge period.They are compared to the pre-surge shuttle radar

Materials and Methods
The study of the recent glacier surge at Bivachny Glacier is based on the interpretation of optical and synthetic aperture radar (SAR) remote sensing data, covering the entire surge phase from 2011 until 2015.Digital elevation models (DEMs) are derived from TanDEM-X data by interferometric SAR processing (InSAR) for the surge period.They are compared to the pre-surge shuttle radar topography mission (SRTM) global elevation model [29].In addition, surface displacements of the glacier are derived from optical data by correlation methods and transferred to surface velocities.

InSAR-Derived Digital Elevation Models
The DEM of the SRTM flown in February 2000 [29] is well established as a reference for glaciological studies.The mission operated two independent radar systems: a C-band SAR, which produced the contiguous DEM widely used, and an X-band SAR producing a DEM with reduced coverage due to a narrower swath width [30].
A number of studies compare newer DEMs from different sources to the state in 2000, in order to derive geodetic mass balance information [13,[31][32][33].There are some aspects, like DEM accuracy and radar signal penetration, that have to be evaluated carefully in order to obtain meaningful results.For Eurasia, the accuracy of the DEM is specified by a 90% absolute height error of 6.2 m [34], which is sufficient for mass balance studies on a decadal scale, but errors can be higher in the extremely rugged terrain in the Pamir Mountains.However, the SRTM DEM-as with any InSAR-derived elevation dataset-refers to the phase scattering centre of the mapped surface that does not necessarily coincide with the real glacier surface due to signal penetration [35].Penetration depth on glaciers depends amongst other factors on radar wavelength and surface properties, with higher penetration for longer wavelength and dry snow surfaces [36,37].
On debris cover, radar penetration is negligible due to strong surface scattering whereas on cold bare ice C-band penetration depth is typically in the order of 1 to 2 m [37].For snow-covered surfaces, penetration depends on the liquid water content.While C-band penetrates several meters into dry snow [37,38], wet snow absorbs the radar energy at the surface, rather resulting in a low backscatter than in a significant penetration depth.There are fewer studies about X-band penetration.Usually penetration of such short wavelengths is assumed negligible [36].
For comparison of elevation data sets from different sources, only the differential effect has to be taken into account.For the evaluation of the difference between C-and X-band penetration depths, the continuous C-band DEM can be compared within the 50-km wide stripes of the SRTM X-band DEM acquired simultaneously [30].Because the SRTM mission took place in winter in presumably dry conditions, penetration at each wavelength can be assumed maximal.For Northern Zulumart and neighbouring glaciers, a suitable location with a common coverage (about 50 km east of Fedchenko Glacier), the elevation difference is less than 2 m in the snow-covered ablation zone and increases to 4 m at 5300 m in the accumulation zone.For comparison, elevation differences on bare ground nearby oscillate around zero, not showing any elevation dependence.
The TanDEM-X bistatic data used in this study were acquired between 2011 and 2015 (Table 1).Data were processed into Raw DEM products by the Integrated TanDEM-X Processor (ITP) [39].This procedure includes all steps necessary to interferometrically process the data into geocoded height maps including co-registration, phase unwrapping and geocoding.Phase unwrapping is especially challenging in the study area, due to the steep terrain which results in a large fraction of shadowing or layover.Additionally, depending on surface conditions, the backscatter of the glacier can be very low in summer months, resulting in a low signal-to-noise ratio.In the operational processing of the global TanDEM-X DEM, the dual baseline phase unwrapping algorithm combines two interferograms with different baselines typically acquired in different years according to the acquisition plan of the global DEM [40].This approach cannot be adopted here because each acquisition represents an independent snap shot in time of the changing glacier topography.
In a standard ITP run for a single interferogram, the stereo-radargrammetric measure based on the geometrical parallaxes from the co-registration step is used to estimate the absolute phase offset of the largest continuously unwrapped image patch.While this approach works well in even terrain, for the study area it can result in the determination of the correct phase offset for a portion of the image not necessarily including the targeted glacier.Therefore, based on the comparison of the automatically unwrapped DEM with the SRTM DEM and supposing a mean height difference of less than height of ambiguity (H a ), the appropriate absolute phase offset was determined for each continuously unwrapped image patch inside the glacier.Using this newly-derived absolute phase offset, the steps of phase-to-height conversion and geocoding were repeated within ITP.Next, the thus unwrapped DEM patches were mosaicked to derive the DEM.Steep slopes causing layover had to be masked out based on the flag map provided by the processor that was edited semi-automatically to exclude erratic height values.To assess the accuracy of the derived DEMs, pairwise elevation differences on bedrock adjacent to the glacier were calculated and vertically adjusted to the same level.The maximum correction with respect to the mean level of all DEMs was detected for the DEM of 5 January 2014 with a vertical shift of 2.85 m.The root mean square (RMS) of consecutive differences on bedrock after the adjustment lies between 2.7 m and 3.9 m, but increases to 8.1 m for pairs involving the DEM of 15 October 2015 due to the larger baseline.
Although the same instrument acquired all DEMs during the surge, radar penetration has to be considered because of changing snow properties.While the presence of liquid water in wet summer snow will inhibit penetration, the scattering centre can be located several meters below the surface of dry winter snow, resulting in the mapping of the previous summer horizon in the ablation area.In the accumulation area of the neighbouring Fedchenko Glacier, melting advances upglacier during spring and summer, eventually reaching even the upper accumulation basin at 5200 m, documented by low backscatter values during summer.In some places, elevation coincides in the spring/summer acquisitions (August 2011, May 2013 and May 2014) within 0.1 m, but is about 6 m lower in winter acquisitions (January 2012, February 2013 and January 2014).These differences are unlikely to be caused by a real surface elevation change but by deeper radar penetration in winter.Maximum elevation differences occur when a winter DEM with maximum penetration is combined with a spring or summer DEM with melting and thus low backscattering at the ice surface.
Accordingly, winter acquisitions can have a negative bias in the accumulation area in the order of up to 6 m.As most of the surge under consideration here takes place in the ablation area on only in winter snow-covered bare ice or debris and considering the magnitude of the elevation change to be observed, radar wave penetration does not play a significant role in this analysis.

Feature Tracking of Optical and SAR Data
Velocity information is derived from feature tracking, a pairwise cross-correlation that can be applied to optical as well as SAR data (e.g., [41][42][43]).Optical data were added to the study in order to derive velocity information with a higher temporal resolution during the surge and to extend the series backwards in time (Table 2).The main prerequisite for feature tracking is the stability of surface features.For optical data, cloud-free scenes acquired in summer/autumn with minimal snow cover extent were selected from the U.S. Geological Survey Landsat archive (https://earthexplorer.usgs.gov/).These images can be downloaded in an orthorectified and terrain-corrected format based on the SRTM DEM [44].To minimise geometrical distortions, only scenes from the same path/row were used.For SAR data the restriction to summer and autumn does not strictly apply because the dry winter snow is penetrated and surface features below the winter snow layer are discernible in the amplitude images.Correlating SAR scenes from different seasons can nevertheless be hindered by strong backscatter differences due to changes in absorption caused by the presence of liquid water in summer.
The displacement fields were derived using the image cross-correlation algorithm (e.g., [42]) implemented in "System for Automated Geoscientific Analyses" SAGA GIS [45].The algorithm performs a normalised cross-correlation in the space domain.The reference window size and the search window size, within which the reference window is incrementally shifted to determine maximum correlation, can be defined independently.Thus, small patches can be located in a wider window, allowing to better resolve the shear margins between stable rock and fast-moving ice.For Landsat 5, near infrared band 4 was selected while for Landsat 7, Landsat 8 and EO-1 panchromatic data were used to make use of the higher spatial resolution.A reference window of 32 × 32 pixels and a search window of 64 × 64 or 128 × 128 pixels were used, depending on the displacement magnitude within the time span.In the upper part of Bivachny Glacier, feature tracking fails due to the lack of distinct features and the temporary presence of shadows from the nearby mountain.This is especially true for the medium-resolution optical data.In this region features where tracked manually to provide additional sampling points and to check the results.
Due to the higher resolution of the TerraSAR-X data, smaller displacements and larger gradients can be mapped.The applied reference window size of 64 × 64 pixels corresponds to a patch of 80 × 80 m for the TerraSAR-X data.Geocoded SAR amplitude images of the active satellite of the TanDEM-X acquisitions were also used for feature tracking with a reference window size of 64 × 64 pixels.
The accuracy of the resultant velocity fields depends on the correlation quality of the respective patch, but also on pixel resolution and time span between the acquisitions.To assess the achieved accuracies, correlation results on adjacent stable bedrock were evaluated.Table 2 lists the 90th percentile of the derived bedrock displacements converted into velocities in m per day as a robust measure of uncertainties in velocity determination.As expected, errors are larger for shorter time spans and larger pixel size.This error already includes possible coregistration residuals of the standard orthorectification.Another error source that is not included in the error estimation on bedrock is inadequate height information used for georeferencing.Inadequate here means referring to another date.The large temporal elevation change caused by the surge results in a georeferencing error on the glacier when the SRTM DEM of 2000 is used.While an absolute height error simply leads to a spatial shift of the derived velocity vector in the order of few meters, an elevation change between the two acquisitions used for a velocity field produces a systematic velocity error.For the geometric position of Bivachny Glacier within the Landsat scenes, an elevation difference of 100 m would result in a dislocation of 11 m perpendicular to the centre line of the scene which is then projected onto the flow direction.With a maximum elevation change rate of 0.5 m day −1 observed during the surge, the velocity error introduced is less than 22 m year −1 and thus within the error budget of the short-term Landsat pairs used during the surge.To test the influence of elevation changes on the velocity estimation in the case of SAR data, the amplitude images of the consecutive TanDEM-X acquisitions containing the largest elevation change (13 May 2014 and 11 September 2014) were used for the velocity estimation once geocoded with their respective DEMs and a second time with SRTM.The difference in velocity is roughly proportional to the elevation change between the two acquisitions.For a 60-m increase in elevation, velocities were underestimated by 10 per cent for this image pair.This particular case can be regarded as a worst case as elevation changes are smaller in most observed cases and the small angle between SAR look direction and flow direction maximises the effect of the wrong elevation used for geocoding on the velocity.Decreasing elevation was observed less frequently during this time span and the correlation with velocity differences is less clear.Larger elevation change rates presumably lead to a noncoherent displacement, preventing successful velocity determination.

Previous Surge
From comparison of the sparsely available satellite imagery from 1991 and 2003, Kotlyakov et al. [5] suspected that a surge happened in the early 1990s.Our analysis of Landsat images from this period that have become available since then reveals that this surge took place in 1996/1997.It was of minor magnitude with an activated zone of about 7.5 km in length and annual velocities reaching 1000 m year −1 .As the most prominent feature of the surge, the MGU lobe advanced by 1700 m.In general, the surge was confined to the upper part of the ablation area of Bivachny Glacier down to the confluence with Oshanin Glacier.After the end of the surge, ice flow velocities dropped by an order of magnitude.

Pre-Surge Conditions
By the year 2000-when the SRTM DEM was acquired-Bivachny Glacier turned almost inactive.Maximum annual velocities derived from Landsat 7 scenes (August 2000 and September 2001) amounted to about 95 m year −1 in the first two kilometres of the central flowline starting at the confluence of the four accumulation basins (see red profile in Figure 1).At 15 km along the flowline (about eight kilometres upstream of the confluence with Fedchenko Glacier) there was a second velocity maximum of about 88 m year −1 .Along the rest of the profile, velocities usually reached less than 40 m year −1 .The SRTM DEM can thus be regarded as reference topography during the quiescent phase before the onset of a new surge.The longitudinal profile of the glacier during this phase shows an even surface, continuously decreasing in elevation, without pronounced elevation steps (black line in Figure 3).
Based on this stable geometry and the corresponding velocities we estimated the ice thickness along the central flowline, using the same approach already applied for the analysis of elevation and velocity changes at the neighbouring Fedchenko Glacier [46].The bedrock elevation, derived from the ice thicknesses is later used for the analysis of ice flux changes during the surge.The model assumes that the driving stress is only balanced by the basal shear stress during the quiescent phase.Then, ice thickness can be derived from the surface slope α and the deformation velocity u def (the difference between surface and basal velocity) with the parameters of Glenn's flow law for temperate ice A and n, density of ice ρ and the gravitational acceleration g: h = (n+1) u de f ( The implicit assumptions of constant glacier width and negligible inflow from the tributaries are approximately met.Neglecting basal motion, i.e., subglacial sediment deformation and basal sliding, during the quiescent phase, the deformation velocity u def can be set to the surface velocity.Surface velocities and surface slopes were smoothed over a distance of 1 km with a running mean filter, before calculating the ice thicknesses.The resulting ice thicknesses vary between 190 m and 370 m along the flow line, with a mean value of 280 m.These values can be regarded as maximum estimations, because potential basal motion will lower the results.For example, if half of the measured surface velocity is due to basal processes, ice thickness estimates will be 15 percent lower.A peculiar feature in the resulting bedrock profile is the noticeable bump around 6 km (see brown line in Figure 3).A sensitivity analysis showed that this rise is robust against velocity uncertainties and the choice of spatial filtering.
Remote Sens. 2017, 9, 388 8 of 17 (the difference between surface and basal velocity) with the parameters of Glenn's flow law for temperate ice A and n, density of ice ρ and the gravitational acceleration g: 1 . ( The implicit assumptions of constant glacier width and negligible inflow from the tributaries are approximately met.Neglecting basal motion, i.e., subglacial sediment deformation and basal sliding, during the quiescent phase, the deformation velocity udef can be set to the surface velocity.Surface velocities and surface slopes were smoothed over a distance of 1 km with a running mean filter, before calculating the ice thicknesses.The resulting ice thicknesses vary between 190 m and 370 m along the flow line, with a mean value of 280 m.These values can be regarded as maximum estimations, because potential basal motion will lower the results.For example, if half of the measured surface velocity is due to basal processes, ice thickness estimates will be 15 percent lower.A peculiar feature in the resulting bedrock profile is the noticeable bump around 6 km (see brown line in Figure 3).A sensitivity analysis showed that this rise is robust against velocity uncertainties and the choice of spatial filtering.

Surge Build-Up
Landsat 5 and 7 images acquired until 2010 show little ice flow activity of Bivachny Glacier, indicating that the glacier was in a "normal mode" during this period (Figure 4b).Downstream of the confluence of the MGU Glacier (7 km on the central flowline profile), velocities stayed below 100 m year −1 .During this quiescent phase, the MGU tributary advanced into the depleted, stagnant Bivachny Glacier until occupying half of the glacier width on the surface in 2011.On the lowest part of the glacier, velocities were even reduced in the pre-surge period 2008-2010 compared to 2000 (20 m year −1 mean value downstream of 12 km on the profile).However, the region upstream of the MGU

Surge Build-Up
Landsat 5 and 7 images acquired until 2010 show little ice flow activity of Bivachny Glacier, indicating that the glacier was in a "normal mode" during this period (Figure 4a).Downstream of the confluence of the MGU Glacier (7 km on the central flowline profile), velocities stayed below 100 m year −1 .During this quiescent phase, the MGU tributary advanced into the depleted, stagnant Bivachny Glacier until occupying half of the glacier width on the surface in 2011.On the lowest part of the glacier, velocities were even reduced in the pre-surge period 2008-2010 compared to 2000 (20 m year −1 mean value downstream of 12 km on the profile).However, the region upstream of the MGU Glacier confluence (0-7 km on the profile) already showed a higher activity with velocities of 150-180 m year −1 .The MGU confluence marked the transition between elevation increase upstream (mean elevation increase 18 m, 8% of the estimated ice thickness) and elevation decrease downstream from February 2000 to August 2011, when the first TanDEM-X DEM was acquired (Figure 4b).Although the exact moment of the initiation of the surge cannot be precisely inferred from the available data, it can be assumed that the DEM from August 2011 depicts a very early stage of the surge.
Remote Sens. 2017, 9, 388 9 of 17 Glacier confluence (0-7 km on the profile) already showed a higher activity with velocities of 150-180 m year −1 .The MGU confluence marked the transition between elevation increase upstream (mean elevation increase 18 m, 8% of the estimated ice thickness) and elevation decrease downstream from February 2000 to August 2011, when the first TanDEM-X DEM was acquired (Figure 4a).Although the exact moment of the initiation of the surge cannot be precisely inferred from the available data, it can be assumed that the DEM from August 2011 depicts a very early stage of the surge.In the region between the MGU and the Oshanin Glacier confluence (between 9 km and 13 km on the profile), a strong elevation reduction was observed, with a mean of 43 m across the glacier and up to 90 m in the debris-free parts of this region in the period 2000-2011.This was also the region of the lowest surface velocities of about 10 m year −1 where the glacier seemed to be almost cut off from upstream ice supply.In contrast, the surface elevation decreased by only 15 m in the lower, completely debris-covered part downstream of the Oshanin confluence.This is three times less than 3 km further upstream.Considering the rather low ice velocities in this region, it can be speculated that a major part of the strong surface lowering results from ice loss by the surface mass balance.The mean surface ablation at the elevation of 4000 m (10 km on the profile) would then be at least −3.7 m year −1 for a surface with 81% supraglacial debris cover and up to 8 m year −1 for the clean ice areas.

The Surge
The surge started at some time between the summers 2010 and 2011 when velocities on the upper reaches of the central flowline increased to 350 m year −1 and the active zone of clearly increased velocities expanded past the MGU confluence.Additional manual tracking of single features distinguishable in all available Landsat images indicated a transition from a gradual velocity increase above the MGU confluence in the period 2006-2010 to a strongly intensified acceleration.
The surge then became visible in optical satellite images as the MGU lobe started being moved noticeably down-glacier.For a quantitative analysis, the evolution of surface topography and ice flow velocities throughout the surge can be monitored using the TanDEM-X DEMs available at irregular intervals and velocity fields that correspond approximately to the same dates (Figure 4).Detailed maps of elevation change and velocity fields can be found in the Supplementary Materials Figure S1. Figure 5 shows the total elevation change comprising the whole surge.
From autumn 2011 until January 2012, the surface elevation decreased along the first 5 km of the profile, corresponding to the upper part of the reservoir area, with a mean elevation decrease of 8 m during those five months.The entire flowline upstream of 12.2 km was already activated in the period autumn 2010 to autumn 2011, and the velocities in late summer 2011 were considerably increased compared to this annual mean (700 m year −1 versus 440 m year −1 ).There was a strong surface velocity gradient between 8.8 and 12.2 km along the profile, where the velocity decreased from 520 m year −1 to practically stagnant.Simultaneously, a strong elevation increase occurred on the same segment with a mean elevation gain of 22 m between August 2011 and January 2012.The maximum elevation change occurred where the strain rate was highest.Until January 2012, the remaining part of the glacier downstream of 12.2 km was not affected at all by the activation of the upper ablation zone.Except within the advancing compression front, velocities in autumn 2011 and the annual velocities 2011-2012 were comparable, indicating quite constant velocities throughout the whole year without pronounced seasonality.Maximum annual mean velocities reached up to 750 m year −1 .No elevation model is available for summer or autumn 2012, but in February 2013, the large negative ice thickness anomaly of the quiescent phase between 9 and 13 km was already compensated.Ice thicknesses were then about 20 m larger than in 2000; compared to a 40-m lower surface in August 2011.The maximum elevation increase amounted to 75 m at 13 km, with a minor bulge of a 15-m increase at 16 km.In the following three months (February 2013 until May 2013), the surface elevation increase maintained the same moderate rate of 0.2 m day −1 and spread downstream until 19 km.While there was no significant change in surface elevation observed on the Bivachny profile upstream 13 km, the elevation on the lowermost 2 km of Oshanin Glacier decreased by 10 m.Ice flow velocities in spring 2013 stayed at the same level of about 600 m year −1 above the MGU confluence and reached about 800 m year −1 at 15 km before dropping to zero at 19 km.On Oshanin Glacier itself, velocities rose to about 440 m year −1 .From January to May 2014, maximum elevation change rates decreased to 0.37 m day −1 .The region of elevation increase stretched to a length of 9.4 km, while the elevation decrease was still restricted to the upper 14 km of the profile.The highest velocities were again observed in late spring/early summer from May to June 2014 with a maximum of 3400 m year −1 around 16 km.Velocities at least doubled in comparison to the autumn 2013 velocities downstream of 10 km.Compressional strain rate peaked at 22.6 km with a maximum −0.0037 day −1 .
In September 2014, the surface bulge finally reached the confluence with Fedchenko Glacier.Elevation change rates were at the highest level showing a maximum elevation increase of 69 m at the same position (0.50 m day −1 ).The depletion area then extended downstream to the region below Oshanin Glacier (surface lowering of up to 40 m from May 2014 until August 2015).Between September 2014 and October 2015 there was still some activity, with an elevation increase of up to 40 m on the lowest 3.5 km and even some ice flow into Fedchenko Glacier.At the southern side of the confluence, the surge bulge clearly overran the flat lateral moraine of Fedchenko Glacier by about 500 m.A photo taken in August 2015 shows Bivachny Glacier pushing into Fedchenko Glacier (Figure 6).From January to May 2014, maximum elevation change rates decreased to 0.37 m day −1 .The region of elevation increase stretched to a length of 9.4 km, while the elevation decrease was still restricted to the upper 14 km of the profile.The highest velocities were again observed in late spring/early summer from May to June 2014 with a maximum of 3400 m year −1 around 16 km.Velocities at least doubled in comparison to the autumn 2013 velocities downstream of 10 km.Compressional strain rate peaked at 22.6 km with a maximum −0.0037 day −1 .
In September 2014, the surface bulge finally reached the confluence with Fedchenko Glacier.Elevation change rates were at the highest level showing a maximum elevation increase of 69 m at the same position (0.50 m day −1 ).The depletion area then extended downstream to the region below Oshanin Glacier (surface lowering of up to 40 m from May 2014 until August 2015).Between September 2014 and October 2015 there was still some activity, with an elevation increase of up to 40 m on the lowest 3.5 km and even some ice flow into Fedchenko Glacier.At the southern side of the confluence, the surge bulge clearly overran the flat lateral moraine of Fedchenko Glacier by about 500 m.A photo taken in August 2015 shows Bivachny Glacier pushing into Fedchenko Glacier (Figure 6).

Surge Termination
After the observed ice speed maximum in May/June 2014 (3400 m year −1 at 16.1 km), velocities dropped below the level of the previous year by August 2014 (1100 m year −1 versus 1400 m year −1 in 2013 at the same position), but still the entire flow line was active and showed enhanced ice transport.Velocities throughout 2015 showed a continued slowdown without the pronounced spring acceleration observed in the three previous years.Finally, velocities did not exceed 180 m year −1 along the entire length of the profile in September 2015.DEMs of October and December 2015 confirmed the end of the surge without any signs of elevation changes above the error level of the DEMs.There is no abrupt stop of the surge, but rather a tapering off during the summer 2015.

Discussion
The most recent surge of Bivachny affected the whole ablation zone down to the confluence with Fedchenko Glacier (Figure 5).It started in 2011 above the MGU confluence and reached the confluence with Fedchenko Glacier in 2014, displacing its western margin by about 500 m.Velocities increased from below 90 m year −1 during quiescence in 2000 up to 3400 m year −1 in spring 2014.A surface bulge up to 80 m high travelled at mean speed of 12 m day −1 (4400 m year −1 ) downstream for a distance of 13 km.
To separate the elevation change effect of the surge from the mass balance signal, Voroshilov Glacier (see Figure 1), a tributary not directly affected by the surge, was used as an analogue.For this small glacier the annual elevation change rate of −0.84 m year −1 was rather constant during the surge with positive deviations in summer and negative ones in winter due to higher penetration in cold winter snow and ice.At one of the accumulation basins, the mean elevation throughout the whole surge period decreased linearly from −19 m year −1 at the beginning of the Bivachny flow line to −4.4 m year −1 at an elevation of 4900 m.Thus, even in the accumulation area of Bivachny, the rate of elevation decrease was 5 times that of an unaffected glacier.
Kotlyakov et al. [5] postulated that the tributary MGU glacier dammed the upper part of Bivachny and a surge was hence triggered whenever the driving stress exceeded the resistance of the

Surge Termination
After the observed ice speed maximum in May/June 2014 (3400 m year −1 at 16.1 km), velocities dropped below the level of the previous year by August 2014 (1100 m year −1 versus 1400 m year −1 in 2013 at the same position), but still the entire flow line was active and showed enhanced ice transport.Velocities throughout 2015 showed a continued slowdown without the pronounced spring acceleration observed in the three previous years.Finally, velocities did not exceed 180 m year −1 along the entire length of the profile in September 2015.DEMs of October and December 2015 confirmed the end of the surge without any signs of elevation changes above the error level of the DEMs.There is no abrupt stop of the surge, but rather a tapering off during the summer 2015.

Discussion
The most recent surge of Bivachny affected the whole ablation zone down to the confluence with Fedchenko Glacier (Figure 5).It started in 2011 above the MGU confluence and reached the confluence with Fedchenko Glacier in 2014, displacing its western margin by about 500 m.Velocities increased from below 90 m year −1 during quiescence in 2000 up to 3400 m year −1 in spring 2014.A surface bulge up to 80 m high travelled at mean speed of 12 m day −1 (4400 m year −1 ) downstream for a distance of 13 km.
To separate the elevation change effect of the surge from the mass balance signal, Voroshilov Glacier (see Figure 1), a tributary not directly affected by the surge, was used as an analogue.For this small glacier the annual elevation change rate of −0.84 m year −1 was rather constant during the surge with positive deviations in summer and negative ones in winter due to higher penetration in cold winter snow and ice.At one of the accumulation basins, the mean elevation throughout the whole surge period decreased linearly from −19 m year −1 at the beginning of the Bivachny flow line to −4.4 m year −1 at an elevation of 4900 m.Thus, even in the accumulation area of Bivachny, the rate of elevation decrease was 5 times that of an unaffected glacier.
Kotlyakov et al. [5] postulated that the tributary MGU glacier dammed the upper part of Bivachny and a surge was hence triggered whenever the driving stress exceeded the resistance of the dam.We speculate that there were low ice flow velocities in the upper reaches of Bivachny during the quiescent phase, and thus the limited volume flux could not effectively counterbalance the mass gain of the accumulation area above, causing the ice to thicken and enabling the recovery of the upper ablation zone from the strong mass losses of the previous surge.The blocking of the MGU glacier and the modelled bedrock bump at 6 km (Figure 3) could have played an additional role in slowing down the ice flux.The importance of a similar bedrock ridge for the surge build-up of a surge-type glacier in Yukon, Canada has been investigated by authors of [47].
The increased driving stress of this mean thickening of 8% above the MGU confluence was not sufficient to cause the observed increase in flow velocity from about 9 m year −1 to 640 m year −1 between 2000 and 2011.The corresponding increase in deformation velocity (assuming temperate conditions) would only be from 4 m year −1 to 9 m year −1 .However, the driving stress showed a moderate increase from 120 kPa to about 150 kPa between 11 km and 13.5 km on the profile just upstream of the Oshanin Glacier confluence in 2000.During the early surge state, it decreased sharply from about 160 kPa to 40 kPa over the same part of the profile in summer 2011, with a rise to 185 kPa over a distance of just 1.5 km further downstream.This expressed local minimum in driving stress indicates that there existed a (temporal) barrier for the ice flow at the confluence with Oshanin Glacier in summer 2011.The seasonal varying ice flow velocities with the observed spring acceleration for all years of the surge, suggests a hydrological trigger.The region where the increase in ice flow started is in the altitude range of 4500 to 4200 m, several hundred meters below the equilibrium line.Therefore, an abundance of meltwater can be expected in the ablation season, while considerable surface water run-off is not detectable on satellite imagery and crevasses are numerous.This indicates that meltwater usually is drained by a developed subglacial channel system.The increased overburden pressure in the upper reaches due to the ice build-up could cause the collapse of the channelized subglacial water transport.A switch from effective meltwater drainage within well-developed channels to a distributed cavity system could explain an increased sliding velocity.Another explanation for surge initiation could be that the higher surface slope from the region of increasing ice thickness to the area of decreasing ice thickness, potentially caused by MGU blocking the Bivachny Glacier, led to higher ice flow velocities.If the build-up is sustainable, a switch in the drainage system could be reached and thus a further increase in ice flow ( [15,22] p. 328).This evolution propagates then along the glacier, switching the basal regime to more ineffective drainage conditions.
The surge can be divided into two phases.In the first phase, ice velocities in the upper reaches gradually increased and the activated region with higher ice flow velocities slowly progressed down glacier.The typical surge bulge developed where the compression between the enhanced flow and the stagnant ice downstream was highest.Elevation increase reached about 0.2 m day −1 during this first stage.By February 2013, the surface depression upstream of the Oshanin confluence was filled and the pre-surge elevation decrease was levelled out.At this point, the previous surge in 1997 had stopped.Apparently, the recent surge was more powerful than its predecessor: in spring 2013, when ice flow accelerated again the surge bulge passed the confluence with the Oshanin Glacier.We do not have any volume flux information about the previous surge, but it seems justified to assume that the recent surge mobilised more glacier volume, so that the ice volume deficit just upstream of the Oshanin confluence could be more than compensated and subsequently the surge conditions could evolve downstream.In summer 2013, the minor bulge that developed 3 km downstream of the Oshanin confluence intensified due to higher velocities and elevation change rates above 0.4 m day −1 .There is a clear seasonal cycle in the velocities with maximum values in early summer before the onset of melting.Consequently, maximum elevation change rates in late winter/ early spring 2014 were slightly lower.In early summer 2014, the highest ice flow velocities were observed together with the highest elevation change rates, but on a shorter section.Upstream of the only 4-km-long zone of elevation increase, the surface lowered at the highest rates observed.By September 2014, the surge front reached the confluence with Fedchenko Glacier.The supraglacial lake at the northern side of the confluence with Fedchenko Glacier drained between 1 and 17 October, as witnessed on cloud-free Landsat 8 images.Velocities further decreased and the surface bulge was compressed to a length of only 3.5 km, while strong thinning occurred immediately behind the surge bulge.Over the subsequent months, the surge slowly faded off with monotonically decreasing velocities without spring acceleration.The DEM difference from October to December 2015 showed that the active phase completely ended.
To assess the magnitude of the surge and the mass movement involved, ice volume fluxes were estimated (Figure 7).We used a simplistic approach that assumes a constant glacier width with an elliptical cross section and neglects inflows from the tributaries.Ice flow velocities were considered uniform across the cross section due to block sliding.Ice thicknesses were calculated based on the DEM closest in time and the estimated bedrock geometry along the central flowline (Equation (1) and Figure 3).The time difference between velocity and elevation determinations is largest during the quiescent phase when change rates are low.During the surge, the elevation changes within the maximum time difference of two months amount to less than 30 m, i.e., 10 per cent of the ice thickness.Consequently, this time difference is only a minor source of error in this estimate.As expected the maximum ice volume flux is observed during maximum ice flow velocities in May/June 2014 and amounted to 2.9 × 10 6 m 3 day −1 .This is forty times as much as the ice flux during quiescence in 2000 calculated with the same assumptions.Curves are similar to the ones of ice flow velocities, but influenced by the ice thickness distribution.The ice volume in the lower 13 km of the glacier increases by 0.72 km 3 in total and gives a measure for the amount of ice relocated during the surge.This is three times the maximum annual ice volume flux of the considerably larger Fedchenko Glacier [46].
Remote Sens. 2017, 9, 388 14 of 17 of only 3.5 km, while strong thinning occurred immediately behind the surge bulge.Over the subsequent months, the surge slowly faded off with monotonically decreasing velocities without spring acceleration.The DEM difference from October to December 2015 showed that the active phase completely ended.
To assess the magnitude of the surge and the mass movement involved, ice volume fluxes were estimated (Figure 7).We used a simplistic approach that assumes a constant glacier width with an elliptical cross section and neglects inflows from the tributaries.Ice flow velocities were considered uniform across the cross section due to block sliding.Ice thicknesses were calculated based on the DEM closest in time and the estimated bedrock geometry along the central flowline (Equation (1) and Figure 3).The time difference between velocity and elevation determinations is largest during the quiescent phase when change rates are low.During the surge, the elevation changes within the maximum time difference of two months amount to less than 30 m, i.e., 10 per cent of the ice thickness.Consequently, this time difference is only a minor source of error in this estimate.As expected the maximum ice volume flux is observed during maximum ice flow velocities in May/June 2014 and amounted to 2.9 × 10 6 m 3 day −1 .This is forty times as much as the ice flux during quiescence in 2000 calculated with the same assumptions.Curves are similar to the ones of ice flow velocities, but influenced by the ice thickness distribution.The ice volume in the lower 13 km of the glacier increases by 0.72 km 3 in total and gives a measure for the amount of ice relocated during the surge.This is three times the maximum annual ice volume flux of the considerably larger Fedchenko Glacier [46].

Conclusions
Our analysis confirms that the combination of surface velocity and elevation observations from multi-sensor satellite sources allows a detailed investigation of changes in even short-term glacier behaviour.In particular, the availability of multiple DEMs per year enables the monitoring of the progression of the surge and the detailed calculation of mass redistribution for Bivachny Glacier.Surface velocities increased from near-zero during the quiescence phase up to 3400 m year −1 in spring 2014 with superimposed seasonal variations.In contrast to what is observed at other surges (e.g., [18]), a sharp surge front developed during the surge of Bivachny Glacier, which is traceable in elevation as well as strain rate.It travelled with a mean velocity of 4400 m year −1 for a distance of 13 km down to the confluence with Fedchenko Glacier.This indicates that the basal conditions did not change simultaneously along the entire glacier tongue, but the activation of the fast ice flow occurred gradually.The seasonality and the evolution of the surge over several years both suggest that there exists a control on glacier flow that is not easily overcome.Based on our observations, only the upper 12 km of the profile were affected by higher velocities before summer 2012.At the lower end of this

Conclusions
Our analysis confirms that the combination of surface velocity and elevation observations from multi-sensor satellite sources allows a detailed investigation of changes in even short-term glacier behaviour.In particular, the availability of multiple DEMs per year enables the monitoring of the progression of the surge and the detailed calculation of mass redistribution for Bivachny Glacier.Surface velocities increased from near-zero during the quiescence phase up to 3400 m year −1 in spring 2014 with superimposed seasonal variations.In contrast to what is observed at other surges (e.g., [18]), a sharp surge front developed during the surge of Bivachny Glacier, which is traceable in elevation as well as strain rate.It travelled with a mean velocity of 4400 m year −1 for a distance of 13 km down to the confluence with Fedchenko Glacier.This indicates that the basal conditions did not change simultaneously along the entire glacier tongue, but the activation of the fast ice flow occurred gradually.The seasonality and the evolution of the surge over several years both suggest that there exists a control on glacier flow that is not easily overcome.Based on our observations, only the upper 12 km of the profile were affected by higher velocities before summer 2012.At the lower end of this region, the main mass loss occurred during the previous quiescent phase.Only after the according volume deficit was compensated could the surge expand to the region below, and the activated zone quickly evolved and covered another 10 km within just one year.When the entire glacier tongue was in surge motion in the beginning of the 2014 ablation season, the observed mass flux was highest but not sustainable.In the summer 2014 the mass flux had already reduced considerably and the surge terminated during the following winter.
Comparison with observations from other surging glaciers shows that this surge is similar to glacier surges in the Karakoram (e.g., [15,17]) with respect to development and duration.Initiation as well as termination of the surge happened gradually, over a period of several months.The active phase itself lasted about four years, with a gradual build-up of the ice velocities from year to year and an accompanied downward motion of the surge front.This is rather similar to the surge of the North Gasherbrum Glacier in the Karakoram from 2003 to 2006 [15].The seasonal slow-down during late summer and the regular speed-up during the onset of the melt period is a strong indication for a hydrological control of the surge, where meltwater availability influences the ice flow.In contrast, the Variegated Glacier in Alaska for example showed a much more vigorous surge evolution with a decidedly shorter activity period [7].Because the Bivachny Glacier is a tributary glacier of the much larger Fedchenko Glacier, no obvious advance of the tongue could be observed.The surging glacier did protrude into the Fedchenko main valley, overriding the western margin of Fedchenko Glacier before in slowed down, but the conditions were not set to initiate a speed-up of Fedchenko Glacier.Our investigation documents another glacier surge in High Asia in detail and indicates that surges on flat valley glaciers in dry High Asia with warm summer temperatures typically build up over a period of several years with a clear seasonality.The mass redistribution is mainly confined within the initial glacier boundaries with some terminus advances where the conditions allow (e.g., [17]).

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/9/4/388/s1, Figure S1: Maps of elevation change rate (or absolute elevation change for May 2014-October 2015) and velocities of Bivachny Glacier at discrete time steps.

Figure 1 .
Figure 1.Bivachny Glacier with its tributaries, Pamir Mountains.Background image: Landsat 7 (24 August 2000), glacier outline in black, longitudinal profile in red, contour lines from the shuttle radar topography mission (SRTM) 2000 in grey.Inset shows location map.

Figure 2 .
Figure 2. The displacement of surface features during the recent surge of Bivachny Glacier as seen in optical satellite images: (a) EO-1 image of the year 2011; (b) Landsat 8 image of the year 2015.The looped moraines downstream of MGU Glacier are marked in white and black, respectively.

Figure 3 .
Figure 3. Longitudinal elevation profile of the Bivachny Glacier central flowline in 2000 (SRTM in black), 2011 (TanDEM-X in blue) and 2015 (TanDEM-X in red) and the reconstructed bedrock topography based on equation 1 in brown.The brown arrow marks the bedrock bump mentioned in the text, black arrows mark the position of the MGU Glacier and Oshanin Glacier confluences, respectively.

Figure 3 .
Figure 3. Longitudinal elevation profile of the Bivachny Glacier central flowline in 2000 (SRTM in black), 2011 (TanDEM-X in blue) and 2015 (TanDEM-X in red) and the reconstructed bedrock topography based on equation 1 in brown.The brown arrow marks the bedrock bump mentioned in the text, black arrows mark the position of the MGU Glacier and Oshanin Glacier confluences, respectively.

Figure 4 .
Figure 4. (a) Surface velocity and (b) elevation change along the central flowline (see Figure 1) for selected periods.

Figure 4 .
Figure 4. (a) Surface velocity and (b) elevation change along the central flowline (see Figure 1) for selected periods.

Figure 5 .
Figure 5.Total elevation difference from August 2011 to October 2015.Colour scale is from blue for elevation loss to red for elevation gain.Background image: Landsat 8.The period from May until October 2013 showed again a strong redistribution of mass, this time affecting the entire flow line from the uppermost part to about 20.1 km.Elevation change rates peaked at 19 km (dh/dtmax = 0.43 m day −1 ), more than double in comparison to the three previous months.In early summer 2013 (May to July), the ice velocities reached maximum values of more than 2000 m year −1 between 10 and 16.5 km along the profile, before dropping below 1500 m year −1 in autumn.These seasonal variations are in accordance with velocity observations for Fedchenko as well as Medvezhiy Glacier [9,46].Even the upper part of the flow line experienced a strong increase in surface velocities, leading to dynamic thinning with a surface elevation drop of 15-20 m along the upper 10 km.Again, the lowermost part of Oshanin Glacier experienced an elevation decrease by 20 to 25 m.During winter 2013/14 the surge passed the Kalinin confluence at 19 km.Behind the advancing surge front, elevation change rates were similar to the ones in the months before, but with a higher spatial fluctuation.From January to May 2014, maximum elevation change rates decreased to 0.37 m day −1 .The region of elevation increase stretched to a length of 9.4 km, while the elevation decrease was still restricted to the upper 14 km of the profile.The highest velocities were again observed in late spring/early summer from May to June 2014 with a maximum of 3400 m year −1 around 16 km.Velocities at least doubled in comparison to the autumn 2013 velocities downstream of 10 km.Compressional strain rate peaked at 22.6 km with a maximum −0.0037 day −1 .In September 2014, the surface bulge finally reached the confluence with Fedchenko Glacier.Elevation change rates were at the highest level showing a maximum elevation increase of 69 m at the same position (0.50 m day −1 ).The depletion area then extended downstream to the region below Oshanin Glacier (surface lowering of up to 40 m from May 2014 until August 2015).Between September 2014 and October 2015 there was still some activity, with an elevation increase of up to 40 m on the lowest 3.5 km and even some ice flow into Fedchenko Glacier.At the southern side of the confluence, the surge bulge clearly overran the flat lateral moraine of Fedchenko Glacier by about 500 m.A photo taken in August 2015 shows Bivachny Glacier pushing into Fedchenko Glacier (Figure6).

Figure 5 .
Figure 5.Total elevation difference from August 2011 to October 2015.Colour scale is from blue for elevation loss to red for elevation gain.Background image: Landsat 8.

Figure 6 .
Figure 6.The confluence of Bivachny Glacier flowing from upper right into Fedchenko Glacier flowing from the left towards the observer in August 2015 (Photo: A. Lambrecht).

Figure 6 .
Figure 6.The confluence of Bivachny Glacier flowing from upper right into Fedchenko Glacier flowing from the left towards the observer in August 2015 (Photo: A. Lambrecht).

Figure 7 .
Figure 7.Estimated ice volume flux along the central flowline shown in Figure 1.

Figure 7 .
Figure 7.Estimated ice volume flux along the central flowline shown in Figure 1.

Table 2 .
Optical and synthetic aperture radar (SAR) data used for velocity determination by feature tracking.As a robust measure of their respective accuracy, the 90th percentile on bedrock adjacent to Bivachny Glacier is listed.