Glacial Research in the Karakoram


Introduction to Regional Climate Change, Surging Glaciers and Glacial Lake Outburst Floods in the Central Karakoram Mountains of Pakistan

Regional climate change directly affects glaciers and their hydrological systems (Parry et al., 2007). Besides the reduction in mass, glaciers also exhibit a broad range of responses to climate change. The speed and magnitude of the response to climate sensitivity of a particular glacier depends, among other factors, on its geometry, and can have significant variations (Banerjee and Shankar, 2013). Therefore, it is considerably challenging to relate climate change to the variations in the response to climate of individual glaciers. For example, some glaciers may undergo surging phases, which are periods of sudden advance and/or episodes of exceptionally high-speed flow, caused by factors that are not necessarily related to climate change (Meier and Post, 1969; Kamb et al., 1985, 1987; Kotlyakov et al., 2004, 2008). Given the rapid and severe changes the system undergoes during surges, it is important to observe both how the surface is affected and how the mass of the glacier changes during the quiescent (stagnation) period.

The Karakoram Mountains represent an important region to study these extreme processes especially since the return periods of Karakoram glacier surges are poorly quantified. Surging activity needs to be better understood if accurate mass balance assessments of Hindu-Kush–Karakoram–Himalaya glaciers are to be made. The distinctive surface morphology of Yukshin Gardan and Khurdopin glaciers in Shimshal Valley such as large, looped moraines is indicative of glaciers that have surged many times in the past, and reports of periodic ice dams blocking the Shimshal River (and associated outbursts from Virjerab Lake) have historically coincided with surge events, some of which dating back to the late 1800s (Visser, 1926; Todd, 1930; Iturrizaga, 2005).

As well as recording and interpreting the changes in the glacier, the presence of surges has hazard-related implications such as GLOFs due to the formation of ice dams and thus, the study of this behaviour is of paramount importance, particularly given the many lives that are under threat (Kääb et al., 2005). Glacial-dammed lakes play a prominent role in the natural landscape shaping processes of the Karakoram. At present, about 30 potential glacier dams are recorded; over half of them are situated in the western and central Karakoram. Shimshal valley is perhaps the most notorious case reflecting these catastrophic GLOFs.

Furthermore, glaciers in the Karakoram-Himalaya play a significant role in providing water to over 1.3 billion people in Asia. Determining and understanding the changes affecting glaciers and the resulting water in these regions is crucial for the sustainable development of local populations.

          The Glacial Situation in the Shimshal Valley

The Shimshal valley at the Karakoram North side, 60 km in length, is one of the few unglacierized trunk valleys. The parallel-running Hispar valley is taken up by a 62 km long glacier, as well as the neighbouring Batura valley (58 km). But instead, the N-exposed tributary glaciers Khurdopin/Yukshin Gardan (47/17 km), Yazghil (31 km) and Malungutti (23 km) advance into the Shimshal valley and form potential glacier dams over a horizontal distance of merely 25 km.

The three glacier dams Khurdopin/Yukshin Gardan, Yazghil and Malungutti are spread over the Shimshal valley floor, which is up to 2 km broad. In the last century, on average every five years, glacier lake floods occurred in the upper Shimshal valley. The valley flanks in the Shimshal valley are covered with thick slope moraines up to several hundred meters above the valley floor, which are reworked into transglacial debris accumulations since the deglaciation (Iturrizaga, 1999, 2002b). During glacier lake outbursts the unconsolidated debris accumulations are eroded and tremendous amounts of debris is transported downstream. These hyperconcentrated glacial mudflows possess a high-erosion potential. The sediment yield of the Hunza River amounts 4.800 t km2/yr (Ferguson, 1984: pp. 581) and also the Shimshal River is fed by suspension-rich tributary rivers. Peak discharge amounts of glacier lake outbursts in the East-Karakoram reached up to 7100 m3/s (Hewitt, 1982: 267 Feng and Qinghua, 1991).

The Glacial Situation in Shimshal valley. Copyrights: Itturizaga, 2001.


In order to assess the state of health of glaciers in the Central Karakorum required two sets of data: i) glacier terminus position and ii) glacier velocity. To assess the risk of GLOFs the KAP science team comprising Sergiu Jiduc and Oliver Forster, gathered one set of data: iii) field based geomorphological mapping of glacigenic landforms.

