Corona is high-resolution (1.8 m) imagery acquired between 1960 and 1972 under the US Key Hole program, which was subsequently declassified in 1995 (Dashora et al., 2007). Corona imagery provides an invaluable dataset for glaciologists, but contains severe distortions that must be removed during the processing. Two Corona stereopairs were available that covered parts of the Manaslu region. The first stereopair from 1967 covered Ponkar Glacier and was partially affected by cloud. The Corona stereopair covering Ponkar Glacier were therefore only used for calculating the changes in glacier area and not geodetic glacier mass balance. The second Corona stereopair covered Thulagi, Punggeon and Hinang Glaciers and was from 1970. This image was cloud free and contained sufficient image contrast to generate a DEM.

As there is no mathematical model for Corona stereo-imagery, in order to generate a DEM the data must be processed in a non-metric approach which requires many tie-points in order to solve the interior and relative orientations (Schmidt et al., 2001; Galiatsatos et al., 2007). All the Corona imagery was processed in ERDAS Imagine 2015. We solved the relative orientation using many (mostly automatically generated) tie points, while the exterior orientation was solved using Ground Control Points (GCPs). The Ponkar Corona image was processed with 28 GCPs and 359 tie points while the Thulagi Corona image was processed with 75 GCPs and 1837 tie points. We used a combination of a 5 m RapidEye image from 2012 and the SRTM, as well as Google Earth imagery (after Surazakov and Aizen, 2010) as a source of GCP information, which although not ideal, allowed the products to be georeferenced enough for the DEMs to be later co-registered with the other DEMs (section DEM Co-registration). This allowed the images to be triangulated with an RMSE of <5 pixels. Orthomosaics from both datasets were exported at 2 m resolution, while the DEM from the Thulagi Corona image was exported at 8 m resolution using the eATE (enhanced Automatic Terrain Extraction) tool within ERDAS IMAGINE. All the DEMs were projected to UTM zone 45 N and geoid/ellipsoid differences were accounted for.

A combination of correlation masks and manual interpretation using the hillshade models was used to clean the DEMs and remove pixels with poor correlation scores (<0.6) or noticeable interpolation errors. This had the undesired effect of leaving large data voids over the accumulation areas of some of the glaciers. The DEMs were then resampled to the spatial resolution of the coarsest DEM (30 m). This was done using a median block statistics operation to avoid biases related to resolution variations of the DEMs and thus in their ability to depict terrain curvature and roughness (Paul, 2008; Gardelle et al., 2013).

Landsat images from 2013, 2005, and 2001 as well as Corona data from 1967/1970 (henceforth referred to as 1970) were used to assess the glacier area change (Table 1). The 2013 outlines were created by Robson et al. (2015) and were created semi-automatically using a combination of Landsat 8 imagery, a DEM and SAR Coherence data within an Object Based Image Analysis (OBIA) classification with additional manual corrections applied. The 2013 outlines were manually modified in order to create glacier outlines from 2005, 2001, and 1970. Noticeable breaks in slope morphology, visible vegetation change, exposed ice, or collapses in terminal moraines were used as indicators of underlying ice when mapping the debris-covered ice. The boundary between the clean ice and debris-covered ice was delineated for each set of outlines in order to monitor changes in the debris-covered area.

Landsat images from 2013/2014, 1999/2000, and 1993/1994 (Table 1) were used for determining glacier velocities. The lack of image contrast on clean ice meant we restricted our velocity analysis to the debris-covered tongues of the glaciers, from the glacier terminus, to the transition between the clean ice and debris-covered ice. Prominent features on the debris-covered glaciers were tracked between sets of successive satellite images in order to determine the surface velocity. Features were matched using cross-correlation on orientation images (CCF-O) using the free software CIAS, developed by the University of Oslo (Kääb and Vollmer, 2000). Orientation images use the gradients between neighboring pixel values instead of the raw digital numbers, thereby reducing the influence of differences in scene illumination. Heid and Kääb (2012a) demonstrated that cross-correlating orientation images outperforms other image matching methods on debris-covered glaciers.

The image pairs were referenced together within CIAS using a Helmert Transformation with 15 to 20 points on stable terrain. Landsat 5 and 7 images were registered with RMSEs typically of 10–12 m, while Landsat 8 images were registered with RMSEs typically of less than 5 m.

Displacements were filtered by direction and magnitude. Displacements with a signal to noise ratio (SnR) lower than 0.5 were excluded. A 3 × 3 moving window was then used to remove displacement vectors that varied by more than 20% in magnitude or direction to the mean values. Some manual filtering was needed for removing matches associated with snow, shadow, or mass movements from the valley sides. Transverse profiles were plotted at the boundary between the clean ice and debris-covered ice in order to examine the changing debris-covered ice-flux. Additionally, glacier centrelines were generated automatically using the method of Kienholz et al. (2014) based on outlines from the Randolf Glacier Inventory (RGI) version 5.0 (Pfeffer et al., 2014) and the void-filled SRTM. The centrelines were manually edited to adapt them to changing glacier size in the older satellite imagery and the mean velocity extracted along the centrelines.

The hypsometry of the 41 glaciers was computed for 2000 and 2013 using the SETSM and SRTM DEMs. The glacier outlines from 2013 and 1999 were divided into 100 m elevation bins and the area of each was determined in order to create hypsometric curves. A hypsometric index (HI) based on the maximum (H_max), median (H_med), and minimum (H_min), elevations of each glacier was calculated to characterize the glacier hypsometry based on the method of Jiskoot et al. (2009) and King et al. (2017):

Glaciers were characterized by their HI as very top heavy (HI < −1.5), top heavy (−1.5 < HI < −1.2), equidimensional (−1.2 < HI < 1.2), bottom heavy (1.2 < HI < 1.5), or very bottom heavy (HI > 1.5).

The total glacier area of our study area in 2013 was 470.4 km² of which 84 km² (17.8%) was debris-covered ice. Between 2001 and 2013 the glacier area reduced by 42 km² (8.2%). Additionally, over half (52%) of the individual glaciers underwent an increase in debris-cover area over this time-period, with two glaciers (G084354E286669N and G084689E28511N) gaining over 1.5 km² of debris-cover, increases of 22 and 12% respectively.

The largest spatial coverage of elevation changes was obtained between 2000 and 2013 by comparing the SRTM and SETSM DEMs. The glaciers underwent a mean annual surface elevation change of −0.25 ± 0.08 m a⁻¹. Over the whole study area, the geodetic mass balance was found to be negative (−0.21 ± 0.16 m w.e. a⁻¹) with individual glacier mass balances varying from −0.79 ± 0.12 to + 0.30 ± 0.17 m w.e. a⁻¹. There are clear spatial patterns to the glacier change (Figure 3). Generally, the glaciers in the north (both those on the Tibetan Plateau, and those flowing south from Himlung) had less negative mass balances, while the glaciers located further to the south, especially those around Manaslu had more elevation differences.

Glaciers south of the topographic divide are primarily debris-covered below 4600 m. Between 1999 and 2013 these glaciers lost 21.7 km² of ice (6.7% reduction in area). Many of the glacier termini remained in the same position, with the majority of glacier area changes occurring at the glacier margins and at the boundary between the clean ice and debris-covered ice. Some of the glaciers however underwent large terminus retreats. Thulagi glacier retreated 1.37 km between 1970 and 2013. Additionally, several large eastern flowing debris-covered glaciers (Punggeon Glacier: G084595E28536N and Hinang Glacier G084640E28489N : G084640E28489N) underwent large frontal retreats. Some debris-covered glaciers did not undergo sizable changes in area (for example G084539E28644N and G084550E28613N). Individual glacier area changes span from −9.0 to +1.1%.

Between 2000 and 2013, the southern glaciers in our study region had an average surface elevation change of −0.32 ± 0.08 m a⁻¹ which equates to a mean geodetic balance of −0.27 ± 0.09 m w.e. a⁻¹. There is a linear relationship (R² = 0.67) between elevation and the rates of glacier thinning over the entire glacier. The effects of the debris-cover appear to be complicated. Between ~4700 and 5700 m a.s.l rates of elevation change of debris-covered ice are greater than those of the clean ice components (Figure 4). The thinning rate generally increases as elevation decreases, reaching 1.56 ± 0.08 m a⁻¹ at ~4000 m a.s.l. Below 4000 m a.s.l. the thinning rate decreases to 0.78 ± 0.08 m a⁻¹ at ~3600 m a.s.l. before increasing to 1.10 ± 0.08 m a⁻¹ at 3300 m a.s.l.