Glacier terminus position or the physical position of the snout of the glacier was estimated using repeat photography. A series of historic photographs taken by the early 20th century European Expeditions across the Karakoram was acquired from the Royal Geographical Society with IBG in London. Based on landmarks in the photograph, the position of the terminus was approximated and sketched on a map. Photo pairs provided visual representations of glacier change over time and mapped terminus positions were used to estimate changes in glacier length. Repeat photography methodology included: photo site selection and establishment, photography, and terminus plotting. The photographs will be later georeferenced and the terminus position precisely digitised using Geographical Information System (GIS).

Glacier surface velocity was measured using a GPS Real Time Kinematic (RTK) Survey (more details about this technique in Section). Velocity stakes were installed on the glacier surface and their position was repeatedly surveyed using high precision Trimble R10 GNSS receivers. Methodology included: velocity stake installation, data collection, post-processing, and velocity calculation (Karpilo, 2009). The science team looked for till deformation processes and estimated the glacier’s profile in order to determine whether the glacier is in quiescence or surging state. The velocity data collected in field is currently analysed at Imperial College London.

Geomorphological mapping of the glacier terminus provided information about: i) the spatial distribution of glacigenic landforms; ii) their genesis; iii) the distribution of periglacial process; iv) the extent of supraglacial debris cover; v) the meltwater drainage mechanisms; and vi) lake outburst flood probability. We particularly looked for surging evidence such as: concertina or zig-zag eskers, thrust moraine complexes, crevasse-squeezed ridges, multiple stacked diamictons, tectonised sediments and hummocky moraine complexes. The mapping procedure followed the guidelines explained by Hubbard and Gasser, (2005) and Otto and Smith, (2011) and was divided into three steps: 1) pre-field preparatory activities; 2) actual fieldwork; and 3) post-mapping activities. The identification and mapping of glacigenic landforms will use the criteria developed by Hambrey (1994), Bennett and Gasser (1996), Benn and Evans (1998), Kirbride et al. (2001) and Lukas (2002).

The likelihood of occurrence of a GLOF was estimated using data gathered by the geomorphic mapping procedure. By using an empirical model developed by McKillop and Clague, which incorporates parameters such as: lake surface area, dam height, lake volume, moraine width, moraine height, elevation, slope of moraine fan and lithology of moraine, we were able to generate a quantitative probability that a GLOF will occur, as well as estimations of outburst peak discharge, maximum volume of material, travel distance and maximum area of inundation (McKillop et al. 2007).

 Specialist Equipment and Technique

In order to calculate the surface ice flow velocity of Yukshin Gardan glacier, the KAP research team used an innovative differential GNSS survey technique called Real Time Kinematic (RTK). The RTK technique is based on the use of carrier measurements and the transmission of corrections from the GPS base station, whose location is well known, to the rover, so that the main errors that drives the stand-alone positioning cancel out.

GNSS Measurements

 Planning the GNSS surveys

  • GNSS survey goal:

Main goal: To calculate the relative surface velocity of i) Yukshin glacier (terminus + Equilibrium Line Altitude areas) and ii) Khurdopin glacier (terminus area)

Supporting objective: To construct an array of 40+ points across Yukshin and Khurdopin glaciers and measure their precise positions over a predetermined time interval using an accurate real time GNSS technique (RTK).

  • GNSS survey plan:

Planned to carry out RTK measurements by establishing an on-site base station, which can transmit real time differential corrections to the rover through the radio module. The base has been centred on a solid point on static ground (moraine) and its coordinates have been precisely logged (on site reference system). The same coordinates have been used as the reference for all the surveys. The key objective was to see how a point on the glacier moves throughout time at RTK precision.

  • Schedule:

4 days drilling to establish the velocity points followed by 4 additional days marking boulders, 20 days survey time interval, with 2 surveys each requiring 4 days to perform. In total we have worked 16 days on the GPS survey. The survey lasted for 36 days including all preparatory activities and waiting time.

  • Contingency plan for the fieldwork:

Use tagged boulders instead of velocity stakes if they fail or too difficult to carry out; Reduce coverage area

  • Check all equipment prior to mobilisation to site:

Due to the many constituent parts of GNSS receiver systems, especially RTK systems: batteries, antennae, receivers, radios, data cards or loggers, associated cables and sundries, it was good practice to check all equipment prior to mobilisation.

  • Order/check fundamental control station data:

The coordinates for the base control station were the same for every survey. Backed up data.

  • Geodetic control