The northern glaciers on the Tibetan Plateau are predominantly clean ice glaciers and have a higher mean elevation (5991 m a.s.l.) than the southern glaciers (5558 m a.s.l.) (Table 3). The glaciers lost an average of 6.7% of their area between 2000 and 2013. There appears to be an east-western trend in area change with G084132E28911NE and G084132E28911NW, in the far northwest of the study area, shrinking by 3.1 km² (13.5%) and 0.9 km² (16.5%). On the other hand G084424E28839N and G084461E28797N, two glaciers in the northeast of the study area and at approximately the same elevation, lost 0.2 km² (2.9%) and 0.3 km² (1.8%) of their area, respectively.

Despite losing 6.7% of their area, the northern glaciers were approximately in balance with a mean geodetic mass balance of −0.03 ± 0.11 m w.e. a⁻¹ and a mean surface elevation change of −0.04 ± 0.08 m a⁻¹ (Figures 3, 5). Individual geodetic mass balances ranged from −0.25 ± 0.15 to +0.30 ± 0.17 m w.e. a⁻¹. Some glaciers thickened over the majority of their surface, however these glaciers are situated at very high altitudes (minimum elevations > 6500 m a.s.l.). Thinning rates are inversely proportional to elevation with an R² of 0.87 (Figure 4). Ice situated at ~ 5400 m a.s.l. thinned by 1.66 ± 0.08 m a⁻¹ and ice at ~ 6700 m a.s.l. thickening by 0.99 ± 0.08 m a⁻¹.

Our results show changing rates of both glacier area as well as glacier mass. Eight of the glaciers (~200 km²) were covered by Corona imagery from 1970. Between 1970 and 2013 these glaciers lost 37.9 km² of ice (17.9% area reduction), equating to a mean rate of −0.36 km² a⁻¹ between 1970 and 2001, −0.97 km² a⁻¹ between 2001 and 2005, and −2.28 km² a⁻¹ between 2005 and 2013. Over the entire study region, between 2001 and 2005 the total glacier area reduced at a rate of 6.4 km² a⁻¹ equating to 1.2% a⁻¹ while less glacier area per year was lost between 2005 and 2013 (2.06 km² a⁻¹ equating to 0.4% a⁻¹).

Twenty-one glaciers were covered by the 2013, 2005, and 2000 DEMs (Figure 1). The subset of these glaciers had a more negative mass balance (−0.44 ± 0.09 w.e. a⁻¹) than the regional average (−0.21 ± 0.16 w.e. a⁻¹). Considering that the subset excluded the more in-balance glaciers on the Tibetan Plateau this is not surprising. The highest regional glacier mass losses were found between 2000 and 2005 (−0.83 ± 0.74 m w.e. a⁻¹) although this varied through the region, with the three large glaciers flowing out of Manaslu losing most mass between 2005 and 2013 (−1.20 ± 0.87 m w.e. a⁻¹). The glaciers that were covered by the Corona data show changing rates of glacier loss. The rates of glacier downwasting have increased over the 43-year period, from −0.46 ± 0.03 m a⁻¹ between 1970 and 2000 to −0.99 ± 0.08 m a⁻¹ between 2000 and 2013. The geodetic mass balance for the three glaciers between 1970 and 2013 was −0.31 ± 0.12 m w.e. a⁻¹ (Figure 6) and −0.84 ± 0.11 m w.e. a⁻¹ between 2000 and 2013 (Table S2).

In 2013 the southern glaciers had a mean HI of −1.29 (Table S1) indicating top heavy hypsometries with more area above the median glacier elevation. The northern glaciers on the other hand had a mean HI of 1.03 relating to a bottom heavy hypsometry. Figure 7 demonstrates the difference in hypsometries of northern and southern glaciers. Between 2000 and 2013 the southern glaciers have become more top heavy (HI change from −1.16 to −1.29) while the northern glaciers have become more bottom heavy (−1.14 to 1.03).

From the 2013 hypsometry, 70% of the total glacier area of the entire study region exists between 5200 and 6500 m a.s.l. This is also where the ELA approximated by the median elevation lies (~5748 m a.s.l.) meaning that a small rise in the ELA will expose a significant amount of ice to melting. Similar findings were also discussed by Shea et al. (2015) and King et al. (2017), who foresee a possibility of the debris-covered tongues becoming detached from the accumulation areas resulting in a bi-model hypsometry. In contrast, the bottom-heavy hypsometry of the northern glaciers implies that although a rise in ELA would expose ice to melting, future losses are likely to be limited by smaller areas of ice stored at higher elevations.

There is no statistical relationship between the glacier hypsometric index and glacier change, although some statements can be made. Punggeon Glacier (G084584E28535N) for example is very bottom heavy (HI = 1.81). The glacier has amongst the highest rates of surface elevation change in the region (−0.74 ± 0.08 m a⁻¹). Some of the glaciers storing large amounts of ice at higher elevations also have faster velocities, for example Ponkar Glacier: G084459E28695N; HI: −2.15 and Thulagi Glacier: G084549E28515N; HI: −1.25 velocities of 47.0 ± and 40.6 m a⁻¹ respectively.

There is a large variation in glacier speed between individual glaciers (Figures 8, 9; Table S3), Ponkar Glacier (ID: G084456E28737N) and Thulagi Glacier (G084595E28536N) flowed at mean velocities in 2013/2014 of 47 ± 8.2 and 41 ± 8.2 m a⁻¹, respectively, while Baudha Himal Glacier (G084684E28407N) and Northern Fukang Glacier (G084425E28837N) have mostly stagnant debris-covered tongues with mean velocities of <3 m a⁻¹.

There is no clear signal in glacier velocity changes in the region, with acceleration, deceleration, and constant velocities observed over the study period. Kechakyu Khola Glacier (G084425E28700N), Suti Glacier (G084354E28669NW), and Himal Chuli Glacier (G084759E28427N) have slowed by 12 m a⁻¹ (42.9%), 15.6 m a⁻¹ (29.0%), and 26.7 m a⁻¹ (45.7%) respectively between 1993/1994 and 2013/2014. Other glaciers have decelerated by a greater proportion, however these glaciers are flowing at such slow velocities (<3 m a⁻¹) that these changes in velocity are of a small magnitude. The termini of some glaciers such as Hinang Glacier (G084640E28489N) have noticeably stagnated (Figure 10). G084759E28427N has stagnated to the point of having near-zero velocities in 2013 (Figure 6D). Thulagi and Ponkar glaciers appear to show alternating patterns of glacier velocity, with slight accelerations in glacier velocity between 1993/1994 and 2013/2014 (5 and 6%, respectively). However, if the velocities are examined only between 1999/2000 and 2013/2014 then Ponkar Glacier decelerated by 10.4% while Thulagi Glacier accelerated by 27.3%. Regionally, the mean velocity varies only slightly between 1993/1994 and 1999/2000 (34.9 ± 5.6 and 32.4 ± 3.2 m a⁻¹ respectively) before decreasing by a mean of 5.2 m a⁻¹ (16.0%) between 1999/2000 and 2013/2014.

There is a noticeable spatial variability of glacier mass change between the northern and southern glaciers, as well as within the northern and southern regions. We now explore some of the potential influencing factors of glacier mass balance within the Manaslu region.

Dona Lake expanded from 0.50 km² in 1970 to 0.93 km² in 2013. During this time the glacier downwasted at a mean rate of 0.49 ± 0.03 m a⁻¹, the highest rates in the Manaslu region. Between 1970 and 2001 Thulagi Glacier terminated in Dona Lake, however by 2005 the glacier appears to have retreated onto land. While moraine subsidence between 1996 and 2009 occurred at a rate of ~1 m a⁻¹, lake surface elevations lowered by 0.3 to −0.5 m a⁻¹ between 1996 and 2009 (Mool, 2011). This implies that the lake discharge through the old moraine is greater than the water added to the lake from the melting of the glacier. With the glacier now terminating on land it seems unlikely that lake will continue to expand, even with increasing glacier melt rates.

In addition to losing mass, Thulagi glacier exhibits complex patterns in velocity. The glacier velocity decelerated between 1993 and 1999 by 16.7%, followed by an acceleration in velocity of 29.7% between 1999 and 2013. A similar trend is seen in the glacier hypsometry. Between 1970 (HI = −1.18) and 2000 (HI = −1.37) Thulagi glacier become more top heavy, while between 2000 and 2013 Thulagi glacier become less top heavy (HI = −1.25).

We calculated the DEM differencing uncertainties in terms of stochastic error and systematic error. The stochastic error was based on the standard error and the mean elevation difference on stable terrain. Magnusson et al. (2016) showed how standard methods for quantifying DEM differencing uncertainties neglect the spatial dependence of DEM errors, and as such typically overestimate uncertainties. It should therefore be noted that our uncertainty analysis is conservative.