Geodetic control has been established on day one at a fixed location near our base camp (3678m) on the top of the 100m high lateral moraine. The on site reference base station had clear views across the glacier for a distance of over 20km. The base station location and its coordinates have been determined in VGS system and precisely registered in the field notebook using an ordinary GPS as well as the Trimble R10 receiver. The base station had remained centered on the fixed point (painted dot) throughout the length of the survey (30 days). The same reference point has been used for both surveys.

Base Point name: 1000

Location: Lat: 36°19’37.01400″N / Long: 75°26’27.62734″E

Elevation: 3678.039m

Coordinate system & Projection: UTM (Universal Transverse Mercator)

Zone: 43 North

Datum: WGS 1984

Antenna height: 2.039

  •  Survey Style & Method
  1. Equipment used: 2x Trimble R10 (one rover and one base) and 1 x TSC3 controller
  2. High-Precision Real-Time kinematic (RTK)
  3. On the Fly (initialisation type)

RTK does not require an initialisation on a known baseline. Thus, the survey can tolerate periodic loss of lock during the survey, as the ‘integer ambiguity’ can be determined whilst the rover is moving to the next point. In such surveys, as the receiver at which the baseline solution is being calculated is moving, a short epoch setting is needed. Typically an epoch setting of one or two seconds is used for standard detailing for land surveys. In the on-the-fly technique, the most advanced baseline solution methods are used to determine the baseline initialisation and then to test it, to ensure it is the correct one when initialising on-the-fly surveys it is important to move to a location, which has an open view of the sky. Once the initial baseline is computed, the coordinates of the roving receiver are computed for each epoch. The integer ambiguity of the baseline is known from the initialisation process and thus the new baseline solution at each epoch can simply be determined from the change in the satellite observations. 43 points GPS have been measured on the glacier surface of which 3 have been compromised

  • Precision and quality control on the GNSS observations

We have used the most advanced GNSS commercial solution with a precision of 8mm +0.5 ppm horizontally and 15mm +0.5 ppm vertically. Carrier phase measurements are therefore extremely precise (down to fractions of a millimetre), but they contain an unknown integer initialisation constant, the so-called “phase ambiguity”. RTK positioning resolves the integer ambiguities to achieve the high level of precision (called initialisation). Convergence time is needed to fix the phase ambiguities. This time depends on the processing algorithm and the distance between base and rover, and ranges from a few seconds to a few minutes. RTK technique is based on the following high-level principles:

  1. In clear-sky locations, the main errors in the GNSS signal processing are constant, and hence they cancel out when differential processing is used. This includes the error in the satellite clock bias, the satellite orbital error, the ionospheric delay and the tropospheric delay.
  2. The noise of carrier measurements is much smaller than that of the pseudo-code measurements. However, the processing of carrier measurements is subject to the so-called carrier phase ambiguity, an unknown integer number of times the carrier wave length, that needs to be fixed in order to rebuild full range measurements from carrier ones.
  3. The phase ambiguity can be fixed for dual-frequency differential measurements for two close receivers.
  • Quality Control

The survey measurements, field records, sketches, and other documentation were examined to verify compliance with the specifications for the intended accuracy of the survey (i.e. 1cm horizontal, 2cm vertical). This examination may lead to a modification of the intended accuracy. The survey accuracy was checked by comparing minimally constrained adjustment results against established control. This comparison takes into account the network accuracy of the existing control, as well as systematic effects such as crustal motion or datum distortion. Accuracy is the measure of the difference between a particular measured coordinate and its true value, often quoted as the root mean square error.

  • Key Quality Control: Occupied a control point with known coordinates. This ensured that the displayed co-ordinates are correct.
  •  Executing Survey Measurements 
  1. GNSS rover antenna was mounted on top of a 2m pole, the operator used a hand-held data logger (controller)
  2. Operator performed 2 -3 observed control point measurements per survey, per point each, consisting of 360 Epochs (6min) per point. We have carried out 2 surveys in total.
  3. Surveys have been carried out from higher to lower elevations (e.g. 4000m to 3500m) on high stable landforms such as ice pinnacles, debris mounds as well as in depressions such as supraglacial valleys.
  4. A higher resolution of points was used at the confluence of Yukshin – Khurdopin in order to assess the ice flow interaction between Yukshin and Khurdopin glaciers.
  5. Logged each measurement in a field notebook such as name, location, elevation, date, time, author, comments
  6. Used the same coordinate values when performed the two surveys
  7. Did not mix receiver and antenna types
  8. Re-measured the base station antenna height at the end of the survey.
  9. Download data daily
  10. Recharged all batteries overnight.
  • Processing and Results
  1. We have gathered the GPS data using Trimble Access Software. Post processing of data has been carried using Trimble Business Centre, Excel, & AutoCAD software. Processing activities included exportation of dxf, csv files, conversion and calculation of displacement (DGPS) using AutoCAD at GISCAD Trimble Dealer, Romania. Our working precision was less than 1 cm.
  2. We have used the manufacturer’s defaults for processing parameters.
  3. GIS, processing will follow shortly with the creation of a DEM, glacier polygons and velocity vectors. In addition, we will use remote sensing to extrapolate the vectors across the glacier.

      Preliminary Results

  • Khurdopin glacier is moving relatively slowly or at normal quiescence velocities of approx. 1 to 5cm/day or 3 to 20 m/year); Yukshin glacier on the other hand is moving faster with recorded velocity values of 5 to 35cm/day or 20 to – 130m/year depending on location. It is worth mentioning the Yukshin glacier is an avalanche fed type glacier situated in a steeper sloping valley, hence the higher velocities.
Surface Velocity across Yukshin glacier.
Surface Velocity across Yukshin glacier.
  • However, measured surface displacements of Khurdopin during surging reached > 5 km/year. Therefore we can conclude, that the glacier is not surging and it is in a quiescent phase. Therefore the people of Shimshal are safe for now! From the topographic relief point of view, the glacier is changing every year due to continuous ice degradation. The surface of Yukshin and Khurdopin glaciers showed an increase in the number of thermokarst features such as rapidly size and form-changing depressions and supraglacial ponds. The surge surface topography became less chaotic and smoother as compared with previous observations. However, the issue of debris cover may complicate the story.
  • Beyond the glacier response to changing climatic conditions, glacier topographic factors (e.g. aspect, altitude, slope, morphology) will influence glacial dynamics. In the case of debris-covered glaciers, the influence of the debris cover must be deeply considered. Previous studies revealed that glaciers with extensive debris cover have a qualitative difference in their response to a warming climate as compared with bare ice glaciers (Banerjee and Shankar, 2013). Depending on its thickness, the debris cover might enhance or alternatively hamper glacier ablation (Mihalcea 2006; Benn and Evans 2010). At similar latitudes and elevation, but in the Northern Hemisphere, Mihalcea et al. (2006) found ablation rates that varied between 3 to 6 cm/day with a debris layer of 18 cm to 0 cm thickness for the debris-covered Baltoro glacier in the Karakoram (4178m, 35° N) with GNSS surveys. In the Bagrot Valley (Karakoram, 36° N), the lowest ablation rates were ~2.3 cm/day, with a debris cover thickness averaging 37.5 cm, between 2500 m to 3350 m (Mayer et al., 2010).
Passu Glacier in 1887 (right) and 2015 (left).
Passu Glacier in 1887 (right) and 2015 (left). Tim Taylor, Photography

       Preliminary Interpretations

  • At present there is no significant risk of GLOF. The Shimshal River drains trans-glacially along the length of its contact with Khurdopin glacier. There are small patches of the sub-glacial drainage but they are limited. Khurdopin-Yukshin glacier would need to advance at least 20m in order to damn the river. Glaciers in the region are retreating overall as our repeat photos demonstrate (see above image pair showing Passu glacier). Moreover, everyone we have spoken to in Pakistan has also indicated that in the last few years the glaciers have gone backwards and their levels have gone down suggesting recession and thinning.

         Preliminary Recommendations

  • The next five years are critical in region because Khurdopin surge periodicity is believed to be approximately 20yrs, with the last surge occurring in the late 1990’s.
  • The influence of Yukshin Gardan glacier remains uncertain. It could have a role in creating a GLOF (more research is needed)
  • Continuous monitoring of the two glaciers over the next 5 years is paramount
  • Installation of a water level, a permanent GPS and a weather monitoring station on Khurdopin and Yukshin glaciers could provide valuable information about ice mass balance fluctuations and weather patterns, which in turn can aid mitigation and adaptation strategies and government action in the region.
  • Close monitoring of the area is recommended during critical times (e.g. summer & surging phase) in order to inform an early warning system