Even with our conservative error estimation, 64% of our surface elevation change results are significantly different from zero, as well as eight of the 10 regional averages. Fifty-eight percent of our geodetic mass balance results however were significant, with only the regional average from 1970 to 2013 being larger than the associated error bars. As stated by Gardelle et al. (2013), the largest component of the uncertainty is the estimation of the SRTM penetration. Values range from 1 to 10 m (Gardelle et al., 2013; Kääb et al., 2015; Berthier et al., 2016) which is considerably greater than the stochastic error. Additionally sources of error include the assumed constant density of 850 kg m⁻³ across the entire glacier (Huss, 2013).

In the absence of a more rigorous method, we estimated uncertainties for our velocity results based on ± one standard deviation. This assumes that no displacements have occurred on non-glacier terrain. Rockfalls and movements are common in the High Himalayas meaning our uncertainty estimates are most likely overestimated. Nevertheless, 94% of our velocity measurements are significant, in addition to all three regional averages.

Our findings add to the consensus of strong yet spatially variable glacier mass losses across the Himalayas. In recent years, studies looking at geodetic mass balances of Himalayan glaciers have become more common (Gardelle et al., 2013; Pellicciotti et al., 2015; Bhattacharya et al., 2016; Ragettli et al., 2016; Brun et al., 2017). Due to the differing climatic regimes across the Himalayas, we restrict our comparisons to the Central Himalayas and the Everest Region. Our results are comparable to surface elevation change rates in the Everest Region of −0.32 ± 0.08 m a⁻¹ between 1970 and 2007 (Bolch et al., 2011) and −0.28 ± 0.42 m a⁻¹ between 1974 and 2006 in the Langtang Himal (Ragettli et al., 2016). This is very comparable to our result of −0.36 ± 0.03 m a⁻¹ between 1970 and 2013. Both our results and the findings of Bolch et al. (2012) indicate increases in the rate of mass loss between 2000 and 2013 relative to other time periods.

When looking at shorter time periods, Kääb et al. (2012) found a mean annual surface elevation change of −0.21 ± 0.05 m a⁻¹ over the entire Hindu Kush Himalaya Region between 2003 and 2008, with some regions losing as much as −0.66 m a⁻¹. Brun et al. (2017) found mean annual geodetic mass balances between 2000 and 2016 of −0.33 ± 0.20 and −0.34 ± 0.09 m w.e. a⁻¹ for East and West Nepal, respectively. The losses we found between 2000 and 2013 of −0.25 ± 0.08 m a⁻¹ therefore compares well with these regional averages. Our results of geodetic mass balance for the southern glaciers between 2000 and 2013 (−0.27 ± 0.09 m w.e. a⁻¹) are less negative than the results found by Gardelle et al. (2013) for western Nepal (−0.32 ± 0.13 m w.e.a⁻¹) but is in line with the result from the Everest region (−0.26 ± 0.13 m w.e.a⁻¹) over roughly the same time period. Bolch et al. (2011) found a more negative mass balance of −0.79 ± 0.52 m w.e.a⁻¹ between 2002 and 2009 in the Everest Region. Our geodetic mass balance between 2005 and 2013 (−0.74 ± 0.55 m w.e.a⁻¹) is more comparable to their findings.

Our estimate of the ELA based on the median glacier elevation of 5748 m is in broad agreement with the estimate of ~5600 m for the ELA in the central Himalayas reported by (Bolch et al., 2012) and for Western Nepal of 5590 ± 138 m by Gardelle et al. (2013).

Quincey et al. (2009) found that the debris-covered glaciers in the Everest region flowed at ~20 m a⁻¹, which is slower than our mean result of 27 ± 8.2 m a⁻¹ in 2013. Kääb (2005) found that the debris-covered glaciers in Bhutan were flowing at higher rates of between 20 and 50 m a⁻¹. There is less work concerning changes in glacier velocity over time, and the work that does exist often find contradictory results of how glacier velocity has varied. Heid and Kääb (2012b) observed reductions in glacier velocity in areas with negative mass balances, including a 43% reduction in velocity in their nearest study site to the Manaslu region, the Pamir Mountains between 2001 and 2010. Haritashya et al. (2015) found reductions in glacier velocity of Khumbu glacier between 2003 and 2008 of between 10 and 20% while Bhattacharya et al. (2016) found a 6.7% decrease from 1996 to 2013 in the Indian Himalaya. The reduction we measured in this study of 15.0% between 1999 and 2013 is in line with these values.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.